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

1、2020/9/24,常微分方程（組）數(shù)值求解,吳鵬（rocwoods） MATLAB從零到進階,2020/9/24,主要內(nèi)容 數(shù)值求解常微分方程（組）函數(shù)概述 非剛性/剛性常微分方程問題求解 隱式微分方程（組）求解 微分代數(shù)方程(DAE)與延遲微分方程(DDE)求解 邊值問題求解,2020/9/24,第一節(jié)數(shù)值求解常微分方程（組）函數(shù)概述,2020/9/24,一、 概述,第9章介紹了符號求解各類型的微分方程組，但是能夠求得解析解的微分方程往往只是出現(xiàn)在大學(xué)課堂上，在實際應(yīng)用中，絕大多數(shù)微分方程（組）無法求得解析解。這就需要利用數(shù)值方法求解。MATLAB以數(shù)值計算見長，提供了一系列數(shù)值求解微分方

2、程的函數(shù)。 這些函數(shù)可以求解非剛性問題，剛性問題，隱式微分方程，微分代數(shù)方程等初值問題，也可以求解延遲微分方程，以及邊值問題等。,2020/9/24,二、初值問題求解函數(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,t0 tf,y0.),輸入?yún)?shù)說明： odefun 表示微分方程（組）的句柄。 tspan 微分

3、方程（組）的求解時間區(qū)間，有兩種格式t0,tf或者t0,t1,tf，兩者都以t0為初值點，根據(jù)tf自動選擇積分步長。前者返回實際求解過程中所有求解的時間點上的解，而后者只返回設(shè)定的時間點上的解。后者對計算效率沒有太大影響，但是求解大型問題時，可以減少內(nèi)存存儲。,2020/9/24,二、初值問題求解函數(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)體的形式返回解。,2020/9/

4、24,二、初值問題求解函數(shù),2. 函數(shù)介紹,2020/9/24,二、初值問題求解函數(shù),2020/9/24,三、 延遲問題以及邊值問題求解函數(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

5、 = bvpinit(x, yinit, params) sol = bvpsolver(odefun,bcfun,solinit) sol = bvpsolver(odefun,bcfun,solinit,options),2020/9/24,四、 求解前準備工作,微分方程的形式是多種多樣的，一般來說，很多高階微分方程可以通過變量替換轉(zhuǎn)化成一階微分方程組，即可以寫成下面的形式： （1） 稱為質(zhì)量矩陣，如果其非奇異的話,上式可以寫成： （2） 將（2）式右半部分用odefun表示出來（具體表現(xiàn)形式可以采用匿名函數(shù)、子函數(shù)、嵌套函數(shù)、單獨m文件等形式），就是ode45，ode23等常微分方程初值

6、問題求解的輸入?yún)?shù)odefun。 如果質(zhì)量矩陣奇異的話，(1)稱為微分代數(shù)方程組（differential algebraic equations, DAEs.），可以利用求解剛性微分方程的函數(shù)如ode15s,ode23s等來求解，從輸入形式上看，求解DAEs和求解普通的ODE很類似，主要區(qū)別是需要給微分方程求解器指定質(zhì)量矩陣。 隱式微分方程無法寫成(1)或者(2)的形式，其求解方法本章也有討論。,2020/9/24,第二節(jié) 非剛性/剛性常微分方程初值問題求解,2020/9/24,一、 概述,所謂剛性、非剛性問題最直觀的判別方法就是從 解在某段時間區(qū)間內(nèi)的變化來看。非剛性問題變化相對緩慢，而剛

