第四講Matlab求解微分方程〔組〕理論介紹：Matlab求解微分方程〔組〕命令求解實(shí)例：Matlab求解微分方程〔組〕實(shí)例實(shí)際應(yīng)用問(wèn)題通過(guò)數(shù)學(xué)建模所歸納得到的方程，絕大多數(shù)都是微分方程，真正能得到代數(shù)方程的時(shí)機(jī)很少.另一方面，能夠求解的微分方程也是十分有限的，特別是高階方程和偏微分方程〔組〕.這就要求我們必須研究微分方程〔組〕的解法：解析解法和數(shù)值解法.一．相關(guān)函數(shù)、命令及簡(jiǎn)介1.在Matlab中，用大寫(xiě)字母D表示導(dǎo)數(shù)，Dy表示y關(guān)于自變量的一階導(dǎo)數(shù)，D2y表示y關(guān)于自變量的二階導(dǎo)數(shù)，依此類(lèi)推.函數(shù)dsolve用來(lái)解決常微分方程〔組〕的求解問(wèn)題，調(diào)用格式為：X=dsolve(‘eqn1’,’eqn2’,…)函數(shù)dsolve用來(lái)解符號(hào)常微分方程、方程組，如果沒(méi)有初始條件，那么求出通解，如果有初始條件，那么求出特解.注意，系統(tǒng)缺省的自變量為t2.函數(shù)dsolve求解的是常微分方程的精確解法，也稱為常微分方程的符號(hào)解.但是，有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無(wú)法求出其解析解，此時(shí)，我們需要尋求方程的數(shù)值解，在求常微分方程數(shù)值解方面，MATLAB具有豐富的函數(shù)，我們將其統(tǒng)稱為solver，其一般格式為：[T,Y]=solver(odefun,tspan,y0)說(shuō)明：(1)solver為命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i之一.(2)odefun是顯示微分方程在積分區(qū)間tspan上從到用初始條件求解.(3)如果要獲得微分方程問(wèn)題在其他指定時(shí)間點(diǎn)上的解，那么令tspan(要求是單調(diào)的).(4)因?yàn)闆](méi)有一種算法可以有效的解決所有的ODE問(wèn)題，為此，Matlab提供了多種求解器solver，對(duì)于不同的ODE問(wèn)題，采用不同的solver.表1Matlab中文本文件讀寫(xiě)函數(shù)求解器ODE類(lèi)型特點(diǎn)說(shuō)明ode45非剛性單步算法：4、5階Runge-Kutta方程；累計(jì)截?cái)嗾`差大局部場(chǎng)合的首選算法ode23非剛性單步算法：2、3階Runge-Kutta方程；累計(jì)截?cái)嗾`差使用于精度較低的情形ode113非剛性多步法：Adams算法；上下精度可達(dá)計(jì)算時(shí)間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法：Gear’s反向數(shù)值微分；精度中等假設(shè)ode45失效時(shí)，可嘗試使用ode23s剛性單步法：2階Rosebrock算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短說(shuō)明：ode23、ode45是極其常用的用來(lái)求解非剛性的標(biāo)準(zhǔn)形式的一階微分方程〔組〕的初值問(wèn)題的解的Matlab常用程序，其中：ode23采用龍格-庫(kù)塔2階算法，用3階公式作誤差估計(jì)來(lái)調(diào)節(jié)步長(zhǎng)，具有低等的精度.ode45那么采用龍格-庫(kù)塔4階算法，用5階公式作誤差估計(jì)來(lái)調(diào)節(jié)步長(zhǎng)，具有中等的精度.3．在matlab命令窗口、程序或函數(shù)中創(chuàng)立局部函數(shù)時(shí)，可用內(nèi)聯(lián)函數(shù)inline，inline函數(shù)形式相當(dāng)于編寫(xiě)M函數(shù)文件，但不需編寫(xiě)M-文件就可以描述出某種數(shù)學(xué)關(guān)系.調(diào)用inline函數(shù)，只能由一個(gè)matlab表達(dá)式組成，并且只能返回一個(gè)變量，不允許[u,v]這種向量形式.因而，任何要求邏輯運(yùn)算或乘法運(yùn)算以求得最終結(jié)果的場(chǎng)合，都不能應(yīng)用inline函數(shù)，inline函數(shù)的一般形式為：FunctionName=inline(‘函數(shù)內(nèi)容’,‘所有自變量列表’)例如：〔求解F(x)=x^2*cos(a*x)-b,a,b是標(biāo)量；x是向量〕在命令窗口輸入:Fofx=inline(‘x.^2*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)對(duì)象函數(shù)inline不需要另外建立m文件，所有使用比擬方便，另外在使用ode45函數(shù)的時(shí)候，定義函數(shù)往往需要編輯一個(gè)m文件來(lái)單獨(dú)定義，這樣不便于管理文件，這里可以使用inline來(lái)定義函數(shù).二．實(shí)例介紹1.幾個(gè)可以直接用Matlab求微分方程精確解的實(shí)例例1求解微分方程程序：symsxy;y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)例2求微分方程在初始條件下的特解并畫(huà)出解函數(shù)的圖形.程序：symsxy;y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x’);ezplot(y)例3求解微分方程組在初始條件下的特解并畫(huà)出解函數(shù)的圖形.程序：symsxyt[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]);axisauto2.用ode23、ode45等求解非剛性標(biāo)準(zhǔn)形式的一階微分方程〔組〕的初值問(wèn)題的數(shù)值解〔近似解〕例4求解微分方程初值問(wèn)題的數(shù)值解，求解范圍為區(qū)間[0,0.5].程序：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);plot(x,y,'o-')例5求解微分方程的解，并畫(huà)出解的圖形.分析：這是一個(gè)二階非線性方程，我們可以通過(guò)變換，將二階方程化為一階方程組求解.令，那么編寫(xiě)M-文件vdp.mfunctionfy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];end在Matlab命令窗口編寫(xiě)程序y0=[1;0][t,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改寫(xiě)成inline函數(shù)程序？3.用Euler折線法求解Euler折線法求解的根本思想是將微分方程初值問(wèn)題化成一個(gè)代數(shù)(差分)方程，主要步驟是用差商替代微商，于是記從而于是例6用Euler折線法求解微分方程初值問(wèn)題的數(shù)值解〔步長(zhǎng)取0.4〕，求解范圍為區(qū)間[0,2].分析：本問(wèn)題的差分方程為程序：>>clear>>f=sym('y+2*x/y^2');>>a=0;>>b=2;>>h=0.4;>>n=(b-a)/h+1;>>x=0;>>y=1;>>szj=[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))說(shuō)明：替換函數(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)特別說(shuō)明：本問(wèn)題可進(jìn)一步利用四階Runge-Kutta法求解，Euler折線法實(shí)際上就是一階Runge-Kutta法，Runge-Kutta法的迭代公式為相應(yīng)的Matlab程序?yàn)椋?gt;>clear>>f=sym('y+2*x/y^2');>>a=0;>>b=2;>>h=0.4;>>n=(b-a)/h+1;>>x=0;>>y=1;>>szj=[x,y];%數(shù)值解>>fori=1:n-1l1=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=[szj;x,y];end>>szj>>plot(szj(:,1),szj(:,2))練習(xí)與思考：(1)ode45求解問(wèn)題并比擬差異.(2)利用Matlab求微分方程的解.(3)求解微分方程的特解.(4)利用Matlab求微分方程初值問(wèn)題的解.提醒：盡可能多的考慮解法三．微分方程轉(zhuǎn)換為一階顯式微分方程組Matlab微分方程解算器只能求解標(biāo)準(zhǔn)形式的一階顯式微分方程〔組〕問(wèn)題，因此在使用ODE解算器之前，我們需要做的第一步，也是最重要的一步就是借助狀態(tài)變量將微分方程〔組〕化成Matlab可接受的標(biāo)準(zhǔn)形式.當(dāng)然，如果ODEs由一個(gè)或多個(gè)高階微分方程給出，那么我們應(yīng)先將它變換成一階顯式常微分方程組.下面我們以兩個(gè)高階微分方程組構(gòu)成的ODEs為例介紹如何將它變換成一個(gè)一階顯式微分方程組.Step1將微分方程的最高階變量移到等式左邊，其它移到右邊，并按階次從低到高排列.形式為：Step2為每一階微分式選擇狀態(tài)變量，最高階除外注意：ODEs中所有是因變量的最高階次之和就是需要的狀態(tài)變量的個(gè)數(shù)，最高階的微分式不需要給它狀態(tài)變量.Step3根據(jù)選用的狀態(tài)變量，寫(xiě)出所有狀態(tài)變量的一階微分表達(dá)式練習(xí)與思考：(1)求解微分方程組其中(2)求解隱式微分方程組提示：使用符號(hào)計(jì)算函數(shù)solve求，然后利用求解微分方程的方法四．偏微分方程解法Matlab提供了兩種方法解決PDE問(wèn)題，一是使用pdepe函數(shù)，它可以求解一般的PDEs,具有較大的通用性，但只支持命令形式調(diào)用；二是使用PDE工具箱，可以求解特殊PDE問(wèn)題，PDEtoll有較大的局限性，比方只能求解二階PDE問(wèn)題，并且不能解決片微分方程組，但是它提供了GUI界面，從復(fù)雜的編程中解脫出來(lái)，同時(shí)還可以通過(guò)File—>SaveAs直接生成M代碼.1.一般偏微分方程〔組〕的求解(1)Matlab提供的pdepe函數(shù)，可以直接求解一般偏微分方程〔組〕，它的調(diào)用格式為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun是PDE的問(wèn)題描述函數(shù)，它必須換成標(biāo)準(zhǔn)形式：這樣，PDE就可以編寫(xiě)入口函數(shù)：[c,f,s]=pdefun(x,t,u,du)，m,x,t對(duì)應(yīng)于式中相關(guān)參數(shù)，du是u的一階導(dǎo)數(shù)，由給定的輸入變量可表示出c,f,s這三個(gè)函數(shù).@pdebc是PDE的邊界條件描述函數(shù)，它必須化為形式：于是邊值條件可以編寫(xiě)函數(shù)描述為：[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a表示下邊界，b表示上邊界.@pdeic是PDE的初值條件，必須化為形式：，故可以使用函數(shù)描述為：u0=pdeic(x)sol是一個(gè)三維數(shù)組，sol(:,:,i)表示的解，換句話說(shuō)，對(duì)應(yīng)x(i)和t(j)時(shí)的解為sol(i,j,k)，通過(guò)sol，我們可以使用pdeval函數(shù)直接計(jì)算某個(gè)點(diǎn)的函數(shù)值.(2)實(shí)例說(shuō)明求解偏微分其中，且滿足初始條件及邊界條件解：(1)對(duì)照給出的偏微分方程和pdepe函數(shù)求解的標(biāo)準(zhǔn)形式，原方程改寫(xiě)為可見(jiàn)%目標(biāo)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)邊界條件改寫(xiě)為：下邊界上邊界%邊界條件函數(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)初值條件改寫(xiě)為：%初值條件函數(shù)functionu0=pdeic(x)u0=[1;0];end(4)編寫(xiě)主調(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í)與思考：Thisexampleillustratesthestraightforwardformulation,computation,andplottingofthesolutionofasinglePDE.Thisequationholdsonanintervalfortimes.ThePDEsatisfiestheinitialconditionandboundaryconditions2.PDEtool求解偏微分方程(1)PDEtool〔GUI〕求解偏微分方程的一般步驟在Matlab命令窗口輸入pdetool，回車(chē)，PDE工具箱的圖形用戶界面(GUI)系統(tǒng)就啟動(dòng)了.從定義一個(gè)偏微分方程問(wèn)題到完成解偏微分方程的定解，整個(gè)過(guò)程大致可以分為六個(gè)階段Step1“Draw模式〞繪制平面有界區(qū)域，通過(guò)公式把Matlab系統(tǒng)提供的實(shí)體模型：矩形、圓、橢圓和多邊形，組合起來(lái)，生成需要的平面區(qū)域.Step2“Boundary模式〞定義邊界，聲明不同邊界段的邊界條件.Step3“PDE模式〞定義偏微分方程，確定方程類(lèi)型和方程系數(shù)c,a,f,d，根據(jù)具體情況，還可以在不同子區(qū)域聲明不同系數(shù).Step4“Mesh模式〞網(wǎng)格化區(qū)域，可以控制自動(dòng)生成網(wǎng)格的參數(shù)，對(duì)生成的網(wǎng)格進(jìn)行屢次細(xì)化，使網(wǎng)格分割更細(xì)更合理.Step5“Solve模式〞解偏微分方程，對(duì)于橢圓型方程可以激