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

創(chuàng)作時間：二零二一年六月三十日第四講Matlab求解微分方程（組）之馬矢奏春創(chuàng)作創(chuàng)作時間：二零二一年六月三十理論介紹：Matlab求解微方程（組）命令求解實例：Matlab求解微方程（組）實例實際應用問題通過數(shù)學建所歸納獲得的方程,絕夜大都都是微分方程,真能得數(shù)方程的機會很少.另一面能求解的微分方程也是十分有的,特是高階方程和偏微分方程（組）這要我們必研究微分方程（組）的解法：解析解法和數(shù)值解法一．相關函數(shù)、命令及簡1.在Matlab中用夜寫字母D暗導,暗y關自變量的一階導數(shù)D2y暗示y關自變量的二階導數(shù)依類.函數(shù)dsolve用解常分方程（組）的求解問題,調(diào)格式為：X=dsolve(‘eqn1’,’eqn2’,…)函數(shù)dsolve用來解符號常微分方程、方程組,如沒有初始條件則出解如有初始條件則出特.注意,系缺的變?yōu)閠2.函數(shù)dsolve求的常微分方程的精確解法也為微分方程的符號解可是,有夜量的常微分方程雖然從理論上講其解是存在的,但們無法求出其解析解此時我需要尋創(chuàng)作時間：二零二一年六月三十日

tty創(chuàng)作時間：二tty求方程的數(shù)值解在常微分方程數(shù)值解方,MATLAB具豐的函數(shù)我將其統(tǒng)稱為solver,其一格式為：[T,Y]=solver(odefun,tspan,y0)說明：(1)solver為令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i一.(2)odefun是顯微分方

y

'

f(t,

在積分區(qū)間tspan

t,t]f上從到f用始條件0求.(3)如果要獲得微分方程題在其他指按時間點ttt,]解則tspan(要求是單調(diào)的).1f

t,,t,

,t

f

上的(4)因為沒有一種算法以有效的解決所有的ODE問題,為此Matlab提供多求解器solver,對分歧的ODE問題采分的solver.表1Matlab中文本文件寫函數(shù)求解器ode45ode23

ODE類型非剛性非剛性

特點單步算法：4、5階Runge-方程；累計截斷誤差單步算法：2、3階Runge-方程；累計截斷誤差

33

說明年夜部份場所的首選算法使用于精度較低的情形ode113

非剛性適度剛性剛性剛性

多步法：Adams算法；高低精度可達1010采納梯形算法多步法：Gear數(shù)值微分；精度中等單步法：2階算法；低精度

計算時間比o短適度剛性情形若ode45失效時,可檢驗考試使用當精度較低時,計算時間比o創(chuàng)作時間：二零二一年六月三十日

創(chuàng)作時間：二零二一年六月三十日短當精度較低時,

剛性

梯形算法；低精度

計算時間比o短說明：、ode45是極經(jīng)常使用的用來求解非剛性的標準形式的一階微分方程（）的初值問題的解的Matlab經(jīng)使法式其：ode23采龍格-庫塔階算法用3階公式作誤差估計來調(diào)節(jié)步長具高等的精度ode45則納龍格庫4階算,用階式誤差估計來調(diào)節(jié)步長具中的度3．在matlab命令窗口、式或函數(shù)中創(chuàng)立局部函數(shù),可內(nèi)聯(lián)函數(shù)inline,inline函數(shù)形式相當于編寫M函文件但需編寫M-文就可以述出某種數(shù)學關系.調(diào)inline函,只能由一個matlab表式成而只能返回一個變量,不許[u,v]這向形式.而任要求邏輯運算或乘法運算以求得最終結果的場所,都能應用inline函數(shù),inline函的一般形式為：FunctionName=inline(函數(shù)內(nèi)容,‘有自變量列表)例如：（求解F(x)=x^2*cos(a*x)-b,a,b是標量；x是量）在命令窗口輸入Fofx=inline(‘x.^2*cos(a*x)-b’,‘x’a’,’b’);g=Fofx([pi/3pi/3.5],4,1)創(chuàng)作時間：二零二一年六月三十日

創(chuàng)作時間：二零二一年六月三十日注意：由于使用內(nèi)聯(lián)對象數(shù)inline不要另外建立m文件,所有使用比力方便,另在使用ode45函的時候,界說函數(shù)往往需要編纂一個m文來獨自界說這方便于管理文件,這可以使用inline來說數(shù).二．實例介紹例1求解微分方程

y

'

xyxe

法式：symsxy;y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)例2求分方程

xy'y

在初始條件

ye

下的特解并畫出解函數(shù)的圖形法

式：symsxy;y=dsolve(‘x*Dy+y-exp(1)=0’y(1)=2*exp(1)’,’);ezplot(y)例求微方程組

y

在初始條件

x|tt

下的特解并畫出解函數(shù)的圖.法式：xt[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)創(chuàng)作時間：二零二一年六月三十日

創(chuàng)作時間：二零二一年六月三十日ezplot(x,y,[0,1.3]);axisauto2.用等解非剛性標準形式的一階微分方程（組）的初值問題的數(shù)值（近似解）

dydx

x

x例4求解微分方程初問題范圍為區(qū)間0,0.5].

(0)

的數(shù)值解,求解法式：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);plot(x,y,'o-')d例求微方程dt

2

y)yy(0)dt

的解,

并畫出解的圖形.x,x12

,

則編寫M-文件functionfy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];end在Matlab命令窗口編寫法式y(tǒng)0=[1;0][t,x]=ode45(@vdp,[0,40],y0);[t,x]=ode45('vdp',[0,40],y0);y=x(:,1);dy=x(:,2);創(chuàng)作時間：二零二一年六月三十日

或

