Matlab求解微分方程組及偏微分方程組

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

2、qn2 ,）函數(shù)dsolve用來解符號常微分方程、方程組，如果沒有初始條件，則求出通解，如果 有初始條件，則求出特解注意，系統(tǒng)缺省的自變量為 t2函數(shù)dsolve求解的是常微分方程的精確解法，也稱為常微分方程的符號解但是，有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無法求出其解析解，此時，我們需要尋求方程的數(shù)值解，在求常微分方程數(shù)值解方面，MATLAB具有豐富的函數(shù)，我們將其統(tǒng)稱為solver，其一般格式為：T,Y=solver（odefu n,tspa n,yO）說明：（1）solver 為命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、od

3、e23tb、ode15i 之一.（2）odefun是顯示微分方程在積分區(qū)間tspan上從到用初始條件求解.（3）如果要獲得微分方程問題在其他指定時間點上的解，則令tspan（要求是單調(diào)的）.（4）因為沒有一種算法可以有效的解決所有的ODE問題，為此，Matlab提供了多種求解器solver，對于不同的 ODE問題，采用不同的 solver.表1 Matlab中文本文件讀寫函數(shù)求解器ODE類型特點說明ode45非剛性單步算法：4、5階Runge-Kutta方程；累計截斷誤差（x）3大部分場合的首選算法ode23非剛性單步算法：2、3階Runge-Kutta方程；累計截斷誤差(x)3使用于精度較低

4、的情形ode113非剛性多步法：Adams算法；高低精度可達 10 3 10 6計算時間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法：Gears反向數(shù)值微分；精度中等若ode45失效時，可嘗試使用ode23s剛性單步法：2階Rosebrock算法；低精度當精度較低時，計算時間比ode15s短ode23tb剛性梯形算法；低精度當精度較低時，計算時間比ode15s短說明：ode23、ode45是極其常用的用來求解非剛性的標準形式的一階微分方程(組)的初值問題的解的 Matlab常用程序，其中：ode23采用龍格-庫塔2階算法，用3階公式作誤差估計來調(diào)節(jié)步長，具有

5、低等的精度ode45則采用龍格-庫塔4階算法，用5階公式作誤差估計來調(diào)節(jié)步長，具有中等的精度3.在matlab命令窗口、程序或函數(shù)中創(chuàng)建局部函數(shù)時，可用內(nèi)聯(lián)函數(shù)inline , inline函數(shù)形式相當于編寫 M函數(shù)文件，但不需編寫 M-文件就可以描述出某種數(shù)學(xué)關(guān)系 調(diào)用inline 函數(shù)，只能由一個 matlab表達式組成，并且只能返回一個變量，不允許 u,v這種向量形式. 因而，任何要求邏輯運算或乘法運算以求得最終結(jié)果的場合，都不能應(yīng)用inline函數(shù)，inline函數(shù)的一般形式為：FunctionName=inline(函數(shù)內(nèi)容所有自變量列表)例如：(求解F(x)=xA2*cos(a*x

6、)-b,a,b 是標量；x是向量)在命令窗口輸入：Fofx=inline( x .A2*cos(0*x), x , a , b);g= Fofx(pi/3 pi/3.5,4,1)系統(tǒng)輸出為：g=-1.5483 -1.7259注意：由于使用內(nèi)聯(lián)對象函數(shù) inline不需要另外建立 m文件，所有使用比較方便，另 外在使用ode45函數(shù)的時候，定義函數(shù)往往需要編輯一個m文件來單獨定義，這樣不便于管理文件，這里可以使用 inline來定義函數(shù).實例介紹1幾個可以直接用Matlab求微分方程精確解的實例例1求解微分方程程序：syms x y; y=dsolve( Dy+2*x*y=x*expX2) ,

7、x)例2求微分方程在初始條件下的特解并畫出解函數(shù)的圖形程序：syms x y; y=dsolve(x*Dy-exp(1)=0 , y(1)=2*exp(1) , x );ezplot(y)例3求解微分方程組在初始條件下的特解并畫出解函數(shù)的圖形程序：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等求解非剛性標準形式的一階微分方程(組)的初值問題的數(shù)值解(近似解)例4求解微分方程初值問題的數(shù)值解，求解

