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

leohust

请问，我解9元一次微分方程组时
m文件：function f=coal2(t,y,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18)
f=[6*a1*a2*y(6)*exp(-a4/y(1))/(a3*y(3))+6*a5*(a6^4-y(1)^4)/(a3*y(3))-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*a15*a7*y(3)*(y(1)-y(2))/(a9*a16^3*y(5))+6*a15*a3(y(1)-a6)*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3(y(1)-a6)*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3(y(1)-a6)*y(3)^2*a1*exp(-a4/y(1))*y(4)*y(6)/(a9*a16^3*y(5)+12*a7*(a6-y(2))/(a9*a16^2*y(5)));
2*a1*y(6)*exp(-a4/y(1));
-6*0.3*a17*exp(-a11/(a12*y(1))*y(9)/(3.14*y(3)^3)+4.2*a18*y(9)*exp(-a11/(a12*y(1)))/(3.14*y(3)^3));
    6*a15*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(3.14*a16^3)+6*a15*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(3.14*a16^3)+6*a15*a1*y(3)^2*y(4)*y(6)/(a16^3);
     -16*a15*y(3)^2*a1*y(4)*y(6)*exp(-a4/y(1))-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)));
    1.8*a17*y(9)*a15*exp(-a11/(a12*y(1)))/(3.14*a16^3*y(5))+a8*y(7)^(-0.3)*y(6)^1.3*exp(-a11/(a12*y(2)));
    4.2*a18*a151*y(9)*exp(-a11/(a12*y(1)))/(3.14*a16^3*y(5))+a13*y(8)^(-0.1)*y(6)^1.85*exp(-a11/(a12*y(2)));
    -a17*exp(-a11/(a12*y(1)))*y(9)-a18*exp(-a11/(a12*y(1)))];
