第12章常微分方程(組)數(shù)值求解方程與方程組的數(shù)值解

2023/2/3常微分方程（組）數(shù)值求解吳鵬（rocwoods）rocwoods@126.comMATLAB從零到進(jìn)階2023/2/3主要內(nèi)容數(shù)值求解常微分方程（組）函數(shù)概述非剛性/剛性常微分方程問題求解

隱式微分方程（組）求解微分代數(shù)方程(DAE)與延遲微分方程(DDE)求解邊值問題求解2023/2/3第一節(jié)數(shù)值求解常微分方程（組）函數(shù)概述2023/2/3一、概述

第9章介紹了符號求解各類型的微分方程組，但是能夠求得解析解的微分方程往往只是出現(xiàn)在大學(xué)課堂上，在實(shí)際應(yīng)用中，絕大多數(shù)微分方程（組）無法求得解析解。這就需要利用數(shù)值方法求解。MATLAB以數(shù)值計(jì)算見長，提供了一系列數(shù)值求解微分方程的函數(shù)。 這些函數(shù)可以求解非剛性問題，剛性問題，隱式微分方程，微分代數(shù)方程等初值問題，也可以求解延遲微分方程，以及邊值問題等。2023/2/3二、初值問題求解函數(shù)1.提供的函數(shù) ode23,ode45,ode113,ode15s,ode23s,ode23t,ode23tb，這些函數(shù)統(tǒng)一的調(diào)用格式如下：

[T,Y]=solver(odefun,tspan,y0) [T,Y]=solver(odefun,tspan,y0,options)sol=solver(odefun,[t0tf],y0...)

輸入?yún)?shù)說明：

odefun表示微分方程（組）的句柄。

tspan微分方程（組）的求解時(shí)間區(qū)間，有兩種格式[t0,tf]或者[t0,t1,…,tf]，兩者都以t0為初值點(diǎn)，根據(jù)tf自動(dòng)選擇積分步長。前者返回實(shí)際求解過程中所有求解的時(shí)間點(diǎn)上的解，而后者只返回設(shè)定的時(shí)間點(diǎn)上的解。后者對計(jì)算效率沒有太大影響，但是求解大型問題時(shí)，可以減少內(nèi)存存儲。

2023/2/3二、初值問題求解函數(shù)y0：微分方程（組）的初值，即所有狀態(tài)變量在t0時(shí)刻的值。

options結(jié)構(gòu)體，通過odeset設(shè)置得到的微分優(yōu)化參數(shù)。返回參數(shù)說明：T：時(shí)間點(diǎn)組成的列向量Y：微分方程（組）的解矩陣，每一行對應(yīng)相應(yīng)T的該行上時(shí)間點(diǎn)的微分方程（組）的解。sol：以結(jié)構(gòu)體的形式返回解。2023/2/3二、初值問題求解函數(shù)2.函數(shù)介紹函數(shù)問題類型精確度說明ode45非剛性中等采用算法為4-5階Runge-Kutta法，大多數(shù)情況下首選的函數(shù)ode23非剛性低基于Bogacki-Shampine2-3階Runge-Kutta公式，在精度要求不高的場合，以及對于輕度剛性方程，ode23的效率可能好于ode45。ode113非剛性低到高基于變階次Adams-Bashforth-MoutlonPECE算法。在對誤差要求嚴(yán)格的場合或者輸入?yún)?shù)odefun代表的函數(shù)本身計(jì)算量很大情況下比ode45效率高。ode113可以看成一個(gè)多步解算器，因?yàn)樗鼤?huì)利用前幾次時(shí)間節(jié)點(diǎn)上的解計(jì)算當(dāng)前時(shí)間節(jié)點(diǎn)的解。因此它不適應(yīng)于非連續(xù)系統(tǒng)。2023/2/3二、初值問題求解函數(shù)ode15s剛性低到中基于數(shù)值差分公式(后向差分公式，BDFs也叫Gear方法)，因此效率不是很高。同ode113一樣，ode15s也是一個(gè)多步計(jì)算器。當(dāng)ode45求解失敗，或者非常慢，并且懷疑問題是剛性的，或者求解微分代數(shù)問題時(shí)可以考慮用ode15sode23s剛性低基于修正的二階Rosenbrock公式。由于是單步解算器，當(dāng)精度要求不高時(shí)，它效率可能會(huì)高于ode15s。它可以解決一些ode15s求解起來效率不太高的剛性問題。

ode23t適度剛性低ode23t可以用來求解微分代數(shù)方程。ode23tb剛性低當(dāng)方程是剛性的，并且求解要求精度不高時(shí)可以使用。2023/2/3三、

延遲問題以及邊值問題求解函數(shù)1.延遲問題 MATLAB提供了dde23和ddesd函數(shù)用來求解。前者用來求解狀態(tài)變量延遲為常數(shù)的微分方程（組），后者用來求解狀態(tài)變量延遲不為常數(shù)的微分方程（組）。調(diào)用格式以及參數(shù)意義大部分類似ode系列求解函數(shù)，不同的是要輸入延遲參數(shù)以及系統(tǒng)在時(shí)間小于初值時(shí)的狀態(tài)函數(shù)。2.邊值問題 兩個(gè)求解函數(shù)函數(shù)bvp4c和bvp5c,后者求解精度要比前者好。以bvpsolver表示bvp4c或者bvp5c，那么這兩個(gè)函數(shù)有著統(tǒng)一的調(diào)用格式：solinit=bvpinit(x,yinit,params)sol=bvpsolver(odefun,bcfun,solinit)sol=bvpsolver(odefun,bcfun,solinit,options)2023/2/3四、

求解前準(zhǔn)備工作

