关于ode45求解方程的问题

chinamiracle

我想问下用ode45可不可以解决由26个方程组成的微分方程组。每单个方程都有参数，这些参数由界面框调入。初值也要从界面框的编辑框中调用。请问各位如何实现。

[ 本帖最后由 eight 于 2008-5-12 19:12 编辑 ]

sigma665

应该可以
你可以把方程用图片的格式附件上来
让大家一起看

chinamiracle

global e20 e19 e21 e22 e23 e24 e25 e26 e27 e29 e30 e31 e28 e33 e32 e34 e35 e36 e37 e38 e39 e40 e42 e43 e44 e41
%  定义各初值全局变量
global e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 e14 e15 e16 e17 e18
%  定义各固定参数全局变量
qi=e1;
q1=e2;
q2=e3;
qw=e4;
qr=e5;
r=e6;
qo=e7;
V1=e8;
V2=e9;
Sa_i=e10;
Sf_i=e11;
Snh4_i=e12;
Sno3_i=e13;
Spo4_i=e14;
Sn2_i=e15;
XI_i=e16;
Xs_i=e17;
Xh_i=e18;
%给方程中各变量付固定参数值
t0=0; tf=60*60*5;
tspan=[t0,tf];
Y0=[e20 e19 e21 e22 e23 e24 e25 e26 e27 e29 e30 e31 e28 e33 e32 e34 e35 e36 e37 e38 e39 e40 e42 e43 e44 e41];
[t,Y] = ode45(@fc,tspan,Y0,options);
plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')
这样写有什么问题。这只是主程序部分，还有一部分是函数文件

sigma665

有什么错误提示

chinamiracle

函数文件大概是
function dy=fc(t,y)
global qi q1 q2 qw qr r qo  V1 V2 Sa_i Sf_i Snh4_i Sno3_i Spo4_i Sn2_i XI_i Xs_i Xh_i
dy=zeros(26,1);
下面是26个方程（我只列出一个）
dy(1)=qi/V1*Sf_i-q1/V1*y(1)+[(0.15+1.8*y(4))/(0.1*y(9)+y(8))*y(8)-7.632*y(4)*y(1)/(4+y(1))*y(1)/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))-1.5*y(1)/(20+y(1))]*y(9)/(0.5+y(4));
问这两个里面到底要如何才可以进行求解。小弟的设想都用过了，出问题又不知错在哪，说不出的不爽请达人给点建议

chinamiracle

提示：??? Error using ==> funfun\private\odearguments
When the first argument to ode45 is a function handle, the tspan and y0 arguments must be supplied.

Error in ==> <a href="error:D:\MATLAB7\toolbox\matlab\funfun\ode45.m,173,1">ode45 at 173</a>
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

Error in ==> <a href="error:D:\MATLAB7\work\bishe\zhengpin\compute.m,27,1">compute at 27</a>
[t,Y] = ode45(@fc,tspan,Y0,options);
时间区间不是有了么？第二个不懂，第三个也不知错哪了。能有位朋友加QQ聊下就好了。那样比较方便点我QQ 358079855

[ 本帖最后由 chinamiracle 于 2008-5-11 16:00 编辑 ]

xjzuo

这种问题以前讨论过，原因可能是ode45不一定能解这么多微分方程；当然，你也可以尝试，有时候是可以解的．

基本上你的语法没有明显问题，若要细究，则最好将方程贴全，尤其是公式要贴出来，这样别人才可以较快地提出有效建议．

另：options你似乎没有设置（也可先置为默认）．

[ 本帖最后由 xjzuo 于 2008-5-11 16:04 编辑 ]

chinamiracle

