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

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

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

3、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i 之一.odefun 是顯示微分方程y=f（t,y）在積分區(qū)間tspan =b,tf上從b到tf用初始條件y0求 解.（3）如果要獲得微分方程問題在其他指定時(shí)間點(diǎn)t0,t1,t2,ni,tf上的解，則令 tspan =心142,11匕（要求是單調(diào)的）.（4）因?yàn)闆]有一種算法可以有效的解決所有的ODE問題，為此，Matlab提供了多種求解器solver ,對于不同的ODE'可題，采用不同的solver.表1Matlab中文本文件讀寫函數(shù)求解器OD#型特點(diǎn)說明ode45非剛性單步算法：4、5階R

4、unge-Kutta大部分場合的首選算法方程；累小戢斷誤差(Ax)3ode23非剛性單步算法：2、3階Runge-Kutta方程；累小戢斷誤差(Ax)3使用于精度較低的情形ode113非剛性多步法：Adams#法；高低精度可達(dá)10與10-6計(jì)算時(shí)間比ode45短ode23t適度ffl性采用梯形算法適度剛性情形ode15s剛性多步法：Gear' s反向數(shù)值微分；精度中等若ode45失效時(shí)，可 嘗試使用ode23s剛性單步法：2階Rosebrock算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短說明：ode23、od

5、e45是極其常用的用來求解非剛性的標(biāo)準(zhǔn)形式的一階微分方程(組)的初值 問題的解的Matlab常用程序，其中：ode23采用龍格-庫塔2階算法，用3階公式作誤差估計(jì)來調(diào)節(jié)步長，具有低等的精度.ode45則采用龍格-庫塔4階算法，用5階公式作誤差估計(jì)來調(diào)節(jié)步長，具有中等的精度.3 .在matlab命令窗口、程序或函數(shù)中創(chuàng)建局部函數(shù)時(shí)，可用內(nèi)聯(lián)函數(shù)inline , inline 函數(shù)形式相當(dāng)于編寫M函數(shù)文件，但不需編寫M-文件就可以描述出某種數(shù)學(xué)關(guān)系.調(diào)用inline函數(shù)， 只能由一個(gè)matlab表達(dá)式組成，并且只能返回一個(gè)變量，不允許 u,v這種向量形式.因而，任何 要求邏輯運(yùn)算或乘法運(yùn)算以求得最

6、終結(jié)果的場合，都不能應(yīng)用inline函數(shù)，inline函數(shù)的一般形式為：J' I I"""-FunctionName=inline('函數(shù)內(nèi)容，所有自變量列表) ,1 例如：(求解F(x)=xA2*cos(a*x)-b,a,b是標(biāo)量；x是向量)在命令窗口輸入：Fofx=inline( 'x.A2*cos(a*x) - b' , 'x' , ' a' , ' b');g=Fofx(pi/3pi/3.5,4,1)系統(tǒng)輸出為：g=-1.5483-1.7259注意：由于使用內(nèi)聯(lián)對象函數(shù)inli

7、ne 不需要另外建立m文件，所有使用比較方便，另外在使 用ode45函數(shù)的時(shí)候，定義函數(shù)往往需要編輯一個(gè) m文件來單獨(dú)定義，這樣不便于管理文件，這里 可以使用inline來定義函數(shù).實(shí)例介紹4 .幾個(gè)可以直接用Matlab求微分方程精確解的實(shí)例2例1求解微分方程y +2xy=xe程序：symsxy;y=dsolve(' Dy+2*x*y=x*exp(-xA2) ' , ' x')例2求微分方程xy'+y-ex=0在初始條件y(1) = 2e下的特解并畫出解函數(shù)的圖形.程序：symsxy;y=dsolve('x*Dy+y-exp(1)=0 '

8、; , ' y(1)=2*exp(1) ' , ' x' );ezplot(y)dxt5x y =e例3求解微分方程組1dt在初始條件x|y = 1,y|y = 0下的特解并畫出解函數(shù)的圖dy-x-3y =0一 一dt y形.程序：symsxytx,y=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t'):r- ：1' / j .1 :simple(x);_simple(y)ezplot(x,y,0,1.3)