微分方程的形式是多種多樣的，一般來說，很多高階微分方程可以通過變量替換轉(zhuǎn)化成一階微分方程組，即可以寫成下面的形式： （1） 稱為質(zhì)量矩陣，如果其非奇異的話,上式可以寫成： （2） 將（2）式右半部分用odefun表示出來（具體表現(xiàn)形式可以采用匿名函數(shù)、子函數(shù)、嵌套函數(shù)、單獨(dú)m文件等形式），就是ode45，ode23等常微分方程初值問題求解的輸入?yún)?shù)odefun。 如果質(zhì)量矩陣奇異的話，(1)稱為微分代數(shù)方程組（differentialalgebraicequations,DAEs.），可以利用求解剛性微分方程的函數(shù)如ode15s,ode23s等來求解，從輸入形式上看，求解DAEs和求解普通的ODE很類似，主要區(qū)別是需要給微分方程求解器指定質(zhì)量矩陣。 隱式微分方程無法寫成(1)或者(2)的形式，其求解方法本章也有討論。2023/2/3第二節(jié)非剛性/剛性常微分方程初值問題求解2023/2/3一、概述

所謂剛性、非剛性問題最直觀的判別方法就是從解在某段時(shí)間區(qū)間內(nèi)的變化來看。非剛性問題變化相對緩慢，而剛性問題在某段時(shí)間內(nèi)會(huì)發(fā)生劇烈變化，即很短的時(shí)間內(nèi)，解的變化巨大。對于剛性問題不適合用ode45來求解，如果硬要用ode45來求解的話，達(dá)到指定精度所耗費(fèi)的時(shí)間往往會(huì)非常長。2023/2/3二、非剛性問題舉例

問題見書中【例12.2-1】，左圖微分方程的解，右圖平面相軌跡2023/2/3三、剛性問題舉例

問題見書中【例12.2-2】，【例12.2-3】。下圖是【例12.2-2】不同求解器得到解的圖像對比。2023/2/3三、剛性問題舉例

下圖是【例12.2-3】得到解的圖像，以及兩個(gè)解的和的圖像2023/2/3第三節(jié)隱式微分方程（組）求解2023/2/3一、概述

一些微分方程組在初始給出的時(shí)候是不容易顯示的表示成上面提到的標(biāo)準(zhǔn)形式的。這時(shí)候就需要想辦法表示成上述的形式。一般來說有三種思路，一種是利用solve函數(shù)符號求解出高階微分的顯式表達(dá)式，一種是利用fzero/fsolve函數(shù)求解狀態(tài)變量的微分值，還有一種是利用MATLAB自帶的ode15i函數(shù)。2023/2/3二、利用solve函數(shù)

問題見書中【例12.3-1】，下圖是求解出的結(jié)果曲線0510152025300.511.522.533.5t

y1(t)y2(t)2023/2/3三、利用fzero/fsolve函數(shù)

問題見書中【例12.3-2】，【例12.3-3】，【例12.3-4】。下圖是【例12.3-2】結(jié)果圖像。0246810121416182000.81.8t2023/2/3三、利用fzero/fsolve函數(shù)下圖是【例12.3-3】結(jié)果圖像?！纠?2.3-4】是利用ode15i求解【例12.3-3】算例，速度明顯增快，結(jié)果一致，圖像也是下圖。2023/2/3第四節(jié)微分代數(shù)方程(DAE)與延遲微分方程(DDE)求解2023/2/3一、微分代數(shù)方程（DAE）舉例

DAE的求解一般有三種辦法，一種是變量替換法，一種是用ode15s函數(shù)還有一種是用12.3節(jié)中提到的ode15i函數(shù)

【例12.4-1】是利用上述三種方法求解的普通微分代數(shù)方程。【例12.4-2】是變量替換后用fsolve函數(shù)求解出每一計(jì)算節(jié)點(diǎn)的值，然后再調(diào)用ode45、ode23tb等函數(shù)求解，另一種方法就是直接利用ode15i函數(shù)求解。2023/2/3一、微分代數(shù)方程（DAE）舉例

【例12.4-1】的結(jié)果圖：2023/2/3一、微分代數(shù)方程（DAE）舉例

【例12.4-2】的結(jié)果圖：2023/2/3二、延遲微分方程（DDE）舉例

延遲微分方程是微分方程表達(dá)式要依賴某些狀態(tài)變量過去一些時(shí)刻的狀態(tài)，即形如： 其中，是時(shí)間延遲項(xiàng)。既可以是常數(shù)也可以是關(guān)于和的函數(shù)。當(dāng)是常數(shù)的時(shí)候可以用dde23和ddesd來求解，另一種情況可以用ddesd求解。

【例12.4-3】是延遲為常數(shù)的求解示例?！纠?2.4-4】是延遲不為常數(shù)的求解示例。2023/2/3二、延遲微分方程（DDE）舉例 【例12.4-3】的結(jié)果圖：2023/2/3二、延遲微分方程（DDE）舉例 【例12.4-4】的結(jié)果圖：2023/2/3第五節(jié)邊值問題求解2023/2/3一、概述

前面討論的ode系列函數(shù)只能用來求解初值問題，但是在實(shí)際中經(jīng)?？梢杂龅揭恍┻呏祮栴}。譬如熱傳導(dǎo)問題，初值時(shí)候的熱源狀態(tài)已知，一定時(shí)間后溫度達(dá)到均勻。再比如弦振動(dòng)問題，弦兩端端點(diǎn)的位置是固定的。像這種知道自變量在前后兩端時(shí)系統(tǒng)狀態(tài)的問題被稱為邊值問題，可以使用下面方程來描述： 定解條件：從 中兩端點(diǎn)0和的兩個(gè)表達(dá)式中各選一個(gè)組成定界條件