【原创】关于Matlab求解微分方程的那些事！！！

Nicky_小牛

根据自己的个人经验，以及结合相关论坛、博客，将Matlab求解微分方程的相关内容和大家分享一下！！（PS：该贴来自中国振动联盟论坛中我发表的原创帖子）。

首先，我们讲一下关于可以求得解析解的微分方程。大家应该都知道，函数dsolve可以解决一部分能够解析求解的微分方程、微分方程组。

关于该函数的调用格式，如下：

y=dsolve(f1,f2,...,fm,'x');

先举个例子，如图：

程序如下：

1. syms t;
   
2. u=exp(-5*t)*cos(2*t-1)+5;
   
3. uu=5*diff(u,t,2)+4*diff(u,t)+2*u;
   
4. syms t y;
   
5. 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'])

复制代码

接下来，介绍一下关于数值解的解法。

大家应该都有这种意识：能够解析求解的微分方程（组）占少数，对于大多数微分方程（组）而言不能得到解析解，这时只能进行数值求解。关于数值分析方法，比较著名的当属Eular法、 Runge-Kutta了！在matlab中内置的ode求解器可以实现不同求解方法的相同格式的调用，而不必太关心matlab究竟是用什么算法完成的。在这里，我主要和大家分享关于ode45等求解器的具体使用方法，其他的求解器的使用有相似之处，在此不再赘述。

先不介绍其具体调用格式，先来个例子，看看他的庐山真面目：

求解方程组

Dx=y+x(1-x^2-y^2);

Dy=-x+y*(1-x^2-y^2)

初值x=0.1;y=0.2;

程序：

1. function weifen1
   
2. clear;clc
   
3. x0=[0.1;0.2];
   
4. [t,x]=ode45(@jxhdot,[0,100],x0);
   
5. plot(x(:,1),x(:,2))
   
6. 
   
7. function dx=jxhdot(t,x)
   
8. dx=[x(2)+x(1).*(1-x(1).^2-x(2).^2); -x(1)+x(2).*(1-x(1).^2-x(2).^2)];

复制代码

结果如图1所示：

看完例子，我说明一下最常用的ode45调用方式，和相应的函数文件定义格式。

[t,x]=ode45(odefun,tspan,x0);

其中，Fun就是导函数，tspan为求解的时间区间（或时间序列，如果采用时间序列，则必须单调），x0为初值。

这时，函数文件可以采用如下方式定义：

function dx=odefun(t,x)

接着， （1）刚性方程的求解：刚性方程就是指各个自变量的变化率差异很大，会造成通常的求解方法失效。下面是matlab中自带的一个例子，使用ode15s求解，如果用ode45求解就会出现错误。

1. function myode15study
   
2. [T,Y] = ode15s(@vdp1000,[0 3000],[2 0]);
   
3. plot(T,Y(:,1),'-o')
   
4. figure;plot(Y(:,1),Y(:,2))
   
5. function dy = vdp1000(T,y)
   
6. dy = zeros(2,1);   
   
7. dy(1) = y(2);
   
8. dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);
   

复制代码

结果如图2和图3.

（2）高阶微分方程的求解！通常的方法是进行变量替换，将原方程降阶，转换成更多变量的一阶方程组进行求解。在这个例子里我们求解一个动力学系统里最常见的一个运动方程。

程序：

1. function myhighoder
   
2. clear;clc
   
3. x0=zeros(6,1);
   
4. [t,x]=ode45(@myhigh,[0,100],x0);
   
5. plot(t,x(:,1))
   
6. function dx=myhigh(t,x)
   
7. f=[sin(t);0;0];
   
8. M=eye(3);
   
9. C=eye(3)*0.1;
   
10. K=eye(3)-0.5*diag(ones(2,1),1)-0.5*diag(ones(2,1),-1);
    
11. dx=[x(4:6);inv(M)*(f-C*x(4:6)-K*x(1:3))];

复制代码

结果如图4.

（3）延迟微分方程：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求解延迟微分方程的例子：

1. function mydde23study
   
2. %   The differential equations
   
3. %
   
4. %        y'_1(t) = y_1(t-1)
   
5. %        y'_2(t) = y_1(t-1)+y_2(t-0.2)
   
6. %        y'_3(t) = y_2(t)
   
7. %
   
8. %   are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for
   
9. %   t <= 0.
   
10. clear;clc
    
11. lags=[1,0.2];
    
12. history=[1;1;1];
    
13. tspan=[0,5];
    
14. sol = dde23(@myddefun,lags,history,tspan)
    
15. plot(sol.x,sol.y)
    
16. function dy = myddefun(t,y,Z)
    
