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

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

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

第9章介紹了符號求解各類型的微分方程組，但是能夠求得解析解的微分方程往往只是出現(xiàn)在大學(xué)課堂上，在實際應(yīng)用中，絕大多數(shù)微分方程（組）無法求得解析解。這就需要利用數(shù)值方法求解。MATLAB以數(shù)值計算見長，提供了一系列數(shù)值求解微分方程的函數(shù)。 這些函數(shù)可以求解非剛性問題，剛性問題，隱式微分方程，微分代數(shù)方程等初值問題，也可以求解延遲微分方程，以及邊值問題等。2023/11/14二、初值問題求解函數(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微分方程（組）的求解時間區(qū)間，有兩種格式[t0,tf]或者[t0,t1,…,tf]，兩者都以t0為初值點，根據(jù)tf自動選擇積分步長。前者返回實際求解過程中所有求解的時間點上的解，而后者只返回設(shè)定的時間點上的解。后者對計算效率沒有太大影響，但是求解大型問題時，可以減少內(nèi)存存儲。

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

options結(jié)構(gòu)體，通過odeset設(shè)置得到的微分優(yōu)化參數(shù)。返回參數(shù)說明：T：時間點組成的列向量Y：微分方程（組）的解矩陣，每一行對應(yīng)相應(yīng)T的該行上時間點的微分方程（組）的解。sol：以結(jié)構(gòu)體的形式返回解。2023/11/14二、初值問題求解函數(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ù)本身計算量很大情況下比ode45效率高。ode113可以看成一個多步解算器，因為它會利用前幾次時間節(jié)點上的解計算當(dāng)前時間節(jié)點的解。因此它不適應(yīng)于非連續(xù)系統(tǒng)。2023/11/14二、初值問題求解函數(shù)ode15s剛性低到中基于數(shù)值差分公式(后向差分公式，BDFs也叫Gear方法)，因此效率不是很高。同ode113一樣，ode15s也是一個多步計算器。當(dāng)ode45求解失敗，或者非常慢，并且懷疑問題是剛性的，或者求解微分代數(shù)問題時可以考慮用ode15sode23s剛性低基于修正的二階Rosenbrock公式。由于是單步解算器，當(dāng)精度要求不高時，它效率可能會高于ode15s。它可以解決一些ode15s求解起來效率不太高的剛性問題。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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