【新手求助】求解微分方程组时为什么会出现这样的错误？

tyut_syz

求解微分方程组时为什么会出现这样的错误？Warning: Failure at t=8.773617e-021.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.117015e-035) at time t.
(Type "warning off MATLAB:ode23tb:IntegrationTolNotMet" to suppress this warning.)
> In C:\MATLAB6p5\toolbox\matlab\funfun\ode23tb.m at line 511。
我反复查验了程序，没有发现哪里有错误，请高手不吝赐教。

tyut_syz

自己顶。
难道高手们不愿意赐教吗？
请大家照顾照顾新手吧。

eight

原帖由 tyut_syz 于 2007-4-23 10:45 发表

                               
登录/注册后可看大图

自己顶。
难道高手们不愿意赐教吗？
请大家照顾照顾新手吧。

你这个是警告，不是错误。如果你知道英文表达的意思，但是对其产生警告的原因不明，我帮不到你，因为我也不懂。建议多看看置顶贴：聚宝盆

tyut_syz

具体程序如下：
% yc=y(1);dycdt=x=y(2);beta=y(3);dbetadt=w=y(4);
function ff=xuan1(t,y,flag)
ro1=7.8*10^3;
ro2=2.13*10^3;
ro3=2.2*10^3;
L=1;
Rc=0.1;
u=0.2;
Alpha=30*pi/180;
Theta=atan(u);
m=pi*Rc^2*ro1*(L+Rc/(3*tan(Alpha)));
k1=(pi*ro2*ro3*cos(Theta))/((ro3-ro2)*(1+u^2));
k2=u*cos(Alpha)^2;
Zc=(6*L^2*tan(Alpha)^2-Rc^2)/(4*tan(Alpha)*(Rc+3*L*tan(Alpha)));
R=Rc+(y(1)+Zc)*tan(Alpha);
Ip=(1/2)*(3*L*m*tan(Alpha)*Rc^2)/(3*L*tan(Alpha)+Rc)+(3/10)*(m*Rc^3)/(3*L*tan(Alpha)+Rc);
% A1=-(sin(Alpha)^2*k1)/(2*m)*R^2;与yc有关
C1=-(sin(Alpha)^2*k1)/(2*m);
A1=C1*R^2;
A2=-(3*k1*k2*sin(Alpha))/m;
A3=(3*k1*k2*sin(Alpha)*cos(Alpha))/m;
A4=-(3*k1*k2^2)/(2*m);
A5=-(k1*k2^3)/(m*sin(Alpha)*cos(Alpha));
A6=(k1*k2^3)/(m*sin(Alpha));
% B1=-(u*k1*sin(Alpha))/(3*Ip)*R^2; 与yc有关
E1=-(u*k1*sin(Alpha))/(3*Ip);
B1=E1*R^2;
B2=(2*u*k1*sin(Alpha)*cos(Alpha)^2)/(3*Ip);
B3=-(2*u*k1*cos(Alpha)^3)/Ip;
% B4=-(u*k1*k2)/Ip*R^2; 与yc有关
E4=-(u*k1*k2)/Ip;
B4=E4*R^2;
B5=(u*k1*k2*cos(Alpha)^2)/Ip;
B6=-(2*u*k1*k2*cos(Alpha)^2)/Ip;
B7=-(u*k1*k2^2*cos(Alpha)^2)/(Ip*sin(Alpha));
B8=-(u*k1*k2^2)/(Ip*sin(Alpha));
B9=(2*u*k1*k2^2*cos(Alpha))/(Ip*sin(Alpha));
cot(Alpha);
cos(Theta);
sin(Alpha+Theta);
tan(Alpha);
initialy1=-Zc-Rc*cot(Alpha);
%以下为求解
   ff1=y(2);
   ff2=A1*y(2)^2+A2*(y(2)^3*sqrt(y(4)^2*R^2+y(2)^2*cos(Alpha)^2))/y(4)^2+...
   (A3+A5)*y(2)^4/y(4)^2+A4*y(2)^4/y(4)^2*log((y(4)^2*R^2+y(2)^2*cos(Alpha)^2)/(y(2)^2*cos(Alpha)^2))+...
   A6*(y(2)^5/(y(4)^2*sqrt(y(4)^2*R^2+y(2)^2*cos(Alpha)^2)));
   ff3=y(3);
   ff4=(B1*y(2)^2/y(4)+B2*y(2)^4/y(4)^3+B8*y(2)^4/y(4)^3)*sqrt(y(4)^2*R^2+y(2)^2*cos(Alpha)^2)+...
       B4*y(2)^3/y(4)+(B5*log(y(4)^2*R^2+y(2)^2*cos(Alpha)^2)+B6*log(y(2))*cos(Alpha)+B9+B3)*y(2)^5/y(4)^3+...
       B7*y(2)^6/(y(4)^3*sqrt(y(4)^2*R^2+y(2)^2*cos(Alpha)^2));
   ff=[ff1;ff2;ff3;ff4];

  [t,y]=ode23tb(@xuan1,[0 0.1],[-0.64355;1000;0;1]);
w=y(:,1)-(-0.64355);
force=diff(y(:,2));
出现的错误为：
      Warning: Failure at t=8.773617e-021.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.117015e-035) at time t.
(Type "warning off MATLAB:ode23tb:IntegrationTolNotMet" to suppress this warning.)
> In C:\MATLAB6p5\toolbox\matlab\funfun\ode23tb.m at line 511
我搞不清为什么？请高手不吝赐教，先谢过了！！！

tyut_syz

我知道是警告，但是它不算呀，我不知道如何改，原因在哪里？

tyut_syz

自己顶吧。难道没有人知道什么原因？

xjzuo

把你的方程用word上传一下看看．
tspan=[0 0.1]区间太小？
从警告看,似乎是积分精度达不到. 看来还是方程本身的问题.

[ 本帖最后由 xjzuo 于 2007-4-25 09:37 编辑 ]

tyut_syz

好像不是时间的问题，我加长时间还是那个问题。附件中是方程组及其初始条件，请您帮忙看看，谢谢了。

xjzuo

我画了一下图形,只有y(:,4)的图形“尚可”,其余都是直线。
%%%----------------------------------------------------------------%%%
clear all
[t,y]=ode23tb(@xuan0,[0 0.1],[-0.64355;1000;0;1]);  % 改为ode15s也无改善.
plot(t,y(:,1),'b',t,y(:,4),'r')
%%%----------------------------------------------------------------%%%
关于积分精度的问题:我在options中调节了积分精度到很高也不行。
将tspan改为[0 8.773617e-021],虽无警告,但结果并无改善。
所以由于方程的奇异性,恐怕要你自己编写RK程序了。

tyut_syz

谢谢您！！！