方程我可以帖出但是比较烦，而且还比第一个长现全部帖出。
dy(1)=qi/V1*Sf_i-q1/V1*y(1)+[(0.15+1.8*y(4))/(0.1*y(9)+y(8))*y(8)-7.632*y(4)*y(1)/(4+y(1))*y(1)/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))-1.5*y(1)/(20+y(1))]*y(9)/(0.5+y(4));
dy(2)=qi/V1*Sa_i-q1/V1*y(2)+[1.5*y(1)/(20+y(1))-7.632*y(4)*y(2)/(4+y(2))*y(2)/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))]*y(9)/(0.5+y(4))-3*y(2)/(4+y(2))*y(11)/(0.01*y(10)+y(11))*y(10)+0.2*y(12);
dy(3)=qi/V1*Snh4_i-q1/V1*y(3)+(0.15+1.8*y(4))/(50+100*y(4))*y(8)/(0.1*y(9)+y(8))*y(9)-[0.1056*y(1)*y(1)+0.336*y(2)*y(2)]/(0.5+y(4))/(4+y(1))*y(4)/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))*y(9)+0.045/(0.5+y(4))*y(1)/(20+y(1))*y(9)+0.031*[0.4*y(9)+0.2*y(10)+0.15*y(13)];
dy(4)=qi/V1*Sno3_i-q1/V1*y(4)-0.96*y(4)/(0.5+y(4))*(y(2)*y(2)+y(1)*y(1))/(4+y(1))/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))*y(9);
dy(5)=qi/V1*Spo4_i-q1/V1*y(5)-[0.0192*y(1)*y(1)+0.096*y(2)*y(2)]/(4+y(1))*y(4)/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))*y(9)+0.015/(0.5+y(4))*y(1)/(20+y(1))*y(9)+0.004*y(9)+0.002*y(10)+0.2*y(11)+0.0015*y(13)+1.2*y(2)/(4+y(2))*y(11)/(0.01*y(10)+y(11))*y(10);
dy(6)=qi/V1*Sn2_i-q1/V1*y(6)+1.008*y(4)/(0.5+y(4))*(y(1)*y(1)+y(2)*y(2))/(4+y(11))/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.05+y(5))*y(9);
dy(7)=qi/V1*XI_i+3*qr/V1*y(20)-q1/V1*y(7)+0.04*y(9)+0.02*y(10)+0.015*y(13);
dy(8)=qi/V1*Xs_i+3*qr/V1*y(21)-q1/V1*y(8)+(0.15+1.8*y(4))/(0.5+y(4))*y(8)/(0.1*y(9)+y(8))*y(9)+0.36*y(9)+0.18*y(10)+0.135*y(13);
dy(9)=qi/V1*Xh_i+3*qr/V1*y(22)-q1/V1*y(9)+4.8*y(4)/(0.5+y(4))*(y(2)^2+y(1)^2)/(4+y(1))/(y(1)+y(2))*y(3)/(0.05+y(3))*y(5)/(0.01+y(5))*y(9)-0.4*y(9);
dy(10)=3*qr/V1*y(23)-q1/V1*y(10)-0.2*y(10);
dy(11)=3*qr/V1*y(24)-q1/V1*y(11)-1.2*y(2)/(4+y(2))*y(11)/(0.01*y(10)+y(11))*y(10)-0.2*y(11);
dy(12)=3*qr/V1*y(25)-q1/V1*y(12)+3*y(2)/(4+y(2))*y(11)/(0.01*y(10)+y(11))*y(10)-0.2*y(12);
dy(13)=3*qr/V1*y(26)-q1/V1*y(13)-0.15*y(13);

