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

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

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

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

5、格-庫塔4階算法，用5階公式作誤差估計來調節步長，具有中等的精度.3在matlab命令窗口、程序或函數中創建局部函數時，可用內聯函數inline，inline函數形式相當于編寫M函數文件，但不需編寫M-文件就可以描述出某種數學關系.調用inline函數，只能由一個matlab表達式組成，并且只能返回一個變量，不允許u,v這種向量形式.因而，任何要求邏輯運算或乘法運算以求得最終結果的場合，都不能應用inline函數，inline函數的一般形式為：FunctionName=inline(函數內容, 所有自變量列表)例如：（求解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)系統輸出為：g=-1.5483 -1.7259注意：由于使用內聯對象函數inline不需要另外建立m文件，所有使用比較方便，另外在使用ode45函數的時候，定義函數往往需要編輯一個m文件來單獨定義，這樣不便于管理文件，這里可以使用inline來定義函數.二實例介紹1.幾個可以直接用Matlab求微分方程精確解的實例例1 求解微分方程程序：syms x y; y=dsolve(Dy+2*x*y=x*exp(-x2),x)例2 求微分方程在初始條件下的特解并畫出解函數

7、的圖形.程序：syms x y; y=dsolve(x*Dy+y-exp(1)=0,y(1)=2*exp(1),x);ezplot(y)例 3 求解微分方程組在初始條件下的特解并畫出解函數的圖形.程序：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、題的數值解（近似解）例 4 求解微分方程初值問題的數值解，求解范圍為區間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)練習與思考：M-文件vdp.m改寫成inline函數程序？3.用Euler折線法求解Euler折線法求解的基本思想是將微分方程初值問題化成一個代數(差分)方程，主要步驟是用差商替代微商，于是記從而于是例 6 用Euler折線法求解微分方程初值問題的數值解（步長取0.4），求解范圍為區間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;%數值解>> for i=1:n-1 y=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，也可以替

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法的迭代公式為相應的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;%數值解>> for i=1:n-1 l1=subs(f, 'x','y',x,y);替換函數 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)練習與思考：(1)ode45求解問題并比較差異.(2)利用Matlab求微分方程的解.(3)求解微分方程的特解.(4)利用Matlab求微分方程初值問題的解.提醒：盡可能多的考慮解法三微分方程轉換為一階顯式微分方程組 Matlab微分方程解算器只能求解標準形式的一階顯式微分方程（組）問題，因此在使用ODE解算器之前，我們需要做的第一步，也是最重要的一步就是借助狀態變量將微分方程（組）化成Matlab可接受的標準形式.當然，如果ODEs由一個或多個高階微分方程給出，則我們應先將它變換成一階顯式

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

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

16、：c,f,s=pdefun(x,t,u,du)，m,x,t對應于式中相關參數，du是u的一階導數，由給定的輸入變量可表示出c,f,s這三個函數.pdebc是PDE的邊界條件描述函數，它必須化為形式：于是邊值條件可以編寫函數描述為：pa,qa,pb,qb=pdebc(x,t,u,du),其中a表示下邊界，b表示上邊界.pdeic是PDE的初值條件，必須化為形式：，故可以使用函數描述為：u0=pdeic(x)sol是一個三維數組，sol(:,:,i)表示的解，換句話說，對應x(i)和t(j)時的解為sol(i,j,k)，通過sol，我們可以使用pdeval函數直接計算某個點的函數值.(2)實例說明

17、求解偏微分其中，且滿足初始條件及邊界條件解：(1)對照給出的偏微分方程和pdepe函數求解的標準形式，原方程改寫為可見%目標PDE函數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)邊界條件改寫為：下邊界上邊界%邊界條件函數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、改寫為：%初值條件函數function u0=pdeic(x)u0=1;0;end(4)編寫主調函數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)練習與思考： 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)系統就啟動了.從定義一個偏微分方程問題到完成解偏微分方程的定解，整個過程大致可以分為六個階段Step 1 “Draw模式”繪制平面有界區域，通過公式把Matlab系統提供的實體模型：矩形、圓、橢圓和多邊形，組合起來，生成需要的平面區域.Step 2 “Boundary模式”定義邊界，聲明不同邊界段的邊界條件.Step 3 “PDE模式”定義偏微分方程，確定方程類型和方程系數c,a,f,d，根據具體情況，還可以在不同子區域聲明不同系數.Step 4 “Mesh模式”網格化區域，可以控制自動生成網格的參數，對生成的網格進行多次細化，使網格分割更細更合理.Step 5 “Solve模式”解偏微分方程