ode45解多元一次微分方程组问题
请问，我解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',,'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;
%%%%%%%%%%%%%
=ode45('coal2',,,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 编辑 ] =ode45('coal2',,,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18);
少打options
改成=ode45('coal2',,,options，a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18); 谢谢
但是修改后又出些
??? Input argument "t" is undefined.
Error in ==> odeplot at 37
nt = length(t);
>> 你等一下啊
我帮你运行一下，上次回答对不起啊
我没有很仔细的看，对不起 我可以给一下啊我的意见
你没有正确选用解决问题的方法
ode是求解微分方程的，所以面对高阶的问题，我们必须降阶
而在你的文件中，没有降阶的代码，所以是运行不出来的
请告诉我你的原始模型是什么情况,我帮你看看在选用数学方法
相信我可以帮你解决
[ 本帖最后由 ChaChing 于 2009-2-21 23:46 编辑 ]
回复 楼主 leohust 的帖子
楼主的式子很杂, 好好检查一下, 再试看看!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 编辑 ]
回复 5楼 xukai871105 的帖子
原始的模型是描述煤粉燃烧过程的模型，求解着火温度、时间等，利用能量守恒等列出的一系列方程，但是都是一阶的回复 6楼 ChaChing 的帖子
谢谢，里面是有些错误，我继续检查下：）回复 8楼 leohust 的帖子
修改了一下，还会报这样的错？另外，将y的初值代入f中的表达式，出现了Inf和-Inf
??? Error using ==> vertcat
CAT arguments dimensions are not consistent.
Error in ==> coal2 at 2
f=
Error in ==> odearguments at 111
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
楼主的式子很杂, 很容易出错! 昨晚检查好像不只6F列的!
只是个人有点怀疑, 为何式子有误了, 为何matlab未报这个错! 还没时间找这个答案! 有人知道吗?
若如楼上修改後, 都还是Inf和-Inf! 是否表示式子还有错误!
楼主真的要好好检查一下了 检查修改后还是出些
??? Input argument "t" is undefined.
Error in ==> odeplot at 37
nt = length(t);
不知道什么原因
回复 11楼 leohust 的帖子
LS最後的function f=coal2为何? 贴出吧总感觉式这里有问题?
回复 12楼 ChaChing 的帖子
coal2为函数的名字，应该不会有问题，我把式中的常量a1，a2。。。用数字计算代入后仍然出现??? Input argument "t" is undefined.
Error in ==> odeplot at 37
nt = length(t);
真是很郁闷啊，不知道什么原因
回复 13楼 leohust 的帖子
我改了之后的，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',,'OutputFcn',@odeplot,'OutputSel',);
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;
=ode45(@coal2,,,options,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18);
回复 14楼 ch_j1985 的帖子
谢谢！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
页:
[1]
2