编程论坛 → Matlab

求助！！！用 MATLAB 的 ode45 求解微分方程组

叶子啊发布于 2014-04-18 10:59， 18550 次点击

我写的M文件是：
function dx=p(t,x)
dx=zeros(6,1);
dx(1)=-0.05*x(1)*x(6)+0.11*x(2)*x(5);
dx(2)=0.05*x(1)*x(6)-0.11*x(2)*x(3)-0.215*x(2)*x(6)+1.228*x(3)*x(5);
dx(3)=0.215*x(2)*x(6)-1.228*x(3)*x(5)-0.242*x(3)*x(6)+0.007*x(4)*x(5);
dx(4)=0.242*x(3)*x(6)-0.007*x(4)*x(5);
dx(5)=0.5*x(1)*x(6)-0.11*x(2)*x(5)+0.215*x(2)*x(6)+1.228*x(3)*x(5)+0.242*x(2)*x(6)-0.007*x(4)*x(5);
dx(6)=-0.5*x(1)*x(6)+0.11*x(2)*x(5)-0.215*x(2)*x(6)-1.228*x(3)*x(5)-0.242*x(2)*x(6)+0.007*x(4)*x(5);

命令窗口：>>[t,x]=ode45('p',[0 3000],[1 0 0 0 0 6],0)
出现的错误：      error('MATLAB:odearguments:IncorrectSyntax',...
            ['Correct syntax is ', solver, '(', funstring(ode), ...
             ',tspan,y0,options).']);
和：[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, dataType] = ...
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
nfevals = nfevals + 1;

诚请各位帮忙看看我的问题怎么解决，谢谢啦

13 回复

#2

zvgbaishi2014-04-19 11:05

应该是没有自变量吧
function dx=p(t,x)

好像你的方程组里没有自变量t吧

#3

zvgbaishi2014-04-19 11:07

还有方程组里用dx（1,1）=
dx（2,1）=
.
.
.
.

#4

叶子啊2014-04-19 20:57

回复 2 楼 zvgbaishi

我看书上也没有对时间变量t定义啊，我方程里dx（1）到dx（6）是对t求导

#5

zvgbaishi2014-04-19 23:28

你把题目发一下

#6

zvgbaishi2014-04-19 23:29

我也是新手，交流交流

#7

zvgbaishi2014-04-19 23:39

这是我们老师给的例题 [local]2[/local]你看看

#8

zvgbaishi2014-04-19 23:41

只有本站会员才能查看附件，请登录

#9

叶子啊2014-04-20 11:43

回复 5 楼 zvgbaishi

我的题目是对生物柴油的生产方程的动力学微分方程，用MATLAB的ode45解微分方程并仿真
微分方程是：dT(t)/dt=-0.05*T(t)*A(t)+0.11*W(t)*E(t);
            dW(t)/dt=0.05*T(t)*A(t)-0.11*W(t)*M(t)-0.215*W(t)*A(t)+1.228*M(t)*E(t);
            dM(t)/dt=0.215*W(t)*A(t)-1.228*M(t)*E(t)-0.242*M(t)*A(t)+0.007*G(t)*E(t);
            dG(t)/dt=0.242*M(t)*A(t)-0.007*G(t)*E(t);
            dE(t)/dt=0.5*T(t)*A(t)-0.11*W(t)*E(t)+0.215*W(t)*A(t)+1.228*M(t)*E(t)+0.242*M(t)*A(t)-0.007*G(t)*E(t);
            dA(t)/dt=-dE(t)/dt;
初始条件：T(0)=1,A(0)=6,W(0)=0,M(0)=0,G(0)=0,E(0)=0
谢谢你的热心帮忙，我也不懂，希望也能和你一起探讨！同时也谢谢你发过来的列子！！

#10

叶子啊2014-04-21 10:42

回复 7 楼 zvgbaishi

你好，我根据你给的列子写了这个M文件
function z=fly(t,y)
z(1,:)=-0.05*y(1).*y(6)+0.11*y(2).*y(5);
z(2,:)=0.05*y(1).*y(6)-0.11*y(2).*y(3)-0.215*y(2).*y(6)+1.228*y(3).*y(5);
z(3,:)=0.215*y(2).*y(6)-1.228*y(3).*y(5)-0.242*y(3).*y(6)+0.007*y(4).*y(5);
z(4,:)=0.242*y(3).*y(6)-0.007*y(4).*y(5);
z(5,:)=0.05*y(1).*y(6)-0.11*y(2).*y(5)+0.215*y(2).*y(6)+1.228*y(3).*y(5)+0.242*y(2).*y(6)-0.007*y(4).*y(5);
z(6,:)=-0.05*y(1).*y(6)+0.11*y(2).*y(5)-0.215*y(2).*y(6)-1.228*y(3).*y(5)-0.242*y(2).*y(6)+0.007*y(4).*y(5);
运行以后一直没有出结果，一直处于busy状态，你知道这是为什么吗？我的命令是：>> y0=[1;0;0;0;0;6];
>> [x,y]=ode45('fly',[0 1000],y0)

希望你能帮忙看看，不胜感激！！

#11

zvgbaishi2014-04-21 12:23

function test_cy
[t,y]=ode45(@fn,[0 300],[1 0 0 0 0 6])

function dy=fn(t,y)
dy=zeros(6,1);
dy(1)=-0.05*y(1)*y(6)+0.11*y(2)*y(5);
dy(2)=0.05*y(1)*y(6)-0.11*y(2)*y(3)-0.215*y(2)*y(6)+1.228*y(3)*y(5);
dy(3)=0.215*y(2)*y(6)-1.228*y(3)*y(5)-0.242*y(3)*y(6)+0.007*y(4)*y(5);
dy(4)=0.242*y(3)*y(6)-0.007*y(4)*y(5);
dy(5)=0.5*y(1)*y(6)-0.11*y(2)*y(5)+0.215*y(2)*y(6)+1.228*y(3)*y(5)+0.242*y(3)*y(6)-0.007*y(4)*y(5);
dy(6)=-dy(5);

#12

zvgbaishi2014-04-21 12:27

你的[0 1000]区间太大了好像
function test_cy
[t,y]=ode45(@fn,[0 100],[1 0 0 0 0 6])

function dy=fn(t,y)
dy=zeros(6,1);
dy(1)=-0.05*y(1)*y(6)+0.11*y(2)*y(5);
dy(2)=0.05*y(1)*y(6)-0.11*y(2)*y(3)-0.215*y(2)*y(6)+1.228*y(3)*y(5);
dy(3)=0.215*y(2)*y(6)-1.228*y(3)*y(5)-0.242*y(3)*y(6)+0.007*y(4)*y(5);
dy(4)=0.242*y(3)*y(6)-0.007*y(4)*y(5);
dy(5)=0.5*y(1)*y(6)-0.11*y(2)*y(5)+0.215*y(2)*y(6)+1.228*y(3)*y(5)+0.242*y(3)*y(6)-0.007*y(4)*y(5);
dy(6)=-dy(5);

有结果

只有本站会员才能查看附件，请登录

#13

鸥翔鱼游2014-04-25 14:16

我也是新手，交流交流

#14

LDH山科2017-01-13 21:18

这个帖子帮到了我，必须提出来感谢！！！！！

1