chinamiracle 发表于 2008-5-10 22:45

关于ode45求解方程的问题

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

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

sigma665 发表于 2008-5-10 22:47

回复 楼主 的帖子

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

chinamiracle 发表于 2008-5-11 14:38

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=;
Y0=;
= ode45(@fc,tspan,Y0,options);
plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')
这样写有什么问题。这只是主程序部分，还有一部分是函数文件

sigma665 发表于 2008-5-11 15:54

回复 3楼 的帖子

有什么错误提示

chinamiracle 发表于 2008-5-11 15:55

函数文件大概是
function dy=fc(t,y)
global qi q1 q2 qw qr r qoV1 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 发表于 2008-5-11 15:58

错误提示我也看过了，但我做了啊，到底还错在哪，

提示：??? 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>
= ode45(@fc,tspan,Y0,options);
时间区间不是有了么？第二个不懂，第三个也不知错哪了。能有位朋友加QQ聊下就好了。那样比较方便点我QQ 358079855

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

xjzuo 发表于 2008-5-11 16:02

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

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

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

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

chinamiracle 发表于 2008-5-11 16:09

方程我可以帖出但是比较烦，而且还比第一个长现全部帖出。
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)+*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.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*;
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)-/(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)+*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)-/(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)+*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)-/(y(14)+y(15))*y(16)/(0.05+y(16))*y(18)/(0.01+y(18))*y(22)+*y(26)-*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)+*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)+;
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)+*y(23)-0.2*y(24);
dy(26)=q1/V2*y(13)-3*(qw+qr)/V2*y(26)+*y(26);
方程应该在表达上应该没有问题了，我直接运行过 哦对了，同时谢谢各位的回答。本人在线。其中各方程一些固定参数是要在gui做的界面编辑框中输入。 不知是不是参数这样根本传不过来，有问题，我在每个编辑框的callback下都定义该编辑框全局变量了。如前面程序中要用到的e1 到e 44,并获取了该编辑框值给e1到e44

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

sigma665 发表于 2008-5-11 16:13

回复 8楼 的帖子

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

chinamiracle 发表于 2008-5-11 16:21

原帖由 sigma665 于 2008-5-11 16:13 发表 http://www.chinavib.com/forum/images/common/back.gif
最好用word写，然后截图，发附件上来
能看清楚就可以了，不用word吧，单要全部写过并写对全部方程就要个把小时那样。截不截图应该不是很大影响，方程是推出来的，在本子上。好不容易在matlabm文件中全部表示正确出来，把相关变量　用了y(1)....y(n)进行代替。就这工作都花了我好长时间，写得头和脖子都痛。眼睛都看花了，生怕搞错。
页: [1]
查看完整版本: 关于ode45求解方程的问题