使用全局变量后方程无法求解？？

chunshui2003

function uu=differential_equation(t,u)
%%%%%%%%%%%%%%%%%%    不平衡磁拉力参数
m1=1.5*10^4;               %发电机转子重量
m2=1.1*10^4;               %水轮机转轮重量
c1=6.0*10^4;               %发电机转子阻尼系数
c2=7.0*10^4;               %水轮机转轮阻尼系数
k1=8.5*10^7;               %上导轴承刚度系数
k2=6.5*10^7;               %下导轴承刚度系数
k3=3.5*10^7;               %水导轴承刚度系数
delta0=8*10^-3;            %发电机转子不偏心时平均气隙长度
e2=0.3*10^-3;              %水轮机转轮偏心距离
mu0=4*pi*10^-7;            %空气磁导系数
kj=5.2;                    %气隙基波磁动势系数
R=1.2;                     %发电机转子半径
L=0.5;                     %发电机转子长度
omega=10;                  %机组旋转频率   
K1=25*k1+9*k2+k3;
K2=-5*k1+3*k2+3*k3;
K3=k1+k2+9*k3;
% e1=0.5*10-3;           %发电机转子偏心距离
global e1
B1=e1*omega^2*cos(2*pi*omega*t)/delta0;
B2=e1*omega^2*sin(2*pi*omega*t)/delta0;

Ij=1500;
F0=R.*L.*pi.*kj^.2*Ij.^2.*mu0./(m1.*delta0^.3);         %%%%%%%    Ij没有定义呢

%%%%%%%%%%%%%%       线性密封力参数
alpha0=1.0;               % 动量修正系数
rho=1.0;                % 液体密度，单位kg/m^3
lambda=1.0;             % 沿程损失系数
xi=1.5;                 % 局部损失系数
Rm=2.925;               % 密封半径 单位m
lm=0.43;                % 密封长度 单位m
v=3.537;                % 液体轴向流速 单位m/s
cm=2.5*10^-3;           % 密封间隙 单位m
Kxx=lambda*lm.^2*rho*pi*v.^2*Rm*(1+xi)/(8*cm.^2*(1+xi+lambda*lm/2/cm));
Kyy=Kxx;
Kyx=-(Rm*omega/2)*(pi*rho*v*lm.^2/2/cm)*(1+xi+lambda*lm/6/cm);
Kxy=-Kyx;                    %%% 刚度系数
xishu=(2*(1+xi)+lambda*lm/4/cm)/(3*(1+xi)+lambda*lm/2/cm);
Dxx=(rho*v*lm.^2/2/cm)*pi*Rm*(1+xi+lambda*lm/6/cm);
Dyy=Dxx;
Dyx=-(Rm*omega/2)*(pi*rho*alpha0*lm.^3/cm)*xishu;
Dxy=-Dyx;                 %%%%% 阻尼系数
mxx=pi*Rm*rho*alpha0*lm.^3/cm*xishu;
myy=mxx;              %%%% 惯性系数

%%%%%%%%%%          线性密封力合成参数
M2x=m2+mxx;
M2y=m2+myy;
D1=(c2+Dxx)/M2x;
D2=Dxy/M2x;
D3=(c2+Dyy)/M2y;
D4=Dyx/M2y;
Kg1=Kxy/M2x;
Kg2=Kyx/M2y;
B3=m2*e2*omega.^2*cos(2*pi*omega*t)/(M2x*delta0);
B4=m2*e2*omega.^2*sin(2*pi*omega*t)/(M2y*delta0);
%%%%%%%%%%%        涉及到变量的一些参数
Z1=(1-sqrt(1-u(1)^2-u(3)^2))/sqrt(u(1)^2+u(3)^2);
Z2=sqrt(u(5)^2+u(7)^2)/sqrt(u(1)^2+u(3)^2);
Z3=1/Z2;
F0=R*L*pi*kj.^2*Ij.^2*mu0/(m1*delta0.^3);
F1=1/(1-u(1)^2-u(3)^2)*(Z1+Z1^3+Z1^5);
F=F0*F1;
%%%%%%%%%    转子-转轮系统运动微分方程
uu=zeros(8,1);
uu(1)=u(2);
uu(2)=B1+F*u(1)/sqrt(u(1)^2+u(3)^2)-c1/m1*u(2)-1/(16*m1)*(K1+K2*Z2)*u(1);
uu(3)=u(4);
uu(4)=B2+F*u(3)/sqrt(u(1)^2+u(3)^2)-c1/m1*u(4)-1/(16*m1)*(K1+K2*Z2)*u(3);
uu(5)=u(6);
uu(6)=B3-D1*u(6)-D2*u(8)-1/(16*M2x)*(K3+K2*Z3)*u(5)-Kg1*u(7);
uu(7)=u(8);
uu(8)=B4-D3*u(8)-D4*u(6)-1/(16*M2y)*(K3+K2*Z3)*u(7)-Kg2*u(5);



%%%%%%%%
clc
clear
tic
global e1
e1_sub=0.5*10^-3;
y0=[1 1 1 1 1 1 1 1];
period=1/10;         %采样的周期，转速的单位是10Hz,所以周期为“1/10”；如果转速的单位为“rad/s”，那么周期为“2*pi/omega”
    [t,u]=ode45(@differential_equation,[0:period/100:100*period],y0);




当函数里面的e1定义为全局变量时（如上），matlab一直提示“busy”；而当e1赋予值，不使用全局变量时，方程可以求解。
奇怪的是我拿另一个参数做实验“Ij”，无论是否将其定义为全局变量，方程均可求解。
我有点不懂这是为什么，难道是刚性方程？
因为涉及到参数变化的多次计算，不希望手动进行，太繁琐了，做循环计算是我想要的，但第一次都进行不了。请大家给我解答一下。

ChaChing

主程式中没看到有给e1赋值!?

chunshui2003

哦，我修改后忘记粘贴了。应该是
clc
clear
tic
global e1
e1_sub=0.5*10^-3;
e=e1_sub;
y0=[1 1 1 1 1 1 1 1];
period=1/10;         %采样的周期，转速的单位是10Hz,所以周期为“1/10”；
    [t,u]=ode45(@differential_equation,[0:period/100:100*period],y0);


本来是不需要使用“e1_sub”的，但是在实际操作时用的是循环，这里只是单一的一次计算。
现在没什么问题了，但还是提示“buzy”，无法运行。

chunshui2003

上面写错了，应该是
e1=e1_sub;

ChaChing

昨晚也试过这个, 没问题!

chunshui2003

请问是用我的方程计算的吗？如果是的话，那我想应该是我的matlab的问题，我再找其他人的机器用不同版本的试一下。谢谢你了！

ChaChing

非常抱歉!
刚刚再仔细看, 才发现昨晚是使用1F中的数据试的(% e1=0.5*10-3;%发电机转子偏心距离), 与0.5*10^-3;差很多!
试下发现的确如LZ说的一直busy, 原因不知!
但可以确定与使用或不使用全局变量无关! 是与e1的值相关

chunshui2003

ChaChing你好，谢谢你的建议。我的问题已经解决了！主要原因是在龙哥库塔法的初值选取上，经过查看计算过程，发现未知量的数值一直在变化，证明方程的系数传导没有问题。将初值修改后，如[2 2 2 2 2 2 2 2 ]，即可进行计算。估计是初值全部为1时，在某一个计算的环节不收敛，所以一直循环下去。希望今后大家遇到类似的问题将初始值更改一下也许问题就解决了！