1、第四講 Matlab 求解微分方程（組）理論介紹： Matlab 求解微分方程（組）命令求解實例： Matlab 求解微分方程（組）實例 實際應用問題通過數學建模所歸納得到的方程，絕大多數都是微分方程，真 正能得到代數方程的機會很少 .另一方面， 能夠求解的微分方程也是十分有限的， 特別是高階方程和偏微分方程（組） .這就要求我們必須研究微分方程（組）的 解法：解析解法和數值解法 .一相關函數、命令及簡介1. 在 Matlab 中，用大寫字母 D 表示導數， Dy 表示 y 關于自變量的一階導數， D2y 表示 y 關于自變量的二階導數，依此類推 .函數 dsolve 用來解決常微分方程 （組

2、）的求解問題，調用格式為：X=dsolve（ eqn 1' , ' eqn2',）函數dsolve用來解符號常微分方程、方程組，如果沒有初始條件，則求出通 解，如果有初始條件，則求出特解 .注意 ， 系統缺省的自變量為 t2函數dsolve求解的是常微分方程的精確解法，也稱為常微分方程的符號解. 但是，有大量的常微分方程雖然從理論上講， 其解是存在的， 但我們卻無法求出 其解析解，此時，我們需要尋求方程的數值解，在求常微分方程數值解方面， MATLAB具有豐富的函數，我們將其統稱為 solver，其一般格式為：T,Y=solver（odefun,tspan,y0）說明：

3、 （1）solver 為命令 ode45、 ode23、 ode113、 ode15s、 ode23s、 ode23t、 ode23tb、 ode15i 之一.odefun是顯示微分方程y' f（t,y）在積分區間tspan t°,tf上從t°到tf用 初始條件y求解（3）如果要獲得微分方程問題在其他指定時間點t0,t1,t2丄1上的解，則令tspan出赳上丄tf （要求是單調的）.（4）因為沒有一種算法可以有效的解決所有的 ODE問題，為此，Matlab提供 了多種求解器solver，對于不同的ODE問題，采用不同的solver.表1 Matlab中文本文件讀寫函

4、數求解器ODE類型特點說明ode45非剛性單步算法：4、5階Runge-Kutta方程；累計截斷誤差（x）3大部分場合的首選算法ode23非剛性單步算法：2、3階Runge-Kutta方程；累計截斷誤差（x）3使用于精度較低的情形ode113非剛性多步法：Adams算法；高低精度可達 10 3 10 6計算時間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法：Gear's反向數值微分；精度中等若ode45失效時，可嘗試使用ode23s剛性單步法：2階Rosebrock算法；低精度當精度較低時，計算時間比ode15s短ode23tb剛性梯形算法；低精度當精

5、度較低時，計算時間比ode15s短說明：ode23、ode45是極其常用的用來求解非剛性的標準形式的一階微分方程（組）的初值問題的解的Matlab常用程序，其中：ode23采用龍格-庫塔2階算法，用3階公式作誤差估計來調節步長，具有低 等的精度ode45則采用龍格-庫塔4階算法，用5階公式作誤差估計來調節步長，具有 中等的精度3.在matlab命令窗口、程序或函數中創建局部函數時，可用內聯函數inline， inline函數形式相當于編寫M函數文件，但不需編寫 M-文件就可以描述出某種 數學關系 調用inline函數，只能由一個 matlab表達式組成，并且只能返回一個 變量，不允許u,v這種

