ode45、15s解微分方程组疑惑
%%%%% 微分方程组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=;
y0=;
=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 编辑 ]
回复 楼主 chunshui2003 的帖子
“R=0.1; %碰摩点到转子形心的距离”？不是很明白什么意思，是定点碰摩，也在做转子局部碰摩，以后大家相互学习:@) 请问楼上，你对这类刚性微分方程组有什么好的解法吗？可以交流一下。回复 板凳 chunshui2003 的帖子
请问最近搞得怎么样了，可以交流相互学习下283707561 这是刚性微分方程组，ode15s可以算，但是结果不理想。也可以用ode23s算，把精度降低一些，试了一下，是可以的，就是慢，需要等。回复 5楼 chunshui2003 的帖子
如果可以算，精度还可以的话，那就只能这样子了。除非自己再根据目前其它的一些刚性算法编程再继续算 恩，谢谢无水的回复，过年一直没上论坛。 刚性的问题怎么都是不好算的，除非根据具体问题自己编程序 确实是这样，但实际中又避免不了
页:
[1]