17. dy=[Z(1,1);Z(1)+Z(2,2); y(2) ];
    

复制代码

结果如图5所示。

（4）隐式微分方程：求解此类方程需要使用ode15i函数。调用格式：[T,Y] = ode15i(odefun,tspan,y0,yp0)

yp0为y'的初值。odefun的格式如下  dy =odefun(t,y,yp),yp表示y'，而方程中应该使得f(t,y,y')=0

1. function myodeIMP
   
2. %   The problem is
   
3. %
   
4. %         y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)
   
5. %         y(2)' =  0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
   
6. %         y(3)' =  3e7*y(2)^2
   
7. %
   
8. %   It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3) = 0
   
9. %   to steady state.
   
10. clear;clc
    
11. y0=[1;0;0];
    
12. fixed_y0=[1;1;1];
    
13. yp0=[0 0 0];
    
14. fixed_yp0=[];
    
15. [y0mod,yp0mod]=decic(@myodefunimp,0,y0,fixed_y0,yp0,fixed_yp0);
    
16. tspan=[0, logspace(-6,6)];
    
17. [t,y] = ode15i(@myodefunimp,tspan,y0mod,yp0mod);
    
18. y(:,2)=1e4*y(:,2);
    
19. semilogx(t,y)
    
20. function res=myodefunimp(t,y,yp)
    
21. res=[ -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 ];

复制代码

运行结果如图6.
在讲了这么多之后，我想大家也看烦了！所以，我结合相关论坛、帖子总结一个新的求解ode的方法——那就是使用simulink的积分器求解！！

还是咱们举的第一个例子：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的计算公式。

在Simulink中建立模块框图，如图7所示。

运行这个仿真，scope中可以看到两个变量的时程如图8.在WorkSpace里可以得到tout和yout，执行plot(yout(:,1),yout(:,1))得到与ode45求解相似的结果！！！

除此之外，使用simulink还可以求解的一类延迟微分方程。（PS：transportDelay,就可以实现对于信号的延迟！！）

实例：x'(t)=A1*x(t-t1)+A2*x'(t-t2)+B*u(t)

t1=0.15;t2=0.5

A1=[-12     3   -3]      A2=[0.02    0     0]    B=[0]

   [106  -116   62]            [0    0.03     0]         [1]

   [207  -207  113]            [0       0  0.04]        [2]   

Simulink中建立的模块框图和仿真结果图，如图9和10所示.

最后，给大家一些matlab求解微分方程的资料！见附件（本人强烈推荐一本书      薛定宇编写的《高等应用数学问题的MATLAB求解》，很不错！链接地址为http://ishare.iask.sina.com.cn/f/7453073.html）。

不足之处，还请大家补充、指正。

来自群组: 重庆大学

小贝

感谢分享……学习学习，貌似很不错！

小蝶

楼主真是太有耐心了！谢谢LZ无私整理分享！！！！

小蝶

自己正在学习Matlab编程基础，是参考LZ发过的帖子：matlab经典算法的程序+实用程序300例+MATlAB中文编程基础！ 这份资料感觉确实不错，等遇到问题，再向LZ请教哈~先谢谢了。

Nicky_小牛

小蝶 发表于 2014-5-22 11:36
自己正在学习Matlab编程基础，是参考LZ发过的帖子：matlab经典算法的程序+实用程序300例+MATlAB中文编程基 ...

恩，好的，初学者要先把Matlab基本编程知识掌握了，有什么问题再讨论~

雨沐枫

正好需要，感谢分享

Nicky_小牛

雨沐枫 发表于 2014-6-18 16:05
正好需要，感谢分享

不客气，有什么问题，可以讨论~

caoyaocui

请教一下，偏微分方程用matlab怎么解决？

Nicky_小牛

caoyaocui 发表于 2014-8-6 16:57
请教一下，偏微分方程用matlab怎么解决？

这个嘛，希望你自己好好看一下我发的帖子，其次，要自己到网上找一些资料，自我先学习一下，不建议这样笼统地问……

caoyaocui

Nicky_小牛 发表于 2014-8-6 21:43
这个嘛，希望你自己好好看一下我发的帖子，其次，要自己到网上找一些资料，自我先学习一下，不建议这样笼 ...

你的帖子看了，不能求解偏微分方程:(

caoyaocui

Nicky_小牛 发表于 2014-8-6 21:43
这个嘛，希望你自己好好看一下我发的帖子，其次，要自己到网上找一些资料，自我先学习一下，不建议这样笼 ...

你发的电子书怎么是从第十章开始的？

追风少年-阿杰

书不全，怎么搞的呢 版主

liuweiji2311

好好学习