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

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

tty創(chuàng)作時(shí)間：二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)如果要獲得微分方程題在其他指按時(shí)間點(diǎn)ttt,]解則tspan(要求是單調(diào)的).1f

t,,t,

,t

f

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

ODE類型非剛性非剛性

特點(diǎn)單步算法：4、5階Runge-方程；累計(jì)截?cái)嗾`差單步算法：2、3階Runge-方程；累計(jì)截?cái)嗾`差

33

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

非剛性適度剛性剛性剛性

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

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

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

剛性

梯形算法；低精度

計(jì)算時(shí)間比o短說明：、ode45是極經(jīng)常使用的用來求解非剛性的標(biāo)準(zhǔn)形式的一階微分方程（）的初值問題的解的Matlab經(jīng)使法式其：ode23采龍格-庫(kù)塔階算法用3階公式作誤差估計(jì)來調(diào)節(jié)步長(zhǎng)具高等的精度ode45則納龍格庫(kù)4階算,用階式誤差估計(jì)來調(diào)節(jié)步長(zhǎng)具中的度3．在matlab命令窗口、式或函數(shù)中創(chuàng)立局部函數(shù),可內(nèi)聯(lián)函數(shù)inline,inline函數(shù)形式相當(dāng)于編寫M函文件但需編寫M-文就可以述出某種數(shù)學(xué)關(guān)系.調(diào)inline函,只能由一個(gè)matlab表式成而只能返回一個(gè)變量,不許[u,v]這向形式.而任要求邏輯運(yùn)算或乘法運(yùn)算以求得最終結(jié)果的場(chǎng)所,都能應(yīng)用inline函數(shù),inline函的一般形式為：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)創(chuàng)作時(shí)間：二零二一年六月三十日

創(chuàng)作時(shí)間：二零二一年六月三十日注意：由于使用內(nèi)聯(lián)對(duì)象數(shù)inline不要另外建立m文件,所有使用比力方便,另在使用ode45函的時(shí)候,界說函數(shù)往往需要編纂一個(gè)m文來獨(dú)自界說這方便于管理文件,這可以使用inline來說數(shù).二．實(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)作時(shí)間：二零二一年六月三十日

創(chuàng)作時(shí)間：二零二一年六月三十日ezplot(x,y,[0,1.3]);axisauto2.用等解非剛性標(biāo)準(zhǔn)形式的一階微分方程（組）的初值問題的數(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)作時(shí)間：二零二一年六月三十日

或

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

替代微記

xxyy(x),kkk

從而

y(x),kk

于是例6用Euler折法求解分方程初值問題的數(shù)值解（步長(zhǎng)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)作時(shí)間：二零二一年六月三十日

創(chuàng)作時(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)特別說明：本問題可進(jìn)一利用四階Runge-Kutta法解,Euler折線法實(shí)際上就一階Runge-Kutta法Runge-Kutta法的迭代公式為相應(yīng)的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)作時(shí)間：二零二一年六月三十日

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

y

''

2(1y)y

'

0,0xy(0)

的特解(4)利用Matlab

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

'|

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

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

)

r(

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

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

組）,

它

的

調(diào)

用

格

式

為：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun是PDE的題述函數(shù)它需換成標(biāo)準(zhǔn)形式：這樣,PDE就可以編寫口函數(shù)：[c,f,s]=pdefun(x,t,u,du),m,x,t對(duì)應(yīng)于式中相關(guān)數(shù)du是u一階導(dǎo)數(shù),由給定的輸入變量可暗示出這個(gè)函數(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對(duì)應(yīng)x(i)和t(j)時(shí)的為sol(i,j,k),通sol,我可以使用pdeval函直計(jì)算某點(diǎn)的函數(shù)值(2)實(shí)例說明求解偏微分其中,

F()

5.73

且滿足初始條件

ux,0)x1

及鴻溝條件

tt)t21

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

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

1)mf1F(u)0.171%目標(biāo)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)作時(shí)間：二零二一年六月三十日

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

fortimes

t

.ThePDEsatisfiestheinitialcondition

u(,0)

x

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

創(chuàng)作時(shí)間：二零二一年六月三十日用戶界面(GUI)系就動(dòng)了從說個(gè)偏微分方程問題到完成解偏微分方程的定解整過程年夜致可以分為六個(gè)階段Step1模”制平面有界區(qū)域

,

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