7、性問題在某段時間內(nèi)會發(fā)生劇烈變化，即很短的時間內(nèi)，解的變化巨大。對于剛性問題不適合用ode45來求解，如果硬要用ode45來求解的話，達到指定精度所耗費的時間往往會非常長 。,2020/9/24,二、 非剛性問題舉例,問題見書中【例12.2-1】，左圖微分方程的解，右圖平面相軌跡,2020/9/24,三、 剛性問題舉例,問題見書中【例12.2-2】， 【例12.2-3】。下圖是【例12.2-2】不同求解器得到解的圖像對比。,2020/9/24,三、 剛性問題舉例,下圖是【例12.2-3】得到解的圖像，以及兩個解的和的圖像,2020/9/24,第三節(jié) 隱式微分方程（組）求解,2020/9/24,

8、一、 概述,一些 微分方程組在初始給出的時候是不容易顯示的表示成上面提到的標準形式的。這時候就需要想辦法表示成上述的形式。一般來說有三種思路，一種是利用solve函數(shù)符號求解出高階微分的顯式表達式，一種是利用fzero/fsolve函數(shù)求解狀態(tài)變量的微分值，還有一種是利用MATLAB自帶的ode15i函數(shù) 。,2020/9/24,二、利用solve函數(shù),問題見書中【例12.3-1】，下圖是求解出的結(jié)果曲線,2020/9/24,三、 利用fzero/fsolve函數(shù),問題見書中【例12.3-2】， 【例12.3-3】， 【例12.3-4】 。下圖是【例12.3-2】結(jié)果圖像。,2020/9/24

9、,三、 利用fzero/fsolve函數(shù),下圖是【例12.3-3】結(jié)果圖像。 【例12.3-4】是利用ode15i求解【例12.3-3】算例，速度明顯增快，結(jié)果一致，圖像也是下圖。,2020/9/24,第四節(jié) 微分代數(shù)方程(DAE)與延遲微分方程(DDE)求解,2020/9/24,一、 微分代數(shù)方程（DAE）舉例,DAE的求解一般有三種辦法，一種是變量替換法，一種是用ode15s函數(shù)還有一種是用12.3節(jié)中提到的ode15i函數(shù) 【例12.4-1】是利用上述三種方法求解的普通微分代數(shù)方程 。 【例12.4-2】是變量替換后用fsolve函數(shù)求解出每一計算節(jié)點的值，然后再調(diào)用ode45、ode2

10、3tb等函數(shù)求解，另一種方法就是直接利用ode15i函數(shù)求解 。,2020/9/24,一、 微分代數(shù)方程（DAE）舉例,【例12.4-1】的結(jié)果圖：,2020/9/24,一、 微分代數(shù)方程（DAE）舉例,【例12.4-2】的結(jié)果圖：,2020/9/24,二、 延遲微分方程（DDE）舉例,延遲微分方程是微分方程表達式要依賴某些狀態(tài)變量過去一些時刻的狀態(tài)，即形如： 其中， 是時間延遲項。既可以是常數(shù)也可以是關(guān)于 和 的函數(shù)。當是常數(shù)的時候可以用dde23和ddesd來求解，另一種情況可以用ddesd 求解。 【例12.4-3】是延遲為常數(shù)的求解示例。 【例12.4-4】是延遲不為常數(shù)的求解示例。,

11、2020/9/24,二、 延遲微分方程（DDE）舉例,【例12.4-3】的結(jié)果圖：,2020/9/24,二、 延遲微分方程（DDE）舉例,【例12.4-4】的結(jié)果圖：,2020/9/24,第五節(jié) 邊值問題求解,2020/9/24,一、 概述,前面討論的ode系列函數(shù)只能用來求解初值問題，但是在實際中經(jīng)常可以遇到一些邊值問題。譬如熱傳導(dǎo)問題，初值時候的熱源狀態(tài)已知，一定時間后溫度達到均勻。再比如弦振動問題，弦兩端端點的位置是固定的。像這種知道自變量在前后兩端時系統(tǒng)狀態(tài)的問題被稱為邊值問題，可以使用下面方程來描述： 定解條件：從中兩端點0和 的兩個表達式中各選一個組成定界條件 。MATLAB中提供了bvp4c和bvp5c求解邊值問題