[转载]MATLAB解微分方程
用matlab时间也不短了，可是一直没有接触过微分方程。这次看看书，学习学习，记点儿笔记。
1.可以解析求解的微分方程。
调用格式为：
y=dsolve(f1,f2,...,fmO;
y=dsolve(f1,f2,...,fm,'x');
如下面的例子,求解了微分方程
u=exp(-5*t)*cos(2*t-1)+5;
uu=5*diff(u,t,2)+4*diff(u,t)+2*u;
y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=87*exp(-5*t)*cos(2*t-1)+92*exp(-5*t)*sin(2*t-1)+10'])
yc=latex(y)
将yc的内容copy到latex中编译，得到结果。
关于Matlab的微分方程，直到今天才更新第2篇，实在是很惭愧的事——因为原因都在于太懒惰，而不是其他的什么。
在上一篇中，我们使用dsolve可以解决一部分能够解析求解的微分方程、微分方程组，但是对于大多数微分方程（组）而言不能得到解析解，这时数值求解也就是没有办法的办法了，好在数值解也有很多的用处。
数值分析方法中讲解了一些Eular法、 Runge-Kutta
法等一些方法，在matlab中内置的ode求解器可以实现不同求解方法的相同格式的调用，而不必太关心matlab究竟是用什么算法完成的。
这一回我们来说明ode45求解器的使用方法。
1.ode45求解的上手例子：
求解方程组
Dx=y+x(1-x^2-y^2);
Dy=-x+y*(1-x^2-y^2)
初值x=0.1;y=0.2;
先说明一下最常用的ode45调用方式，和相应的函数文件定义格式。
[t,x]=ode45(odefun,tspan,x0);
其中，Fun就是导函数，tspan为求解的时间区间（或时间序列，如果采用时间序列，则必须单调），x0为初值。
这时，函数文件可以采用如下方式定义
function dx=odefun(t,x)
对于上面的小例子，可以用如下的程序求解。
function jixianhuan
x0=[0.1;0.2];
[t,x]=ode45(@jxhdot,[0,100],x0);
plot(x(:,1),x(:,2))
function dx=jxhdot(t,x)
x(2)+x(1).*(1-x(1).^2-x(2).^2);
-x(1)+x(2).*(1-x(1).^2-x(2).^2)
&2.终值问题
tspan可以是递增序列，也可以为递减序列，若为递减则可求解终值问题。
[t,x]=ode45(@zhongzhiode,[3,0],[1;0;2]);plot(t,x)
function dx=zhongzhiode(t,x)
dx=[2*x(2)^2-2;
-x(1)+2*x(2)*x(3)-1;
-2*x(2)+2*x(3)^2-4];
options = odeset('name1',value1,'name2',value2,...)
[t,x]=solver(@fun,tspan,x0,options)
通过odeset设置options
第一，通过求解选项的设置可以改善求解精度，使得原本可能不收敛的问题收敛。
options=odeset('RelTol',1e-10);
第二，求解形如M(t,x)x'=f(t,x)的方程。
例如，方程
x'=-0.2x+yz+0.3xy
y'=2xy-5yz-2y^2
可以变形为
0& 0][x'] [-0.2x+yz+0.3xy]
0][y']=[2xy-5yz-2y^2& ]
[x+y+z-2&&&&&&
这样就可以用如下的代码求解该方程
function mydae
M=[1 0 0;0 1 0;0 0 0];
options=odeset('Mass',M);
x0=[1.6,0.3,0.1];
[t,x]=ode15s(@daedot,[0,1.5],x0,options);plot(t,x)
function dx=daedot(t,x)
-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);
2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);
x(1)+x(2)+x(3)-2];
4.带附加参数的ode45
有时我们需要研究微分方程组中的参数对于解的影响，这时采用带有参数的ode45求解会使求解、配合循环使用，可以使得求解的过程更加简捷。
使用方法：只需将附加参数放在options的后面就可以传递给odefun了。
看下面的例子。
function Rossler
a=[0.2,0.2];
b=[0.2,0.5];
c=[5.7,10];
x0=[0 0 0];
for jj=1:2
[t,x]=ode45(@myRossler,[0,100],x0,[],a(jj),b(jj),c(jj));
plot3(x(:,1),x(:,2),x(:,3));
function dx=myRossler(t,x,a,b,c)
-x(2)-x(3);
x(1)+a*x(2);
b+(x(1)-c)*x(3)];
5. 刚性方程的求解
刚性方程就是指各个自变量的变化率差异很大，会造成通常的求解方法失效。
这是matlab中自带的一个例子，使用ode15s求解，如果用ode45求解就会出现错误。
function myode15study
[t,Y] = ode15s(@vdp0],[2 0]);
plot(T,Y(:,1),'-o')
plot(Y(:,1),Y(:,2))
function dy = vdp1000(t,y)
zeros(2,1);&&&
dy(1) = y(2);
dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);
6.高阶微分方程的求解
通常的方法是进行变量替换，将原方程降阶，转换成更多变量的一阶方程组进行求解。
在这个例子里我们求解一个动力学系统里最常见的一个运动方程
,其中f=sin(t)
function myhighoder
x0=zeros(6,1);
[t,x]=ode45(@myhigh,[0,100],x0);
plot(t,x(:,1))
function dx=myhigh(t,x)
f=[sin(t);0;0];;
C=eye(3)*0.1;
K=eye(3)-0.5*diag(ones(2,1),1)-0.5*diag(ones(2,1),-1);
dx=[x(4:6);inv(M)*(f-C*x(4:6)-K*x(1:3))];
7.延迟微分方程
matlab提供了dde23求解非中性微分方程。dde23的调用格式如下：
sol = dde23(ddefun,lags,history,tspan)
lags是延迟量，比如方程中包含y1(t-0.2)和y2(t-0.3)则可以使用lags=[0.2,0.3]。
这里的ddefun必须采用如下的定义方式：
dydt = ddefun(t,y,Z)
其中的Z(:,1)就是y(t-lags(1)),Z(:,2)就是y(t-lags(2)）...
下面是个使用dde23求解延迟微分方程的例子。
function mydde23study
The differential equations
y'_1(t) = y_1(t-1)&
y'_2(t) = y_1(t-1)+y_2(t-0.2)
y'_3(t) = y_2(t)
are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) =
lags=[1,0.2];
history=[1;1;1];
tspan=[0,5];
sol = dde23(@myddefun,lags,history,tspan)
plot(sol.x,sol.y)
function dy = myddefun(t,y,Z)
Z(1)+Z(2,2);
8.ode15i求解隐式微分方程
[T,Y] = ode15i(odefun,tspan,y0,yp0)
yp0为y'的初值。
odefun的格式如下& dy =
odefun(t,y,yp),yp表示y'，而方程中应该使得f(t,y,y')=0
function myodeIMP
The problem is
y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)
y(2)' =& 0.04*y(1) - 1e4*y(2)*y(3) -
3e7*y(2)^2
y(3)' =& 3e7*y(2)^2
It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3)
to steady state.&
y0=[1;0;0];
fixed_y0=[1;1;1];
yp0=[0 0 0];
fixed_yp0=[];
[y0mod,yp0mod]=decic(@myodefunimp,0,y0,fixed_y0,yp0,fixed_yp0);
tspan=[0, logspace(-6,6)];
[t,y] = ode15i(@myodefunimp,tspan,y0mod,yp0mod);
y(:,2)=1e4*y(:,2);
semilogx(t,y)
function res=myodefunimp(t,y,yp)
-yp(1)-0.04*y(1)+1e4*y(2)*y(3);
-yp(2)+0.04*y(1)-1e4*y(2)*y(3)-3e7*y(2)^2;
-yp(3)+3e7*y(2)^2;
这次要接触一个新的求解ode的方法，就是使用simulink的积分器求解。
1.还是做我们研究过的一个例子（在初识matlab微分方程(2)中采用的）。
Dx=y+x(1-x^2-y^2);
Dy=-x+y*(1-x^2-y^2)
初值x=0.1;y=0.2;
积分器中设置初始条件；f(u)中指定Dx,Dy的计算公式。
运行这个仿真，scope中可以看到两个变量的时程如下：
在WorkSpace里可以得到tout和yout，执行plot(yout(:,1),yout(:,1))得到与ode45求解相似的结果如下
&2.这部分解决一个使用ode求解器dde23没法求解的一类延迟微分方程(中性微分方程)。
形如x'(t)=f(x'(t-t1),x(t),x(t-t2),x(t-t3))这类方程。dde23是无法求解的，但是可以借助simulink仿真求解。
看下面的这个例子。
x'(t)=A1*x(t-t1)+A2*x'(t-t2)+B*u(t)
t1=0.15;t2=0.5
A2=[0.02&&&&0&&&&
62]&&&&&&&&
&&&[207&&-207&
113]&&&&&&&&&[0&&&&&&
0.04]&&&&&&[2]&&&
在continuous里找到transport Delay,就可以实现对于信号的延迟，因此可以建立如下仿真模型
从而在scope中可以得到如下仿真结果
OK~初识微分方程到了这里我想应该可以做个终结，因为我想作为零基础的材料来看，到这里也就可以了。以后还可能再有微分方程的内容，还请感兴趣的朋友多捧场吧。
最后，大力推荐一本书薛定宇老师的《高等应用数学问题的Matlab求解》，确实很经典。学习Matlab的时间也不算短了，可是每次翻看这本书总是能让我有温故而知新的感觉，是我目前见过的最好的Matlab书。强烈推荐!（对于从来没有接触过matlab的人来说或许有点儿难，但是如果你以后要用matlab的话买一本绝对不会后悔的。）
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。1被浏览23分享邀请回答暂时还没有回答，开始写第一个回答欢迎加入我们，一同切磋技术 &
用户名： &&&
密　码： &
共有 18197 人关注过本帖
标题：求助！！！用 MATLAB 的 ode45 求解微分方程组
等　级：新手上路
结帖率：100%
&&已结贴√
&&问题点数：10&&回复次数：13&&&
求助！！！用 MATLAB 的 ode45 求解微分方程组
我写的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;
诚请各位帮忙看看我的问题怎么解决，谢谢啦
搜索更多相关主题的帖子:
等　级：论坛游民
专家分：10
&&得分:10&
应该是没有自变量吧
function dx=p(t,x)
好像你的方程组里没有自变量t吧
等　级：论坛游民
专家分：10
还有方程组里用dx（1,1）=
dx（2,1）=
等　级：新手上路
回复 2 楼 zvgbaishi
我看书上也没有对时间变量t定义啊，我方程里dx（1）到dx（6）是对t求导
等　级：论坛游民
专家分：10
你把题目发一下
等　级：论坛游民
专家分：10
我也是新手，交流交流
等　级：论坛游民
专家分：10
这是我们老师给的例题&&[local]2[/local]你看看
等　级：论坛游民
专家分：10
附件: 您没有浏览附件的权限，请
等　级：新手上路
回复 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)/
初始条件：T(0)=1,A(0)=6,W(0)=0,M(0)=0,G(0)=0,E(0)=0
谢谢你的热心帮忙，我也不懂，希望也能和你一起探讨！同时也谢谢你发过来的列子！！
等　级：新手上路
回复 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)
希望你能帮忙看看，不胜感激！！
版权所有，并保留所有权利。
Powered by , Processed in 0.069831 second(s), 7 queries.
Copyright&, BCCN.NET, All Rights Reserved如何用matlab解微分方程：dx/dt=x(t)*(1-X(t-1))._百度知道
如何用matlab解微分方程：dx/dt=x(t)*(1-X(t-1)).
急需，那位高手帮忙一下，谢了！！！
x(0)=0吧。
我有更好的答案
//h;tau=1;history=0;%延迟问题求解&%延迟函数接下来命令求解&&&nbsp.&nbsp；所以需要之到一个初始条件x(0)的值，但是是数值解法,z)tau=z;%定义延迟时间dx=x*(1-tau).hiphotos.baidu.com/zhidao/wh%3D600%2C800/sign=f543fc414a90fe3df8dcd100baa1b2d037cbfc2ec5；%作图下面附上了图片x(0)=0和x(0)=2的情况显然初始值不同结果不同;%初始值&&&tapan[0,10];%求解时间范围&&&sol=dde23(@yanchi,tau,history,tapan).com/zhidao/pic/item/9e3df8dcd100baa1b2d037cbfc2ec5.jpg" target="_blank" title="点击查看大图" class="ikqb_img_alink"><img class="ikqb_img" src="http；你能给出x(0)的值我可以帮你解&首先编写关于延迟函数的M文件；function&nbsp://h;%给定延迟时间&gt,x;&&plot(sol.x,sol.y).jpg" esrc="&nbsp，这就是为什么需要初始值的情况<a href="http://h.dx=yanchi(t这是一个延迟微分方程；MATLAB可以解这类延迟微分方程.hiphotos.baidu.com/zhidao/wh%3D450%2C600/sign=d6b67c52b31bbc034af682/9e3df8dcd100baa1b2d037cbfc2ec5
采纳率：46%
为您推荐：
其他类似问题
微分方程的相关知识
&#xe675;换一换
回答问题，赢新手礼包&#xe6b9;
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。}