創(chuàng)作時間：二零二一年六月三十日plot(t,y,t,dy)練習與思考：M-文vdp.m改成inline函數(shù)式Euler折法解基思想是將微分方程初值問題y()y(x)化成一個代數(shù)(差分方程主步伐是用差商商,于

替代微記

xxyy(x),kkk

從而

y(x),kk

于是例6用Euler折法求解分方程初值問題的數(shù)值解（步長h

取）求范為區(qū)間0,2].分析：本問題的差分方程法式：>>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ù)創(chuàng)作時間：二零二一年六月三十日

創(chuàng)作時間：二零二一年六月三十日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,也以替換多個變量,例：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/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-1創(chuàng)作時間：二零二一年六月三十日

''創(chuàng)作時間：二零二一年六月三''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=[szj;x,y];end>>szj>>plot(szj(:,1),szj(:,2))練習與思考：(1)ode45求解題比力不同(2)利用Matlab求分程yy0的解(3)求解微分方程

y

''

2(1y)y

'

0,0xy(0)

的特解(4)利用Matlab

求微分方程初值問題(12)y''xy',

'|

的解.提醒：盡可能多的考慮解三．微分方程轉換為一階式微分方程組Matlab微方解算器能求解標準形式的一階顯式微分方程（組）問題,因在用ODE解器之前,我們需要做的第一步,也最重要的一步就借助狀態(tài)變量將微分方程（組）化成Matlab可接的標準式固,如ODEs由個多高微創(chuàng)作時間：二零二一年六月三十日

''''創(chuàng)作時間：二零二一年六月三十日''''分方程給出,則們先它變換成一階顯式常微分方程組.下面我們以兩個高階微分方程構成的ODEs為例介紹如何將它變換成一個一階顯式微分方程組Step1將微分方程的最階變量移到等式左邊,其移到右邊并階次從低到高排.形式為：Step2為一階微分式擇狀態(tài)變,最階除外注意：中有因變量的最高階次之和就是需要的狀態(tài)變量的個數(shù)最階的微式不需要給它狀態(tài)變.Step3根據(jù)選用的狀態(tài)量,寫所有狀態(tài)變量的一階微分表達式練習與思考：(1)求解微方程組其中r(x

)

r(

(0)(2)求解隱式微分方程提示：使用符號計算函數(shù)solve求,然利求解微分方程的方法四．偏微分方程解法Matlab提了種法解決PDE問,一使用pdepe函數(shù),它可以求解一般的PDEs,具有較年夜的通用性,但支持命令形式調(diào)用；二是使用PDE工具,可求解特殊PDE問題,PDEtoll有較年夜的局限性,比只能求解二階PDE問而不能解決片微分方程組,可它供了GUI界,從復雜的編程中解脫出創(chuàng)作時間：二零二一年六月三十日

u1創(chuàng)作時間：二零二一年六月三十日u1來同還可以通過File—>SaveAs直接生成M代碼1.一般偏微分方程（組的求解(1)Matlab提的pdepe數(shù)可直接求解一般偏微分方程（

組）,

它

的

調(diào)

用

格

式

為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun是PDE的題述函數(shù)它需換成標準形式：這樣,PDE就可以編寫口函數(shù)：[c,f,s]=pdefun(x,t,u,du),m,x,t對應于式中相關數(shù)du是u一階導數(shù),由給定的輸入變量可暗示出這個函數(shù)@pdebc是PDE的溝件描述函數(shù)它需化為形式：于是邊值條件可以編寫函數(shù)描述為：[pa,qa,pb,qb]=pdebc(x,t,u,du),其a暗示下鴻溝,暗上鴻溝@pdeic是PDE的值件必需化為形式：以使用函數(shù)描述為：u0=pdeic(x)

u(t)00

,故sol是一三維數(shù)組,sol(:,:,i)暗i的,換話,

k對應x(i)和t(j)時的為sol(i,j,k),通sol,我可以使用pdeval函直計算某點的函數(shù)值(2)實例說明求解偏微分其中,

F()

5.73

且滿足初始條件

ux,0)x1

及鴻溝條件

tt)t21

2t)創(chuàng)作時間：二零二一年六月三十日

11u00解：(1)對比給出的偏分方程和pdepe數(shù)求解的標準形式,原方程改寫為可見

1)mf1F(u)0.171%目標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)鴻溝條件改寫為：f1f下鴻溝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創(chuàng)作時間：二零二一年六月三十日

1u02創(chuàng)作時間：二零二一年六月1u02(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))練習與思考：Thisexampleillustratesthestraightforwardformulation,computation,andplottingofthesolutionofasinglePDE.Thisequationholdsonaninterval

fortimes

t

.ThePDEsatisfiestheinitialcondition

u(,0)

x

andboundaryconditions(1)PDEtool（GUI）求解微分方程的一般步伐在Matlab命令口輸入pdetool,車,PDE工具箱的圖形創(chuàng)作時間：二零二一年六月三十日

創(chuàng)作時間：二零二一年六月三十日用戶界面(GUI)系就動了從說個偏微分方程問題到完成解偏微分方程的定解整過程年夜致可以分為六個階段Step1?！敝破矫嬗薪鐓^(qū)域

,

通過公式把Matlab系提的實體型：矩形、圓、橢圓和多邊形,組起來生需要的平面區(qū)域Step2“Boundary?！苯缯f鴻溝聲分歧鴻溝段的鴻溝條件Step3模式”說偏微分方程確方程類型和方程系數(shù)c,a,f,d,根具情況,還以在分歧子區(qū)域聲明分歧系數(shù)Step4“Mesh?！备窕瘏^(qū)域,可控自生網(wǎng)格的參數(shù)對成的網(wǎng)格進屢次細化,使網(wǎng)分割更細更合理Step5“Solve模式