查看: 3501|回复: 24

[共享资源] 【分享】关于Matlab求解微分方程的那些事!!!

[复制链接]
发表于 2014-3-8 21:10 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
本帖最后由 牛小贱 于 2014-5-24 14:02 编辑

看到论坛里,有很多人问关于Matlab求解微分方程的问题。所以,鄙人不才,根据自己的个人经验,以及结合相关论坛、博客,将Matlab求解微分方程的相关内容和大家分享一下!!
首先,我们讲一下关于可以求得解析解的微分方程。大家应该都知道,函数dsolve可以解决一部分能够解析求解的微分方程、微分方程组。
关于该函数的调用格式,如下:
y=dsolve(f1,f2,...,fm,'x');
先举个例子,如图:
搜狗截图20140423191008.png

程序如下:
  1. syms t y;
  2. u=exp(-5*t)*cos(2*t-1)+5;
  3. uu=5*diff(u,t,2)+4*diff(u,t)+2*u;
  4. y=dsolve(['D4y+10*D3y+35*D2y+24*Dy= uu '])
复制代码
接下来,介绍一下关于数值解的解法。
大家应该都有这种意识:能够解析求解的微分方程(组)占少数,对于大多数微分方程(组)而言不能得到解析解,这时只能进行数值求解。关于数值分析方法,比较著名的当属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. function dx=jxhdot(t,x)
  7. 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)
yp0y'的初值。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里可以得到toutyout,执行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)。
不足之处,还请大家补充、指正。






图1

图1

图2

图2

图3

图3

图4

图4
6.png

图6

图6

8

8

图7

图7
11.png
10.png

高等应用数学问题的MATLAB 求解.pdf

1.43 MB, 下载次数: 86

Matlab微分方程的解法.pdf

214.63 KB, 下载次数: 56

点评

赞成: 5.0
赞成: 5
感谢楼主分享,很好的资料  发表于 2014-6-4 13:43
赞成: 5
好资料!!!  发表于 2014-3-25 19:59
如果有不足之处,还请大家补充指正!!  发表于 2014-3-19 09:19
楼主太有耐心了!解答的很给力!!  发表于 2014-3-12 11:04

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2014-3-9 19:19 | 显示全部楼层
关于Matlab求解微分方程的那些事!!!讲的真心不错啊!!楼主真的很有耐心啊!32个赞
发表于 2014-5-23 16:28 | 显示全部楼层
LZ第一个例子好像式子与代码就不一致了
可能误值了
 楼主| 发表于 2014-5-23 20:00 | 显示全部楼层
ChaChing 发表于 2014-5-23 16:28
LZ第一个例子好像式子与代码就不一致了
可能误值了

貌似应该没有吧~呵呵。
发表于 2014-5-24 13:42 | 显示全部楼层
本帖最后由 ChaChing 于 2014-5-24 13:43 编辑

图 u=exp(-5*t)*cos(2*t+1)+5
程序u=exp(-5*t)*cos(2*t-1)+5

图 y''''+10y'''+35y''+24y'=5u''+4u'+2u
程序 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'])

点评

赞成: 5.0
赞成: 5
谢谢您提醒,让我发现错误~谢谢,已更正~  发表于 2014-5-24 14:13
誤植更正下即可 @牛小贱  发表于 2014-5-24 14:03
奥,是的。打错了~我还是看的之前那个“正确的代码”,没太注意现在写的~  发表于 2014-5-24 13:56

评分

1

查看全部评分

发表于 2014-6-17 22:36 | 显示全部楼层
谢谢楼主分享!
发表于 2014-12-23 12:45 | 显示全部楼层
学习啦 谢谢
发表于 2014-12-23 16:43 | 显示全部楼层
楼主,高人!
发表于 2014-12-25 16:29 | 显示全部楼层
感谢楼主分享,学习了
发表于 2015-1-5 21:43 | 显示全部楼层
楼猪好帅!
发表于 2015-5-26 11:01 | 显示全部楼层
很全面 感谢楼主分享
发表于 2015-11-11 16:44 | 显示全部楼层
好书啊 啊啊
发表于 2015-11-11 22:21 | 显示全部楼层
学习了,感谢楼主,感谢老师,给我们发群了
发表于 2015-12-30 17:58 | 显示全部楼层
值得关注,很好的指导
发表于 2016-1-23 13:20 | 显示全部楼层
感谢楼主
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表