使用全局变量后方程无法求解？？
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=;
period=1/10; %采样的周期，转速的单位是10Hz,所以周期为“1/10”；如果转速的单位为“rad/s”，那么周期为“2*pi/omega”
=ode45(@differential_equation,,y0);
当函数里面的e1定义为全局变量时（如上），matlab一直提示“busy”；而当e1赋予值，不使用全局变量时，方程可以求解。
奇怪的是我拿另一个参数做实验“Ij”，无论是否将其定义为全局变量，方程均可求解。
我有点不懂这是为什么，难道是刚性方程？
因为涉及到参数变化的多次计算，不希望手动进行，太繁琐了，做循环计算是我想要的，但第一次都进行不了。请大家给我解答一下。 主程式中没看到有给e1赋值!? 哦，我修改后忘记粘贴了。应该是
clc
clear
tic
global e1
e1_sub=0.5*10^-3;
e=e1_sub;
y0=;
period=1/10; %采样的周期，转速的单位是10Hz,所以周期为“1/10”；
=ode45(@differential_equation,,y0);
本来是不需要使用“e1_sub”的，但是在实际操作时用的是循环，这里只是单一的一次计算。
现在没什么问题了，但还是提示“buzy”，无法运行。 上面写错了，应该是
e1=e1_sub;
回复 板凳 chunshui2003 的帖子
昨晚也试过这个, 没问题! 请问是用我的方程计算的吗？如果是的话，那我想应该是我的matlab的问题，我再找其他人的机器用不同版本的试一下。谢谢你了！ 非常抱歉!刚刚再仔细看, 才发现昨晚是使用1F中的数据试的(% e1=0.5*10-3;%发电机转子偏心距离), 与0.5*10^-3;差很多!
试下发现的确如LZ说的一直busy, 原因不知!
但可以确定与使用或不使用全局变量无关! 是与e1的值相关 ChaChing你好，谢谢你的建议。我的问题已经解决了！主要原因是在龙哥库塔法的初值选取上，经过查看计算过程，发现未知量的数值一直在变化，证明方程的系数传导没有问题。将初值修改后，如，即可进行计算。估计是初值全部为1时，在某一个计算的环节不收敛，所以一直循环下去。希望今后大家遇到类似的问题将初始值更改一下也许问题就解决了！
页:
[1]