8、范圍為區(qū)間0,0.5.程序：fun=i nli ne(-2*y+2*xA2+2*x,x,y);x,y=ode23(fu n,0,0.5,1);plot(x,y,o-)例5求解微分方程的解，并畫出解的圖形.分析：這是一個二階非線性方程，我們可以通過變換，將二階方程化為一階方程組求解.令，則編寫M-文件vdp.mfunction fy=vdp(t,x)fy=x(2);7*(1-x(1)A2)*x (2)-x(1);end在Matlab命令窗口編寫程序y0=i;0t,x=ode45(vdp,0,40,y0);或t,x=ode45(vdp,0,40,y0);y=x(:,1);dy=x(:,2);plo

9、t（t,y,t,dy）練習(xí)與思考：M-文件vdp.m改寫成inline函數(shù)程序？3用Euler折線法求解Euler折線法求解的基本思想是將微分方程初值問題化成一個代數(shù)（差分）方程，主要步驟是用差商替代微商，于是記從而于是例6用Euler折線法求解微分方程初值問題的數(shù)值解（步長取0.4），求解范圍為區(qū)間0,2.分析：本問題的差分方程為程序： clear f=sym(y+2*x/yA2); a=0; b=2; h=0.4; n=(b_a)/h+1; x=0; y=1; szj=x,y;% 數(shù)值解 for i=1: n-1y=y+h*subs(f,x,y,x,y);%subs，替換函數(shù)x=x+h;s