运行命令
options=odeset('RelTol',1e-3,'AbsTol',[1e-3],'OutputFcn',odeplot,'OutputSel',y(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=2243; a16=0.01; a17=3.7*10^5; a18=1.46*10^13;
%%%%%%%%%%%%%
[t,y]=ode45('coal2',[0 t_end],[300 300 7.5*10^(-5) 1250 1.183 0.23 0 0 6.2*10^(-7)],a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18);

检查m文件总是出现 Input argument "y" is undefined.
请问大家是什么原因？？
另外不知道上述程序是否有问题，请大家帮忙看下，万分感谢！

[ 本帖最后由 ChaChing 于 2009-2-21 23:56 编辑 ]

xukai871105

[t,y]=ode45('coal2',[0 t_end],[300 300 7.5*10^(-5) 1250 1.183 0.23 0 0 6.2*10^(-7)],a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18);
少打options
改成[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);

leohust

谢谢
但是修改后又出些
??? Input argument "t" is undefined.
Error in ==> odeplot at 37
  nt = length(t);
>>

xukai871105

你等一下啊
我帮你运行一下，上次回答对不起啊
我没有很仔细的看，对不起

xukai871105

我可以给一下啊我的意见

你没有正确选用解决问题的方法

ode是求解微分方程的，所以面对高阶的问题，我们必须降阶
而在你的文件中，没有降阶的代码，所以是运行不出来的

请告诉我你的原始模型是什么情况,我帮你看看在选用数学方法
相信我可以帮你解决

[ 本帖最后由 ChaChing 于 2009-2-21 23:46 编辑 ]

ChaChing

楼主的式子很杂, 好好检查一下, 再试看看!
ex:
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*a15*a7*y(3)*(y(1)-y(2))/(a9*a16^3*y(5))+6*a15*a3(y(1)-a6)*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3(y(1)-a6)*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3(y(1)-a6)*y(3)^2*a1*exp(-a4/y(1))*y(4)*y(6)/(a9*a16^3*y(5)+12*a7*(a6-y(2))/(a9*a16^2*y(5))); ...

[ 本帖最后由 ChaChing 于 2009-2-22 00:24 编辑 ]

leohust

原始的模型是描述煤粉燃烧过程的模型，求解着火温度、时间等，利用能量守恒等列出的一系列方程，但是都是一阶的

leohust

谢谢，里面是有些错误，我继续检查下：）

ch_j1985

修改了一下，还会报这样的错？
另外，将y的初值代入f中的表达式，出现了Inf和-Inf
??? Error using ==> vertcat
CAT arguments dimensions are not consistent.
Error in ==> coal2 at 2
f=[6*a1*a2*y(6)*exp(-a4/y(1))/(a3*y(3))+6*a5*(a6^4-y(1)^4)/(a3*y(3))-12*a7*(y(1)-y(2))/(y(3)^2*y(4)*a3);
Error in ==> odearguments at 111
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

ChaChing

楼主的式子很杂, 很容易出错! 昨晚检查好像不只6F列的!
只是个人有点怀疑, 为何式子有误了, 为何matlab未报这个错! 还没时间找这个答案! 有人知道吗?
若如楼上修改後, 都还是Inf和-Inf! 是否表示式子还有错误!
楼主真的要好好检查一下了

leohust

检查修改后还是出些
??? Input argument "t" is undefined.
Error in ==> odeplot at 37
  nt = length(t);

不知道什么原因

ChaChing

LS最後的function f=coal2为何? 贴出吧
总感觉式这里有问题?

leohust

coal2为函数的名字，应该不会有问题，我把式中的常量a1，a2。。。用数字计算代入后仍然出现
??? Input argument "t" is undefined.

Error in ==> odeplot at 37
  nt = length(t);

真是很郁闷啊，不知道什么原因

ch_j1985

我改了之后的，LZ可以参考一下，不过会报9楼的错误

function f=coal2(t,y,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18)
f=[6*a1*a2*y(6)*exp(-a4/y(1))/(a3*y(3))+6*a5*(a6^4-y(1)^4)/(a3*y(3))-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*a15*a7*y(3)*(y(1)-y(2))/(a9*a16^3*y(5))+6*a15*a3*(y(1)-a6)*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)...
+6*a15*a3*(y(1)-a6)*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3*(y(1)-a6)*y(3)^2*a1*exp(-a4/y(1))*y(4)*y(6)/(a9*a16^3*y(5)...
+12*a7*(a6-y(2))/(a9*a16^2*y(5)));
2*a1*y(6)*exp(-a4/y(1));
-6*0.3*a17*exp(-a11/(a12*y(1))*y(9)/(3.14*y(3)^3)+4.2*a18*y(9)*exp(-a11/(a12*y(1)))/(3.14*y(3)^3));
    6*a15*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(3.14*a16^3)+6*a15*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(3.14*a16^3)+6*a15*a1*y(3)^2*y(4)*y(6)/(a16^3);
     -16*a15*y(3)^2*a1*y(4)*y(6)*exp(-a4/y(1))-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)));
    1.8*a17*y(9)*a15*exp(-a11/(a12*y(1)))/(3.14*a16^3*y(5))+a8*y(7)^(-0.3)*y(6)^1.3*exp(-a11/(a12*y(2)));
    4.2*a18*a15*y(9)*exp(-a11/(a12*y(1)))/(3.14*a16^3*y(5))+a13*y(8)^(-0.1)*y(6)^1.85*exp(-a11/(a12*y(2)));
    -a17*exp(-a11/(a12*y(1)))*y(9)-a18*exp(-a11/(a12*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=2243; a16=0.01; a17=3.7*10^5; a18=1.46*10^13;
[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);

leohust

谢谢！M文件漏写了点，检查m文件的话会出现??? Input argument "a1" is undefined. 不检查，直接运行命令的话，可以出现图像，但是只有初始值，我想得到y1随时间变化的曲线，不知道该怎么写，本人新手，还请多多指教，非常感谢！！！
function f=coal2(t,y,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18)
f=[6*a1*a2*y(6)*exp(-a4/y(1))/(a3*y(3))+6*a5*(a6^4-y(1)^4)/(a3*y(3))-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*a15*a7*y(3)*(y(1)-y(2))/(a9*a16^3*y(5))+6*a15*a3*(y(1)-a6)*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3*(y(1)-a6)*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(y(5)*3.14*a16^3*a9)+6*a15*a3*(y(1)-a6)*y(3)^2*a1*exp(-a4/y(1))*y(4)*y(6)/(a9*a16^3*y(5))+12*a7*(a6-y(2))/(a9*a16^2*y(5));
2*a1*y(6)*exp(-a4/y(1));
-6*0.3*a17*exp(-a11/(a12*y(1))*y(9)/(3.14*y(3)^3)-4.2*a18*y(9)*exp(-a11/(a12*y(1)))/(3.14*y(3)^3));
    6*a15*0.3*a17*exp(-a11/(a12*y(1)))*y(9)/(3.14*a16^3)+6*a15*0.7*a18*exp(-a11/(a12*y(1)))*y(9)/(3.14*a16^3)+6*a15*a1*y(3)^2*y(4)*y(6)*exp(-a4/y(1))/(a16^3);
     -16*a15*y(3)^2*a1*y(4)*y(6)*exp(-a4/y(1))/(y(5)*a16^3)-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)));
    1.8*a17*y(9)*a15*exp(-a11/(a12*y(1)))/(3.14*a16^3*y(5))+a8*y(7)^(-0.3)*y(6)^1.3*exp(-a11/(a12*y(2)));
    4.2*a18*a15*y(9)*exp(-a11/(a12*y(1)))/(3.14*a16^3*y(5))+a13*y(8)^(-0.1)*y(6)^1.85*exp(-a11/(a12*y(2)));
    -a17*exp(-a11/(a12*y(1)))*y(9)-a18*exp(-a11/(a12*y(1)))];
end