9、;axisauto5 .用ode23、ode45等求解非剛性標(biāo)準(zhǔn)形式的一階微分方程(組)的初值問題的數(shù)值解(近似解),» ,、十、 ,+、 口上曳=-2v+2x2+2x -八,、例4求解微分方程初值問題dx 2y 2x2x的數(shù)值解，求解范圍為區(qū)間0,0.5.y(0) = 1程序：fun=inline('-2*y+2*xA2+2*x,x,y,);x,y=ode23(fun,0,0.5,1);plot(x,y,o-1)d2d例5求解微分萬程d-N(1-y2)dy + y =0,y(0) =1,y(0) =0的解，并回出解的圖形. dtdt分析：這是一個(gè)二階非線性方程，我們可以通過

10、變換，將二階方程化為一階方程組求解.令x1 = y, x2 =dy, N =7 ,則dt編寫M-文件vdp.mfunctionfy=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);plot(t,y,t,dy)練習(xí)與思考：M-文件vdp.m改寫成inline 函數(shù)程序？6 .用Euler折線法求解Euler折線法求解的基本思想是將微分方程初值問題化成一個(gè)代數(shù)(差分)方

11、程，主要步驟是用差商y(x + h) y(x)替代微商曳,于是 hdx記 xk + =xk +h,yk =y(xk),從而 yp = y(xk +h),于是例6用Euler折線法求解微分方程初值問題的數(shù)值解(步長h取0.4),求解范圍為區(qū)間0,2.分析：本問題的差分方程為f I,4程序：>>clear>>f=sym('y+2*x/yA2');>>a=0;>>b=2;1I 1>>h=0.4;''.I>>n=(b-a)/h+1;>>x=0;>>y=1;>>sz

12、j=x,y;% 數(shù)值解>>fori=1:n-1y=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(bǔ)用4替換掉，返回4+b,也可以替換多個(gè)變量，例如：subs(cos(a)+sin(b),a,b,sym('alpha'),2)分別用字符 alpha 替換 a和 2 替換 b,返回 cos(alpha)+sin(2)特別說明：本

