用四阶龙哥库塔法求解常微分方程的问题求助

zhihaikou

用四阶龙哥库塔法求解常微分方程的问题求助
用四阶龙哥库塔法求解常微分方程时，尺度比较小，最小在纳米量级。出现以下错误。减小相对误差 绝对误差 步长后，还是没有改观。

Warning: Failure at t=1.287919e-007.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.235165e-022) at time t.

计算到t=1.287919e-007就断掉了。 这是怎么回事？


函数：
function du=eqns3(x,u)
global Tw Tv Hfg den_l Vl Kl e A
du=zeros(4,1);         
du(1)=u(2);            
du(2)=u(3);            
du(3)=u(4);
t=1+u(2)^2;
Pc=A/u(1)^3+e*u(3)/t^1.5;
Tiv=Tv*(1+Pc/(den_l*Hfg));
dK1=u(4)/t^1.5-3*u(1)*u(3)^2*t^2.5;
dK2=-3*(A*u(2)^2/u(1)^2-A*u(3)/u(1)+e*u(1)^2*u(2)*dK1-Vl*Kl*(Tw-Tiv)/(Hfg*u(1)))/(e*u(1)^3);        
du(4)=t^1.5*dK2+9*u(2)*u(3)*u(4)/t+3*u(3)^3/t+15*u(2)^2*u(3)^3/t^2;


主程序：
global Tw Tv Hfg den_l den_v Vl Kl e A f Ri Rg
Tv=333.15;      
Tw=334.15;      
Hfg=2357600;
den_l=983.2;
Vl=4.78e-7;
Kl=0.659;
e=0.06622;
A=5.05e-20;
H0=((Tw/Tv-1)*den_l*Hfg/A)^(-1/3);      
y0=[H0 0 1e-15 0];                           
[x,num_y1]=ode45('eqns3',[0,1e-6],y0);     
plot(x,num_y1(:,1));

sswwss

可能是其中有分母为零 成为不连续点
试一下隐式求解ode15i

zhihaikou

原帖由 sswwss 于 2009-6-13 17:18 发表

                               
登录/注册后可看大图

可能是其中有分母为零 成为不连续点
试一下隐式求解ode15i

非常感谢 我回去试试看