6、向量形式.因而，任何要求邏輯運算或乘法運算以求得最 終結果的場合，都不能應用inline函數，inline函數的一般形式為：FunctionName=inline（函數內容'所有自變量列表'例如：（求解F（x）=x2*cos（a*x）-b ,a,b是標量；x是向量）在命令窗口輸入：Fofx=inline(x .A2*cos(a(*X) ,x' , ' a , ' b');g= Fofx(pi/3 pi/3.5,4,1)系統輸出為：g=-1.5483 -1.7259注意：由于使用內聯對象函數inline不需要另外建立m文件，所有使用比較 方便，另外

7、在使用ode45函數的時候，定義函數往往需要編輯一個 m文件來單 獨定義，這樣不便于管理文件，這里可以使用 inline來定義函數.二實例介紹1. 幾個可以直接用Matlab求微分方程精確解的實例例1求解微分方程y 2xy xe “程序：syms x y; y=dsolve( Dy+2*x*y=x*exp(-xA2)例2求微分方程xy' y ex 0在初始條件y(1) 2e下的特解并畫出解函數 的圖形.程序：syms x y; y=dsolve( x*Dy+y-exp(1)=0 'y(1)=2*exp(1) '''ezplot(y)5xdxtdt魚dte

8、例3求解微分方程組在初始條件x|t o 1, y |t o 0下的特解3y0并畫出解函數的圖形.程序：syms x y tx,y=dsolve(,Dx+5*x+y=exp(t)',Dy-x-3*y=0,x(0)=1,y(0)=0,t,)simple(x);simple(y)ezplot(x,y,0,1.3);axis auto2. 用 ode23 ode45等求解非剛性標準形式的一階微分方程(組)的初值問題 的數值解(近似解)史 2y 2x2 2x例4求解微分方程初值問題dx 2y 2x2x的數值解，求解范圍為區y(0) 1間0,0.5.程序：fun=inline(,-2*y+2*xA

9、2+2*x,x,y,);x,y=ode23(fu n,0,0.5,1);plot(x,y,'o-')歡迎下載13例5求解微分方程d2ydt2dx2dtdx1dtX2,捲(0)17(1 x2 )x2 x!, x2(0)0(1 y2)dy y 0, y(0) l,y'(0)0的解，并畫出dt解的圖形.分析：這是一個二階非線性方程，我們可以通過變換，將二階方程化為一階 方程組求解令“ y，X2篇7，則編寫M-文件vdp.mfun cti on fy=vdp(t,x)fy=x(2);7*(1-x(1)A2)*x (2) -x(1);end在Matlab命令窗口編寫程序y0=1;

10、0t,x=ode45(vdp,0,40,y0);或t,x=ode45('vdp',0,40,y0);y=x(:,1);dy=x(:,2);plot(t,y,t,dy)練習與思考：M-文件vdp.m改寫成inline函數程序？3用Euler折線法求解Euler折線法求解的基本思想是將微分方程初值問題dydxf (x, y)y(x) y。化成一個代數(差分)方程，主要步驟是用差商y(x h) y(x)替代微商 魚，于是 hdxy(xk hh y(xk)f(xk,y(xk)hy。y(x。)記 Xk 1Xkh, yky(Xk),從而 yk 1 yX h),于是yo y(x。)，Xk 1

11、 Xk h, k 0,1,2丄,n 1 yk 1 yk hf(Xk,yQ.例6用Euler折線法求解微分方程初值問題dydxy(0)2x2y1的數值解(步長h取0.4),求解范圍為區間0,2.分析：本冋題的差分方程為X0 0, y0 1,h 0.4Xk 1 Xk h, k 0,1,2,L , n 1 yk 1yk hf (Xk, yk).程序：>> clear>> f=sym('y+2*X/yA2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;&g

12、t;> y=1;>> szj=x,y;% 數值解>> for i=1: n-1y=y+h*subs(f,'x','y',x,y);%subs，替換函數 x=x+h; szj=szj;x,y;end>>szj>> plot(szj(:,1),szj(:,2)說明：替換函數subs例如：輸入subs(a+b,a,4)意思就是把a用4替換掉, 返回 4+b，也可以替換多個變量，例如：subs(cos(a)+sin(b),a,b,sym('alpha'),2)分別用字符alpha替換a和2替換b,返回

13、cos(alpha)+sin(2)特別說明：本問題可進一步利用四階 Runge-Kutta法求解，Euler折線法實際上就是一階Runge-Kutta法，Runge-Kutta法的迭代公式為y。y(x。)，Xk 1 Xk h,L1f (Xk, yk)1f (xkh,yk如,f (xkh2,ykf(xkh,ykhLs).>> cleark 0,1,2丄，n 1L2L4相應的Matlab程序為:ykiyk (Li 2L2 2L3 L4),6>> f=sym('y+2*x/yA2');>> a=0;>> b=2;>> h=0

14、.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=x,y;% 數值解>> for i=1: n-1l1= subs(f, 'x','y',x,y);替換函數I2=subs(f, 'x','y',x+h/2,y+l1*h/2);I3=subs(f, 'x','y',x+h/2,y+l2*h/2);l4=subs(f, 'x','y',x+h,y+l3*h);y=y+h*(l1+2*l

15、2+2*l3+l4)/6; x=x+h;szj=szj;x,y;end>>szj>> plot(szj(:,1),szj(:,2)練習與思考：(1)ode45求解問題并比較差異.(2) 利用Matlab求微分方程y2yy"0的解.HQI(3) 求解微分方程 y 2(1 y )y y 0,0 x 30,y(0)1,y'(0)0 的特解., QHII(4) 利用Matlab求微分方程初值問題(1 x )y 2xy , y |x 0 1, y |x 0 3的解.提醒：盡可能多的考慮解法三微分方程轉換為一階顯式微分方程組Matlab微分方程解算器只能求解標準形

16、式的一階顯式微分方程(組)問題， 因此在使用ODE解算器之前，我們需要做的第一步，也是最重要的一步就是借 助狀態變量將微分方程(組)化成 Matlab可接受的標準形式.當然，如果ODEs 由一個或多個高階微分方程給出，則我們應先將它變換成一階顯式常微分方程組 下面我們以兩個高階微分方程組構成的ODEs為例介紹如何將它變換成一個一階顯式微分方程組.Step 1將微分方程的最高階變量移到等式左邊，其它移到右邊，并按階次從低到高排列.形式為：(m)'''(m 1)'''(n 1)x( ) f (t,x,x,x 丄,x( ),y, y ,y 丄，y()y

17、() g(t,x,x,x 丄,x( ),y,y ,y 丄,y()Step 2為每一階微分式選擇狀態變量，最高階除外'''(m 1)y(n 1)X1 X,X2 X,X3 x 丄，Xm xXm 1 y,Xm 2 y ,Xm 3 、，L ,Xm n注意：ODEs中所有是因變量的最高階次之和就是需要的狀態變量的個數, 最高階的微分式不需要給它狀態變量Step 3根據選用的狀態變量，寫出所有狀態變量的一階微分表達式x； X2,X2 X3,X3 X4丄,Xmf(t,X1,X2,X3,L ,xg(t,X1,X2,x 3,l練習與思考：(1)求解微分方程組IIIx 2y x(x )3r

18、i*(x )32y 2x yy3A1/82.45, x(0)1.2,其中2,(x 丫 y2，i 、.(x )2 y2,*1y(0) 0, x(0) 0, y'(0)1.049355751(2) 求解隱式微分方程組ninx 2yx 2yn ii nix y 3x y xy y 5提示：使用符號計算函數solve求x", y"，然后利用求解微分方程的方法四.偏微分方程解法Matlab提供了兩種方法解決PDE問題，一是使用pdepe函數，它可以求解 一般的PDEs具有較大的通用性，但只支持命令形式調用；二是使用 PDE工具 箱，可以求解特殊PDE問題，PDEtoll有較大

19、的局限性，比如只能求解二階PDE 問題，并且不能解決片微分方程組，但是它提供了GUI界面，從復雜的編程中解脫出來，同時還可以通過 File>SaveAs直接生成M代碼.1.一般偏微分方程(組)的求解(1) Matlab提供的pdepe函數，可以直接求解一般偏微分方程(組)，它的調 用格式為：sol=pdepe(m,pdef un, pdeic,pdebc,x,t)pdefun是PDE的問題描述函數，它必須換成標準形式：c(x,tuu x mxmf (Xtu,) MX't'U,)x txxx這樣，PDE就可以編寫入口函數：c,f,s=pdefun(x,t,u,du)，m,x

20、,t對應于式中相關 參數，du是u的一階導數，由給定的輸入變量可表示出c,f,s這三個函數.pdebc是PDE的邊界條件描述函數，它必須化為形式：p(x,t,u) q(x,t,u).* f (x,t,uu)0x于是邊值條件可以編寫函數描述為：pa,qa,pb,qb=pdebc(x,t,u,du)其中a表示下 邊界，b表示上邊界.pdeic是PDE的初值條件，必須化為形式:u(x,tc) U0，故可以使用函數描述為：uO=pdeic(x)sol是一個三維數組,sol(:,:,i)表示Ui的解，換句話說，Uk對應x(i)和t(j)時的解為sol(i,j,k)，通過sol，我們可以使用pdeval函

21、數直接計算某個點的函數值.(2) 實例說明求解偏微分uu.tu2t20.024 F(u, u2)x20.17 F(u, u2)x其中，F(x) e5.73x e 11.46x且滿足初始條件u1(x,0) 1,u2(x,0)0及邊界條件巴(0,t) 0,業(0心 0M(1,t) 1,且(1,t) 0xx解：(1)對照給出的偏微分方程和pdepe函數求解的標準形式，原方程改寫為可見m 0,c1u1* I1 . t u20.0240.17u2F(W 比)F(u1 氏)0.0245F(u1 u?)x,s0.17-u2F(u1 吐)xx%目標PDE函數function c,f,s=pdefun(x,t,

22、u,du)c=1；1；f=0.024*du(1);0.17*du (2);temp=u(1)-u(2);s=-1;1.*(exp(5.73*temp)-exp(-11.46*temp) end(2)邊界條件改寫為:下邊界上邊界u2u1 10%邊界條件函數fun cti on pa,qa,pb,qb=pdebc(xa,ua,xb,ub,t)pa=0;ua(2);qa=1;0;pb=ub(1)-1;0;qb=0;1;end(3) 初值條件改寫為：:1%初值條件函數fun cti on u0=pdeic(x)u0=1;0;end編寫主調函數clcx=0:0.05:1;t=0:0.05:2;m=0;s

23、ol=pdepe(m,pdefu n,pdeic,pdebc,x,t);subplot(2,1,1) surf(x,t,sol(:,:,1)subplot(2,1,2) surf(x,t,sol(:,:,2)練習與思考: This example illustrates the straightforward formulation, computation, and plott ing of the soluti on of a sin gle PDE.2專+t x xThis equation holds on an interval 0 x 1 for times t 0. The PD

24、E satisfies theinitial condition u(x,0) sin x and boundary conditionsu(0,t)0; et (1,t)0x2.PDEtool求解偏微分方程(1) PDEtool (GUI)求解偏微分方程的一般步驟在Matlab命令窗口輸入 pdetool,回車，PDE工具箱的圖形用戶界面(GUI) 系統就啟動了 從定義一個偏微分方程問題到完成解偏微分方程的定解，整個過 程大致可以分為六個階段Step 1 “Draw模式”繪制平面有界區域，通過公式把Matlab系統提供的實體模型：矩形、圓、橢圓和多邊形，組合起來，生成需要的平面區域Step 2 “Boundary模式”定義邊界，聲明不同邊界段的邊界條件.Step 3 “ PDE模式”定義偏微分方