13、問題可進(jìn)一步利用四階Runge-Kutta法求解，Euler折線法實(shí)際上就是一階Runge-Kutta法，Runge-Kutta法的迭代公式為相應(yīng)的Matlab程序?yàn)椋?gt;>clear>>f=sym('y+2*x/yA2');>>a=0;>>b=2;>>h=0.4;>>n=(b-a)/h+1;>>x=0;>>y=i;>>szj=x,y;% 數(shù)值解.1X '/-.'J;>>fori=1:n-1一l1=subs(f,'x','

14、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;一、 JHx=x+h;szj=szj;x,y;/ 廣 f rend>>szj>>plot(szj(:,1),szj(:,2)練習(xí)與思考：(1)ode45求解問題并比較差異.(2)利用Matlab求微分

15、方程y-2y+y'' =0的解.(3)求解微分方程 y''2(1 y2)y'+y = 0,0 < x < 30, y(0) =1,y,(0) = 0 的特解.(4)利用Matlab求微分方程初值問題(1+x )y =2xy , y心=1, y |x田=3的解.提醒：盡可能多的考慮解法三.微分方程轉(zhuǎn)換為一階顯式微分方程組ODEW 算Matlab微分方程解算器只能求解標(biāo)準(zhǔn)形式的一階顯式微分方程（組）問題，因此在使用器之前，我們需要做的第一步，也是最重要的一步就是借助狀態(tài)變量將微分方程(組)化成Matlab可接受的標(biāo)準(zhǔn)形式.當(dāng)然，如果ODE前一個(gè)或

16、多個(gè)高階微分方程給出，則我們應(yīng)先將它變換成一階 顯式常微分方程組.下面我們以兩個(gè)高階微分方程組構(gòu)成的ODEs為例介紹如何將它變換成一個(gè)一階顯式微分方程組.Stepl將微分方程的最高階變量移到等式左邊，其它移到右邊，并按階次從低到高排列.形式為：Step2為每一階微分式選擇狀態(tài)變量，最高階除外注意：ODEs中所有是因變量的最高階次之和就是需要的狀態(tài)變量的個(gè)數(shù)，最高階的微分式不 需要給它狀態(tài)變量.Step3根據(jù)選用的狀態(tài)變量，寫出所有狀態(tài)變量的一階微分表達(dá)式練習(xí)與思考：(1)求解微分方程組其中 r2 =7(x-k*)2 +y2, r1 = J(x + N)2 +丫2, N* =1 N, N =1/

17、82.45, x(0) =1.2, (2)求解隱式微分方程組提示：使用符號計(jì)算函數(shù)solve求x,y",然后利用求解微分方程的方法四.偏微分方程解法Matlab提供了兩種方法解決PD日可題，一是使用pdepe函數(shù)，它可以求解一般的PDEs具有較 大的通用性，但只支持命令形式調(diào)用；二是使用 PDEX具箱，可以求解特殊PDE問題，PDEtoll有 較大的局限性，比如只能求解二階PD日可題，并且不能解決片微分方程組，但是它提供了 GUI界面， 從復(fù)雜的編程中解脫出來，同時(shí)還可以通過 File >SaveAs直接生成M代碼.1.一般偏微分方程(組)的求解(1)Matlab提供的pdep

18、e函數(shù)，可以直接求解一般偏微分方程(組),它的調(diào)用格式為：sol=pdepe(m,pdefun,pdeic,pdebc,x,t)pdefun是PDE的問題描述函數(shù)，它必須換成標(biāo)準(zhǔn)形式：這樣，PDEB可以編寫入口函數(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這三個(gè)函數(shù).pdebc! PDE的邊界條件描述函數(shù)，它必須化為形式：于是邊值條件可以編寫函數(shù)描述為：pa,qa,pb,qb=pdebc(x,t,u,du), 其中a表示下邊界，b表示 上邊界.pdeic是PDE的初值條件，必須化為形式：u(x,b)

19、 = u0,故可以使用函數(shù)描述為：u0=pdeic(x)sol是一個(gè)三維數(shù)組，sol(:,:,i) 表示ui的解，換句話說，耳對應(yīng)x(i)和t(j)時(shí)的解為sol(i,j,k),通過sol ,我們可以使用pdeval函數(shù)直接計(jì)算某個(gè)點(diǎn)的函數(shù)值(2)實(shí)例說明求解偏微分 其中，F(xiàn)(x)=e5 .x -3- x且滿足初始條件ui(x,0) =1,u2(x,0) =0及邊界條件5(0,t) =0, U2(0,t) =0,Ul(1,t) =1,以(1,t) =0 .xfx0.024 世 1 ex0.17 也 歙-F (u1 -u2 )” 上(5 叫)_解：(1)對照給出的偏微分方程和pdepe函數(shù)求解的

20、標(biāo)準(zhǔn)形式，原方程改寫為可見 m = 0,c, f一1%g標(biāo)PDEg數(shù)functionc,f,s=pdefun(x,t,u,du) r-_ ：1' / Jc=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)邊界條件改寫為：i I下邊界1+p】.* f =° i上邊界卜-q+pi* f=，?2.| 001 0 一 |o 一1。%邊界條件函數(shù)functionpa,qa,pb,qb=pdebc(xa,ua,xb,ub,t)pa=0;ua(

21、2);Jqa=1；0;pb=ub(1)-1;0;qb=0；1；end(3)初值條件改寫為:%初值條件函數(shù)functionu0=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)I , - - .：； ；，X 練習(xí)與思考：Thisexampleillustratesthestraightforwardformulation,

22、computation,andplottingofthesolut ionofasinglePDE.一Thisequationholdsonaninterval 0< x <1 fortimes t_0 .ThePDEsatisfiestheinitialconditionu(x,0) =sin 二 xandboundaryconditions2.PDEtool求解偏微分方程(1)PDEtool (GUI)求解偏微分方程的一般步驟在Matlab命令窗口輸入pdetool ,回車，PDET具箱的圖形用戶界面(GUI)系統(tǒng)就啟動(dòng)了 .從定 義一個(gè)偏微分方程問題到完成解偏微分方程的定解，整個(gè)過程大致可以分為六個(gè)階段Step1 "Draw模式”繪制平面有界區(qū)域 C ,通過公式把Matlab系統(tǒng)提供的實(shí)體模型：矩形、 圓、橢圓和多邊形，組合起來，生成需要的平面區(qū)域 .Step2 "Boundary模式”定義邊界，聲明不同邊界段的邊界條件. I 1 'Step3 "PDE真式”定義偏微分方程，確定方程類型和方程系數(shù)c,a,f,d ,根據(jù)具體情況，還可以在不同子區(qū)域聲明不同系數(shù).Step4 "Mesh模式