10、zj=szj;x,y;endszj plot(szj(:,1),szj(:,2)說明：替換函數(shù) subs例如：輸入subs(a+b,a,4)意思就是把a用4替換掉，返回 4+b,也可以替換多個變量，例如：subs(cos(a)+sin(b),a,b,sym(alpha),2)分別用字符 alpha替換a 和 2 替換 b，返回 cos(alpha)+sin(2)特別說明：本問題可進一步利用四階Runge-Kutta法求解，Euler折線法實際上就是一階Runge-Kutta法，Runge-Kutta法的迭代公式為相應(yīng)的Matlab程序為： clear f=sym(y+2*x/yA2); a=0

11、; b=2; h=0.4; n=(b_a)/h+1; x=0; y=1; szj=x,y;% 數(shù)值解 for i=1: n-1I仁 subs(f, x,y,x,y);替換函數(shù)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*l2+2*l3+l4)/6;x=x+h;szj=szj;x,y;endszj plot(szj(:,1),szj(:,2)練習(xí)與思考：(1)ode45求解問題并比較差異.利用Matlab求微分方程的解.(3)求解微分方程的特

12、解.利用Matlab求微分方程初值問題的解.提醒：盡可能多的考慮解法三微分方程轉(zhuǎn)換為一階顯式微分方程組Matlab微分方程解算器只能求解標準形式的一階顯式微分方程（組）問題，因此在使用ODE解算器之前，我們需要做的第一步，也是最重要的一步就是借助狀態(tài)變量將微分方 程（組）化成Matlab可接受的標準形式當然，如果ODEs由一個或多個高階微分方程給出， 則我們應(yīng)先將它變換成一階顯式常微分方程組下面我們以兩個高階微分方程組構(gòu)成的ODEs為例介紹如何將它變換成一個一階顯式微分方程組Step 1將微分方程的最高階變量移到等式左邊，其它移到右邊，并按階次從低到高排列.形式為：Step 2為每一階微分式選

13、擇狀態(tài)變量，最高階除外注意：ODEs中所有是因變量的最高階次之和就是需要的狀態(tài)變量的個數(shù)，最高階的微分式不需要給它狀態(tài)變量Step 3根據(jù)選用的狀態(tài)變量，寫出所有狀態(tài)變量的一階微分表達式練習(xí)與思考：（1）求解微分方程組其中（2）求解隱式微分方程組提示：使用符號計算函數(shù)solve求，然后利用求解微分方程的方法四.偏微分方程解法Matlab提供了兩種方法解決 PDE問題，一是使用pdepe函數(shù)，它可以求解一般的 PDEs, 具有較大的通用性，但只支持命令形式調(diào)用；二是使用PDE工具箱，可以求解特殊 PDE問題，PDEtoll有較大的局限性，比如只能求解二階PDE問題，并且不能解決片微分方程組，但是

14、它提供了 GUI界面，從復(fù)雜的編程中解脫出來，同時還可以通過File Save As直接生成M代碼.1. 一般偏微分方程（組）的求解(1) Matlab提供的pdepe函數(shù)，可以直接求解一般偏微分方程(組)，它的調(diào)用格式為：sol=pdepe(m,pdefu n, pdeic,pdebc,x,t)pdefun是PDE的問題描述函數(shù)，它必須換成標準形式：這樣，PDE就可以編寫入口函數(shù)：c,f,s=pdefun(x,t,u,du) ， m,x,t對應(yīng)于式中相關(guān)參數(shù)，du是u的一階導(dǎo)數(shù)，由給定的輸入變量可表示出c,f,s這三個函數(shù).pdebc是PDE的邊界條件描述函數(shù)，它必須化為形式：于是邊值條件可

15、以編寫函數(shù)描述為：pa,qa,pb,qb=pdebc(x,t,u,du),其中a表示下邊界，b表示上邊界.pdeic是PDE的初值條件，必須化為形式：，故可以使用函數(shù)描述為： u0=pdeic(x) sol是一個三維數(shù)組，sol(:,:,i)表示的解，換句話說，對應(yīng)x(i)和t(j)時的解為sol(i,j,k)，通過sol，我們可以使用 pdeval函數(shù)直接計算某個點的函數(shù)值.(2) 實例說明求解偏微分其中，且滿足初始條件及邊界條件解：(1)對照給出的偏微分方程和pdepe函數(shù)求解的標準形式，原方程改寫為可見%目標PDE函數(shù)function c,f,s=pdefu n( x,t,u,du)c=

16、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)邊界條件改寫為：下邊界上邊界%邊界條件函數(shù)7function 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) 初值條件改寫為：%初值條件函數(shù)function uO=pdeic(x)u0=1;0;end(4) 編寫主調(diào)函數(shù)clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,pdefu

17、 n, pdeic,pdebc,x,t);subplot(2,1,1) surf(x,t,sol(:,:,1)subplot(2,1,2) surf(x,t,sol(:,:,2)練習(xí)與思考：This example illustrates the straightforward formulati on, putati on, and plott ing ofthe soluti on of a sin gle PDE.This equation holds on an intervalfor times . The PDE satisfies the initial condition an

18、dboun dary con diti ons2. PDEtool求解偏微分方程（1）PDEtool （GUI）求解偏微分方程的一般步驟在Matlab命令窗口輸入pdetool，回車，PDE工具箱的圖形用戶界面（GUI）系統(tǒng)就啟動 了 .從定義一個偏微分方程問題到完成解偏微分方程的定解，整個過程大致可以分為六個階Step 1“ D模式”繪制平面有界區(qū)域，通過公式把Matlab系統(tǒng)提供的實體模型：矩形、圓、橢圓和多邊形，組合起來，生成需要的平面區(qū)域Step 2“ Bound模式”定義邊界，聲明不同邊界段的邊界條件.Step 3“ P模式”定義偏微分方程，確定方程類型和方程系數(shù) c,a,f,d,根據(jù)具體情況，還可以在不同子區(qū)域聲明不同系數(shù).Step 4 “ Mesh莫式”網(wǎng)格化區(qū)域，可以控制自動生成網(wǎng)格的參數(shù)，對生成的網(wǎng)格進行多次細化，使網(wǎng)格分割更細更合理Step 5 “ Solve模式”解偏微分方程，對于橢圓型方程可以激活并控制非線