ode45解多元一次微分方程组问题

ch_j1985

给的程序运行之后除了第一行，其它的全是NaN！
能否把多元一次微分方程组贴出来？

leohust

ode45解一阶微分方程组.doc (77 KB, 下载次数: 15)

2009-2-25 20:39 上传

点击文件名下载附件
一阶微分方程组

一组9元一阶微分方程组，我化简后写了m文件及命令，计算后老是出错，总是得到NaN，附件是原方程，请大家帮忙看看出了什么问题。因为在原来那个帖子贴不出附件，所以另外发了一个贴，还请大家给些指导，谢谢！
function f=coal2(t,y,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20)
f=[6*a1*a2*y(6)*exp(-a4/y(1))/(a3*y(3))+6*a5*(a6^4-y(1)^4)/(a3*y(3)*y(4))-12*a7*(y(1)-y(2))/(y(3)^2*y(4)*a3);
a8*a10*exp(-a11/(a12*y(2)))*y(7)^(-0.3)*y(6)^1.3/a9+a13*a14*exp(-a11/(a12*y(2)))*y(8)^(-0.1)*y(6)^1.85/a9+12*3.14*a15*a7*y(3)*(y(1)-y(2))/(a9*y(5))+12*a7*(a6-y(2))/(a9*a16^2*y(5))+a15*a3*(y(1)-y(2))*0.3*a17*y(9)*exp(-a19/y(1))/(a9*y(5))+a15*a3*(y(1)-y(2))*0.6*a18*y(9)*exp(-a20/y(1))/(a9*y(5))+a15*a3*(y(1)-y(2))*3.14*a1*y(3)^2*y(4)*y(6)*exp(-a4/y(1))/(a9*y(5));
-2*a1*y(6)*exp(-a4/y(1));
-1.8*a17*exp(-a19/y(1))*y(9)/(3.14*y(3)^3)-3.6*a18*y(9)*exp(-a20/y(1))/(3.14*y(3)^3);
    0.3*a15*a17*exp(-a19/y(1))*y(9)+0.6*a15*a18*exp(-a20/y(1))*y(9)+3.14*a15*a1*y(3)^2*y(4)*y(6)*exp(-a4/y(1));
     -8*3.14*a15*y(3)^2*a1*y(4)*y(6)*exp(-a4/y(1))/(3*y(5))-4*a8*y(7)^(-0.3)*y(6)^1.3*exp(-a11/(a12*y(2)))-3*a13*y(8)^(-0.1)*y(6)^1.85*exp(-a11/(a12*y(2)));
    0.3*a17*y(9)*a15*exp(-a19/y(1))/y(5)+a8*y(7)^(-0.3)*y(6)^1.3*exp(-a11/(a12*y(2)));
    0.6*a18*a15*y(9)*exp(-a20/y(1))/y(5)+a13*y(8)^(-0.1)*y(6)^1.85*exp(-a11/(a12*y(2)));
    -a17*exp(-a19/y(1))*y(9)-a18*exp(-a20/y(1))];
end
options=odeset('RelTol',1e-3,'AbsTol',[1e-3],'OutputFcn',@odeplot,'OutputSel',[1]);
t_end=5; a1=7630; a2=32805; a3=1.1; a4=17615; a5=4.7628*10^(-8); a6=1273; a7=0.05; a8=1.5*10^7; a9=1; a10=3.95*10^4;
a11=125000; a12=8.314; a13=1.3*10^9; a14=3.21*10^4; a15=1808; a16=0.01; a17=3.7*10^5; a18=1.46*10^13;a19=8899;a20=30186;
[t,y]=ode45(@coal2,[0 t_end],[300 300 7.5*10^(-5) 1250 1.183 0.23 0 0 6.2*10^(-7)],options,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20);

sogooda

一点想法，仅供参考。
因为原方程的性质不清楚，可以试试ode15s求解器。因为如果待求解的方程存在比较明显的刚性的话，用ode45可能造成求解失败，而ode15s求解成功的可能性更大一些。