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

1、精選優(yōu)質(zhì)文檔-傾情為你奉上第四講 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=ds

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

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

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

5、格-庫塔4階算法，用5階公式作誤差估計來調(diào)節(jié)步長，具有中等的精度.3在matlab命令窗口、程序或函數(shù)中創(chuàng)建局部函數(shù)時，可用內(nèi)聯(lián)函數(shù)inline，inline函數(shù)形式相當(dāng)于編寫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)=x2*cos(a*x)-b ,a,b是標量；x是向量 ）在命令

6、窗口輸入:Fofx=inline(x .2*cos(a*x)-b , 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*exp(-x2),x)例2 求微分方程在初始條件下的特解并畫出解函數(shù)

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

8、題的數(shù)值解（近似解）例 4 求解微分方程初值問題的數(shù)值解，求解范圍為區(qū)間0,0.5.程序：fun=inline('-2*y+2*x2+2*x','x','y');x,y=ode23(fun,0,0.5,1);plot(x,y,'o-')例 5 求解微分方程的解，并畫出解的圖形.分析：這是一個二階非線性方程，我們可以通過變換，將二階方程化為一階方程組求解.令，則編寫M-文件vdp.mfunction fy=vdp(t,x)fy=x(2);7*(1-x(1)2)*x(2)-x(1);end在Matlab命令窗口編寫程序y0=1;0t,

9、x=ode45(vdp,0,40,y0);或t,x=ode45('vdp',0,40,y0);y=x(:,1);dy=x(:,2);plot(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/y2');&

10、gt;> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=x,y;%數(shù)值解>> for i=1:n-1 y=y+h*subs(f,'x','y',x,y);%subs，替換函數(shù) x=x+h; szj=szj;x,y; end>>szj>> plot(szj(:,1),szj(:,2)說明：替換函數(shù)subs例如：輸入subs(a+b,a,4) 意思就是把a用4替換掉，返回 4+b，也可以替

11、換多個變量，例如：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/y2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>&

12、gt; x=0;>> y=1;>> szj=x,y;%數(shù)值解>> for i=1:n-1 l1=subs(f, 'x','y',x,y);替換函數(shù) l2=subs(f, 'x','y',x+h/2,y+l1*h/2); l3=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=sz

13、j;x,y; end>>szj>> plot(szj(:,1),szj(:,2)練習(xí)與思考：(1)ode45求解問題并比較差異.(2)利用Matlab求微分方程的解.(3)求解微分方程的特解.(4)利用Matlab求微分方程初值問題的解.提醒：盡可能多的考慮解法三微分方程轉(zhuǎn)換為一階顯式微分方程組 Matlab微分方程解算器只能求解標準形式的一階顯式微分方程（組）問題，因此在使用ODE解算器之前，我們需要做的第一步，也是最重要的一步就是借助狀態(tài)變量將微分方程（組）化成Matlab可接受的標準形式.當(dāng)然，如果ODEs由一個或多個高階微分方程給出，則我們應(yīng)先將它變換成一階顯式

14、常微分方程組.下面我們以兩個高階微分方程組構(gòu)成的ODEs為例介紹如何將它變換成一個一階顯式微分方程組.Step 1 將微分方程的最高階變量移到等式左邊，其它移到右邊，并按階次從低到高排列.形式為：Step 2 為每一階微分式選擇狀態(tài)變量，最高階除外注意：ODEs中所有是因變量的最高階次之和就是需要的狀態(tài)變量的個數(shù)，最高階的微分式不需要給它狀態(tài)變量.Step 3 根據(jù)選用的狀態(tài)變量，寫出所有狀態(tài)變量的一階微分表達式練習(xí)與思考：(1)求解微分方程組其中(2)求解隱式微分方程組提示：使用符號計算函數(shù)solve求，然后利用求解微分方程的方法四偏微分方程解法Matlab提供了兩種方法解決PDE問題，一是

15、使用pdepe函數(shù)，它可以求解一般的PDEs,具有較大的通用性，但只支持命令形式調(diào)用；二是使用PDE工具箱，可以求解特殊PDE問題，PDEtoll有較大的局限性，比如只能求解二階PDE問題，并且不能解決片微分方程組，但是它提供了GUI界面，從復(fù)雜的編程中解脫出來，同時還可以通過File>Save As直接生成M代碼.1.一般偏微分方程（組）的求解(1)Matlab提供的pdepe函數(shù)，可以直接求解一般偏微分方程（組），它的調(diào)用格式為：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)pdefun是PDE的問題描述函數(shù)，它必須換成標準形式：這樣，PDE就可以編寫入口函數(shù)

16、：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ù)，它必須化為形式：于是邊值條件可以編寫函數(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)實例說明

17、求解偏微分其中，且滿足初始條件及邊界條件解：(1)對照給出的偏微分方程和pdepe函數(shù)求解的標準形式，原方程改寫為可見%目標PDE函數(shù)function c,f,s=pdefun(x,t,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)邊界條件改寫為：下邊界上邊界%邊界條件函數(shù)function 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)初值條件

18、改寫為：%初值條件函數(shù)function u0=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,pdefun,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 formulation, computation, and plotting of the solution of a singl

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