dy(14)=q1/V2*y(1)-qo/V2*y(14)+[3*y(21)/0.1*y(22)+y(21)-9.54*y(14)/(4+y(14))*y(14)/(y(14)+y(15))*y(18)/(0.01+y(18))]*y(22);
dy(15)=q1/V2*y(2)-qo/V2*y(15)-9.45*y(15)/(4+y(15))*y(15)/(y(14)+y(15))*y(16)/(0.05+y(16))*y(22)-3*y(15)/(4+y(14))*y(24)/(0.01*y(23)+y(24))*y(23)+0.2*y(25);
dy(16)=q1/V2*y(3)-qo/V2*y(16)+0.03*y(21)/(0.1*y(22)+y(21))*y(22)-[0.132*y(14)^2/(4+y(14))+0.42*y(15)^2/(4+y(15))]/(y(14)+y(15))*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))*y(22)+[-0.07*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))*y(25)/(0.01*y(23)+y(25))+0.0062]*y(23)+0.0124*y(22)+[0.00465-4.24*y(16)/(1+y(16))*y(18)/(0.01+y(18))]*y(26);
dy(17)=q1/V2*y(4)-qo/V2*y(17)+4.17*y(16)/(1+y(16))*y(18)/(0.01+y(18))*y(26);
dy(18)=q1/V2*y(5)-qo/V2*y(18)-[0.024*y(14)^2/(4+y(14))+0.12*y(15)^2/(4+y(15))]/(y(14)+y(15))*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))*y(22)+[0.0015-0.02*y(16)/(1+y(16))*y(18)/(0.01+y(18))]*y(26)-[0.02*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))+1.5*y(18)/(0.2+y(18))*(0.34-y(24)/y(23))/(0.36-y(24)/y(23))]*y(25)/(0.01*y(23)+y(25))*y(23)+0.002*y(23)+0.004*y(22)+0.2*y(24)+1.2*y(15)/(4+y(15))*y(24)/(0.01*y(23)+y(24))*y(23);
dy(19)=q1/V2*y(6)-qo/V2*y(19);
dy(20)=q1/V2*y(7)-3*(qw+qr)/V2*y(20)+0.04*y(22)+0.02*y(23)+0.015*y(26);
dy(21)=q1/V2*y(8)-3*(qw+qr)/V2*y(21)+[0.36-3*y(21)/(0.1*y(22)+y(21))]*y(22)+0.18*y(23)+0.135*y(26);
dy(22)=q1/V2*y(9)-3*(qw+qr)/V2*y(22)+[(y(14)^2/(4+y(14))+y(15)^2/(4+y(15)))*6/(y(14)+y(15))*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))-0.4]*y(22);
dy(23)=q1/V2*y(10)-3*(qw+qr)/V2*y(23)+[y(16)/(0.05+y(16))*y(18)/(0.01+y(18))*y(25)/(0.01*y(23)+y(25))];
dy(24)=q1/V2*y(11)-3*(qw+qr)/V2*y(24)+[-1.2*y(15)/(4+y(15))*y(24)/(0.01*y(23)+y(24))+1.5*y(18)/(0.2+y(18))*y(25)/(0.01*y(23)+y(25))*(0.34-y(24)/y(23))/(0.36-y(24)/y(23))]*y(23)-0.2*y(24);
dy(25)=q1/V2*y(12)-3*(qw+qr)/V2*y(25)+[3*y(15)/(4+y(15))*y(24)/(0.01*y(23)+y(24))-0.3*y(18)/(0.2+y(18))*y(25)/(0.01*y(23)+y(25))*(0.34-y(24)/y(23))/(0.36-y(24)/y(23))-1.6*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))*y(25)/(0.01*y(23)+y(25))]*y(23)-0.2*y(24);
dy(26)=q1/V2*y(13)-3*(qw+qr)/V2*y(26)+[y(16)/(1+y(16))*y(18)/(0.01+y(18))-0.15]*y(26);
方程应该在表达上应该没有问题了，我直接运行过 哦对了，同时谢谢各位的回答。本人在线。其中各方程一些固定参数是要在gui做的界面编辑框中输入。 不知是不是参数这样根本传不过来，有问题，我在每个编辑框的callback下都定义该编辑框全局变量了。如前面程序中要用到的e1 到e 44,并获取了该编辑框值给e1到e44

[ 本帖最后由 chinamiracle 于 2008-5-11 16:28 编辑 ]

sigma665

最好用word写，然后截图，发附件上来

chinamiracle

原帖由 sigma665 于 2008-5-11 16:13 发表

                               
登录/注册后可看大图

最好用word写，然后截图，发附件上来

能看清楚就可以了，不用word吧，单要全部写过并写对全部方程就要个把小时那样。截不截图应该不是很大影响，方程是推出来的，在本子上。好不容易在matlabm文件中全部表示正确出来，把相关变量　用了y(1)....y(n)进行代替。就这工作都花了我好长时间，写得头和脖子都痛。眼睛都看花了，生怕搞错。