ode45、15s解微分方程组疑惑

chunshui2003

%%%%% 微分方程组
function y=fangcheng(t,x)    % 求解定值问题

m=20;             %转子质量
e=0.0009;         %圆盘偏心距
xi=0.13;
xi_z=0.15;
omega_0=16*2*pi;      
omega_e=22*2*pi;      
delta=0.0018;     
kc=25*10^5;       %定子的径向刚度
f=0.1;            %转子与定子之间的摩擦系数
R=0.1;            %碰摩点到转子形心的距离
Fe=9;             %轴向推力的幅值
omega_z=20*2*pi;
g=9.80;           %重力加速度

global omega   

Vr=R*omega+(x(2).*x(5)-x(3).*x(4))/sqrt(x(2).^2+x(4).^2);
x(1)=t;

y=zeros(7,1);
y(1)=t;
y(2)=x(3);
y(3)=e*omega.^2*cos(omega.*x(1))-2.*xi*omega_0.*x(3)-omega_0^2.*x(2)-kc/m*(1-delta/sqrt(x(2).^2+x(4).^2)).*(x(2)-f*Vr/sqrt(Vr^2+x(7).^2).*x(4));
y(4)=x(5);
y(5)=e*omega.^2*sin(omega.*x(1))-2.*xi*omega_0.*x(5)-omega_0^2.*x(4)-kc/m*(1-delta/sqrt(x(2).^2+x(4).^2)).*(x(4)+f*Vr/sqrt(Vr^2+x(7).^2).*x(2));
y(6)=x(7);
y(7)=Fe/m*cos(omega_e*x(1))-g-2*xi_z*omega_z*x(7)-omega_z^2*x(6)-kc/m*f*(sqrt(x(2)^2+x(4)^2)-delta)*x(7)/sqrt(Vr^2+x(7)^2);

%%%%%%%

clc
clear

global omega   
tic

omega=2.51*16*2*pi;
period=2*pi/omega;
tspan=[0:period/100:10*period];
y0=[1 1 1 1 1 1 1 ];
[t,x]=ode45(@fangcheng,tspan,y0);

toc

问题出现在y(7)行的红色部分，之前将方程写错了。“*x(7)”写成了“-x(7)”，但是可以算。结果改过来“*x(7)”反倒是不能算了，一直提示buzy，改为ode15s可以计算，但结果不理想，本来应该出现1000个点，可只有600多个，并且还提示：

Warning: Failure at t=1.615971e-001.
Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.741082e-016) at time t.
我又把时间t的步长给进一步缩小至10^17，但提示“Maximum variable size allowed by the program is exceeded.”

不知道是哪里错了，请大家给指点一下。

[ 本帖最后由 无水1324 于 2010-1-4 10:00 编辑 ]

yangxiaokang

“R=0.1;            %碰摩点到转子形心的距离”？不是很明白什么意思，是定点碰摩，也在做转子局部碰摩，以后大家相互学习:@)

chunshui2003

请问楼上，你对这类刚性微分方程组有什么好的解法吗？可以交流一下。

yangxiaokang

请问最近搞得怎么样了，可以交流相互学习下283707561

chunshui2003

这是刚性微分方程组，ode15s可以算，但是结果不理想。也可以用ode23s算，把精度降低一些，试了一下，是可以的，就是慢，需要等。

无水1324

如果可以算，精度还可以的话，那就只能这样子了。除非自己再根据目前其它的一些刚性算法编程再继续算

chunshui2003

恩，谢谢无水的回复，过年一直没上论坛。

yangzhanwen

刚性的问题怎么都是不好算的，除非根据具体问题自己编程序

chunshui2003

确实是这样，但实际中又避免不了