1、9/3/2020星期六, 2008-9- 6, 13:10:16,第7章微分方程問(wèn)題的計(jì)算機(jī)求解,高等應(yīng)用數(shù)學(xué)問(wèn)題的MATLAB求解,清華大學(xué)出版社2008,CAI課件開(kāi)發(fā)：薛定宇、劉瑩瑩、董雯彬,9/3/2020星期六, 2008-9- 6, 13:10:16,第7章 微分方程問(wèn)題的計(jì)算機(jī)求解,常系數(shù)線性微分方程的解析解方法 微分方程問(wèn)題的數(shù)值解法 微分方程轉(zhuǎn)換 特殊微分方程的數(shù)值解 邊值問(wèn)題的計(jì)算機(jī)求解 偏微分方程求解入門(mén) 微分方程的框圖求解,9/3/2020星期六, 2008-9- 6, 13:10:16,7.1 常系數(shù)線性微分方程的解析解方法,線性常系數(shù)微分方程解析解的數(shù)學(xué)描述 微分方

2、程的解析解方法 Laplace變換在線性微分方程求解中的應(yīng)用 線性狀態(tài)空間方程的解析解 特殊非線性微分方程的解析解,9/3/2020星期六, 2008-9- 6, 13:10:16,7.1.1 線性常系數(shù)微分方程解析解的數(shù)學(xué)描述,常系數(shù)線性微分方程的一般描述方法為 其中， 均為常數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,對(duì)零初值問(wèn)題： 對(duì)應(yīng)得出Laplace變換為 設(shè)代數(shù)方程的特征根si均相異，則解析解為 其中Ci為待定系數(shù)，并且g(t)是滿(mǎn)足u(t)輸入的一個(gè)特解,9/3/2020星期六, 2008-9- 6, 13:10:16,簡(jiǎn)單形式，自變量設(shè)為t 指明自變量

3、x,7.1.2 微分方程的解析解方法,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.1,假設(shè)輸入信號(hào)為 試求出下面微分方程的通解,9/3/2020星期六, 2008-9- 6, 13:10:16,假設(shè)有如下初值條件,結(jié)果為：,9/3/2020星期六, 2008-9- 6, 13:10:16,9/3/2020星期六, 2008-9- 6, 13:10:16,假設(shè) 則可以獲得方程的解析解,9/3/2020星期六, 2008-9- 6, 13:10:16,最終的近似結(jié)果為：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.2,給定輸入信號(hào) ，求解

4、 其中，,9/3/2020星期六, 2008-9- 6, 13:10:16,結(jié)果：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.3,試求解線性微分方程組的解析解 直接求解,9/3/2020星期六, 2008-9- 6, 13:10:16,7.1.3 Laplace變換在線性微分方程求解中的應(yīng)用,傳遞函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,那么，輸出信號(hào)可以表示成如下的s-域函數(shù) 其中， 如果U(s) 是一個(gè)有理函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.4,給定輸入信號(hào)u(t)=e-5tcos(2

5、t+1)+5，假定所有的初始條件都為0，試求出下述微分方程的解析解,9/3/2020星期六, 2008-9- 6, 13:10:16,用Laplace變換法,9/3/2020星期六, 2008-9- 6, 13:10:16,在時(shí)域中輸出信號(hào)的解析解：,9/3/2020星期六, 2008-9- 6, 13:10:16,使用dsolve()函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,繪制解曲線： 解析解相同,9/3/2020星期六, 2008-9- 6, 13:10:16,7.1.4 線性狀態(tài)空間方程的解析解,假設(shè)線性狀態(tài)空間模型的一般表示為 其中，A,B,C,D是常數(shù)

6、矩陣，且已知狀態(tài)向量初值 ，該方程的解析解是：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.5,給定輸入信號(hào)為u(t) = 2 + 2e-3t sin 2t，求出下面矩陣描述的狀態(tài)空間方程的解析解,9/3/2020星期六, 2008-9- 6, 13:10:16,直接積分法求出結(jié)果,9/3/2020星期六, 2008-9- 6, 13:10:16,7.1.5 特殊非線性微分方程的解析解,只有少數(shù)非線性微分方程可以通過(guò)dsolve()函數(shù)得出解析解 通過(guò)例子演示非線性方程的解析解求解 同時(shí)還將演示不能求解的例子,9/3/2020星期六, 2008-9- 6, 13:

7、10:16,例 7.6,求解非線性微分方程： 改變?cè)⒎址匠痰男问?9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.7,試求出著名的Van der Pol方程的解析解 嘗試如下的MATLAB命令,9/3/2020星期六, 2008-9- 6, 13:10:16,7.2 微分方程問(wèn)題的數(shù)值解法,微分方程問(wèn)題算法概述 四階定步長(zhǎng)Runge-Kutta算法及MATLAB實(shí)現(xiàn) 一階微分方程組的數(shù)值解 微分方程數(shù)值解的驗(yàn)證,9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.1 微分方程問(wèn)題算法概述,一階顯式的微分方程組 標(biāo)準(zhǔn)形式 其中， 狀態(tài)向量 非線性

8、函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.1.1 微分方程求解的誤差與步長(zhǎng)問(wèn)題,Euler算法： 設(shè)初始時(shí)刻系統(tǒng)狀態(tài)向量的值為 微分方程左側(cè)的導(dǎo)數(shù)近似為： 時(shí)刻微分方程的近似解為:,9/3/2020星期六, 2008-9- 6, 13:10:16,在 時(shí)刻系統(tǒng)狀態(tài)向量的值為： 簡(jiǎn)記為 故可以假設(shè)在tk時(shí)刻系統(tǒng)的狀態(tài)向量為xk在 時(shí)刻Euler算法的數(shù)值解為： h 被稱(chēng)為步長(zhǎng),9/3/2020星期六, 2008-9- 6, 13:10:16,不能無(wú)限制地減小h 的值的兩條原因： 減慢計(jì)算速度 增加累積誤差 在對(duì)微分方程求解過(guò)程中應(yīng)采取的三個(gè)措施： 選擇適當(dāng)?shù)?/p>

9、步長(zhǎng) 改進(jìn)近似算法精度 采用變步長(zhǎng)方法,9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.2 四階定步長(zhǎng)Runge-Kutta算法及MATLAB實(shí)現(xiàn),四階定步長(zhǎng)Runge-Kutta算法的數(shù)學(xué)描述：,9/3/2020星期六, 2008-9- 6, 13:10:16,其中h為計(jì)算步長(zhǎng) 下一個(gè)步長(zhǎng)的狀態(tài)變量值為： MATLAB調(diào)用格式 函數(shù)中使用了循環(huán)結(jié)構(gòu),9/3/2020星期六, 2008-9- 6, 13:10:16,構(gòu)造MATLAB函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.3 一階微分方程組的數(shù)值解,四階五級(jí)Runge-Kutta

10、-Felhberg算法 基于MATLAB的微分方程求解函數(shù) MATLAB下帶有附加參數(shù)的微分方程求解 微分方程數(shù)值解的驗(yàn)證,9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.3.1 四階五級(jí)Runge-Kutta-Felhberg算法,Runge-Kutta-Felhberg算法 假設(shè)當(dāng)前的步長(zhǎng)為hk ，定義6 個(gè)ki變量：,9/3/2020星期六, 2008-9- 6, 13:10:16,下一步的狀態(tài)向量 定義一個(gè)誤差向量：,9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.3.2 基于MATLAB的微分方程求解函數(shù),求解常微分方程的MATLA

11、B函數(shù)調(diào)用格式,9/3/2020星期六, 2008-9- 6, 13:10:16,描述需要求解的微分方程組： 修改控制變量:,9/3/2020星期六, 2008-9- 6, 13:10:16,微分方程的求解步驟,寫(xiě)出微分方程的數(shù)學(xué)形式 編輯方程的MATLAB代碼 M-函數(shù) 匿名函數(shù) Inline函數(shù)，不推薦使用 求解方程 將繪制結(jié)果 證明結(jié)果的正確性,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.8,求解下列Lorenz模型 式中參數(shù)為 初始條件為,9/3/2020星期六, 2008-9- 6, 13:10:16,方程的描述 Inline函數(shù)（不推薦使用）,9/3/

12、2020星期六, 2008-9- 6, 13:10:16,MATLAB求解命令 用comet3()繪制相空間軌跡,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.9,試用匿名函數(shù)描述例7.8的該方程，并求出該微分方程的數(shù)值解，與前面的結(jié)果進(jìn)行比較 定義該方程：,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,引入附加參數(shù)的目的 在前面的例子中，如果參數(shù)改變 ，不需要修改原函數(shù) 可以利用附帶的參數(shù) 新的函數(shù)編輯格式 新的函數(shù)調(diào)用命令,7.2.3.3 MATLAB下

13、帶有附加參數(shù)的微分方程求解,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.10,試編寫(xiě)帶有附加參數(shù)的MATLAB函數(shù)來(lái)描述Lorenz方程，求解 和一組新參數(shù) 下方程的數(shù)值解 編寫(xiě)帶附加參數(shù)的方程，并求解：,9/3/2020星期六, 2008-9- 6, 13:10:16,接上頁(yè),9/3/2020星期六, 2008-9- 6, 13:10:16,選擇新參數(shù) 新的調(diào)用命令,9/3/2020星期六, 2008-9- 6, 13:10:16,7.2.4 微分方程數(shù)值解的驗(yàn)證,仿真算法和控制參數(shù)選擇不當(dāng)，如相對(duì)誤差限，則可能得出不可信的結(jié)果，甚至是錯(cuò)誤的結(jié)果 修改仿真控制參

14、數(shù)，如可以接受的誤差限，再判斷是否和上次得出的結(jié)果一致 選擇不同的微分方程求解算法,9/3/2020星期六, 2008-9- 6, 13:10:16,7.3 微分方程轉(zhuǎn)換,單個(gè)高階常微分方程處理方法 高階常微分方程組的變換方法 矩陣微分方程的變換與求解方法,9/3/2020星期六, 2008-9- 6, 13:10:16,7.3.1 單個(gè)高階常微分方程處理方法,一個(gè)高階常微分方程的一般形式為： 輸出變量的各階導(dǎo)數(shù)初始值為： 選擇一組狀態(tài)變量：,9/3/2020星期六, 2008-9- 6, 13:10:16,原高階常微分方程模型變換為下述一階微分方程組： 初值為：,9/3/2020星期六, 2

15、008-9- 6, 13:10:16,例 7.11,Van der Pol 方程 初值為 選擇初值為 新微分方程組 MATLAB代碼：,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解： 一個(gè)反例,9/3/2020星期六, 2008-9- 6, 13:10:16,一個(gè)反例,給定如下參數(shù) 建議：不要運(yùn)行下述語(yǔ)句 應(yīng)該使用剛性微分方程的算法,9/3/2020星期六, 2008-9- 6, 13:10:16,7.3.2 高階常微分方程組的變換方法,多元高階常微分方程組的處理 狀態(tài)變量的選擇不唯一 建議：選擇如下?tīng)顟B(tài)變量,9/3/2020星期六, 2008-9- 6,

16、 13:10:16,新的狀態(tài)方程是：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.12,Apollo衛(wèi)星的運(yùn)動(dòng)軌跡(x, y)滿(mǎn)足下面的方程 其中， 并且 初始狀態(tài),9/3/2020星期六, 2008-9- 6, 13:10:16,選擇一組狀態(tài)變量： 得出一階常微分方程組： 其中， 并且,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB描述,9/3/2020星期六, 2008-9- 6, 13:10:16,求解： 改變精度：,9/3/2020星期六, 2008-9- 6, 13:10:16,繪制計(jì)算步長(zhǎng)的曲線，求解全程所采用的最小步長(zhǎng)

17、： 評(píng)注 不應(yīng)當(dāng)太依賴(lài)于默認(rèn)值 最好設(shè)置一下控制精度到 步長(zhǎng)值必須達(dá)到0.0001 對(duì)結(jié)果應(yīng)該進(jìn)行驗(yàn)證,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.13,用定步長(zhǎng)的四階Runge-Kutta算法求下式 其中，,9/3/2020星期六, 2008-9- 6, 13:10:16,選擇步長(zhǎng)為0.01： 選擇步長(zhǎng)為0.001：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.14,轉(zhuǎn)換成一階微分方程組 選擇狀態(tài)變量： 并且,9/3/2020星期六, 2008-9- 6, 13:10:16,帶入第二個(gè)等式得出 得出一階微分方程組：,9/3/2020

18、星期六, 2008-9- 6, 13:10:16,使用solve()函數(shù) 和前面完全一致的結(jié)果,9/3/2020星期六, 2008-9- 6, 13:10:16,7.3.3 矩陣微分方程的變換與求解方法,Lagrange方程 其中， 均為 矩陣 均為 矩陣 選擇狀態(tài)變量,9/3/2020星期六, 2008-9- 6, 13:10:16,寫(xiě)出系統(tǒng)的狀態(tài)方程模型可以重新寫(xiě)為 其中， 因此，矩陣方程可以通過(guò)函數(shù)ode45()或 ode15s()來(lái)求解,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.15,考慮二級(jí)倒立擺系統(tǒng)的數(shù)學(xué)模型 其中，q = a,q1,q2T，和a為小

19、車(chē)位置，q1，q2為下擺桿、上擺桿與垂直方向夾角，則矩陣 為：,9/3/2020星期六, 2008-9- 6, 13:10:16,參數(shù)為,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB描述語(yǔ)句,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解代碼： 原系統(tǒng)并不穩(wěn)定,9/3/2020星期六, 2008-9- 6, 13:10:16,若矩陣 與 無(wú)關(guān)，則該微分方程為線性微分方程，相應(yīng)的線性狀態(tài)方程模型為,9/3/2020星期六, 2008-9- 6, 13:10:16,Riccati微分方程,Riccati微分方程的數(shù)學(xué)描述 求解這

20、樣的方程同樣需要轉(zhuǎn)換成向量型一階顯式微分方程組，然后進(jìn)行求解,9/3/2020星期六, 2008-9- 6, 13:10:16,描述微分Riccati方程的MATLAB函數(shù) MATLAB求解語(yǔ)句調(diào)用格式 允許終止時(shí)間小于起始時(shí)間,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.16,若已知某微分Riccati方程中矩陣及初值如下，試求解該方程,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,7.4 特殊微分方程的數(shù)值解,剛性微分方程的求解 隱式微分方程求解 微分

21、代數(shù)方程的求解 延遲微分方程求解 切換微分方程的求解 隨機(jī)線性微分方程的求解,9/3/2020星期六, 2008-9- 6, 13:10:16,7.4.1 剛性微分方程的求解,Van der Pol方程中 與上述方程類(lèi)似的微分方程中，某些狀態(tài)的解變化緩慢，另一些變化快 這類(lèi)方程被稱(chēng)為剛性方程 而應(yīng)該使用函數(shù)ode15s()求解該類(lèi)方程，函數(shù)的調(diào)用格式和ode45()完全一致,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.17,試求解m=1000時(shí)的Van der Pol方程的數(shù)值解 MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:

22、16,例 7.18,求解該微分方程 傳統(tǒng)的教科書(shū)中，都認(rèn)為上式為剛性的 解析解：,9/3/2020星期六, 2008-9- 6, 13:10:16,該問(wèn)題的數(shù)值解：,9/3/2020星期六, 2008-9- 6, 13:10:16,使用定步長(zhǎng)算法： 為了得到更好的結(jié)果，需要減小步長(zhǎng)值，但是，計(jì)算的速度會(huì)隨之降低,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.19,對(duì)下式用數(shù)值解法求解： 初值條件 MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,變步長(zhǎng)解法： 用函數(shù)ode15s()替代函數(shù)ode45()：,9/3/2020星期

23、六, 2008-9- 6, 13:10:16,7.4.2 隱式微分方程求解,所謂隱式微分方程，就是那些不能轉(zhuǎn)換成式一階顯式常微分方程組的微分方程 可以借用顯式微分方程求解的函數(shù)來(lái)求解這類(lèi)問(wèn)題 或者使用函數(shù)ode15i(),9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.20,給定x1(0)=x2(0)=0 ，求下式的數(shù)值解 令x=x1,x2T，可以將原微分方程改寫(xiě)成矩陣形式：,9/3/2020星期六, 2008-9- 6, 13:10:16,其中 MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.21,給定x(0)=1,

24、0,0,1T，求下式的數(shù)值解 定義狀態(tài)變量： 改寫(xiě)原隱式方程：,9/3/2020星期六, 2008-9- 6, 13:10:16,用MATLAB描述原方程 求解方程：,9/3/2020星期六, 2008-9- 6, 13:10:16,使用函數(shù)ode15i(),MATLAB 7.0版本的新函數(shù)ode15i()可以直接用于隱式微分方程的求解. 隱式微分方程的數(shù)學(xué)描述為：,9/3/2020星期六, 2008-9- 6, 13:10:16,編寫(xiě)函數(shù)fun()描述該隱式微分方程，再decic()函數(shù)求出未完全定義的初值條件 函數(shù)decic()的調(diào)用格式 求解隱式微分方程的函數(shù)調(diào)用格式,9/3/2020星

25、期六, 2008-9- 6, 13:10:16,例 7.22,給定初始狀態(tài) x(0)=1,0,0,1T，用隱式微分方程求解的方法解出： 選擇狀態(tài)變量：,9/3/2020星期六, 2008-9- 6, 13:10:16,原方程可以變換成： 考慮下述形式,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,7.4.3 微分代數(shù)方程的求解,所謂微分代數(shù)方程(DAE)，是指在微分方程中，某些變量間滿(mǎn)足某些代數(shù)方程的約束 微分方程的更一般形式可以寫(xiě)成： 微分代數(shù)方程中， 矩陣為奇異矩陣 微分代數(shù)方程的

26、求解中不能使用,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.23,給定初值條件，x1(0)=0.8, x2(0)=x3(0)=0.1，求數(shù)值解： 矩陣形式表示該微分代數(shù)方程：,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句： 解此微分代數(shù)方程,9/3/2020星期六, 2008-9- 6, 13:10:16,可以將這個(gè)微分代數(shù)方程可以轉(zhuǎn)換成常微分方程求解，從約束式子中求出 代入其他兩個(gè)微分方程式子 可以寫(xiě)出匿名函數(shù)描述微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,再次求解,9/3/2020星期六

27、, 2008-9- 6, 13:10:16,通過(guò)的隱式微分方程求解函數(shù)ode15i()求解 描述隱式微分方程：,9/3/2020星期六, 2008-9- 6, 13:10:16,令 ，其中*表示自由值 解出相容的初始條件，并直接求解該微分代數(shù)方程,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.24,用微分代數(shù)方程求解的方式求解下述隱式微分方程 其中，初始條件為,9/3/2020星期六, 2008-9- 6, 13:10:16,描述微分方程與質(zhì)量矩陣 求解微分代數(shù)方程,9/3/2020星期六, 2008-9- 6, 13:10:16,7.4.4 延遲微分方程求解,延遲

28、微分方程組的一般形式為： 其中， 為狀態(tài)變量x(t)的延遲常數(shù) 隱式Runge-Kutta算法dde23(),9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.25,用數(shù)值法求解下述延遲微分方程組 初始狀態(tài)： 選擇狀態(tài)變量：,9/3/2020星期六, 2008-9- 6, 13:10:16,新的標(biāo)準(zhǔn)形式 M-函數(shù)描述方程組：,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.26,中性延遲微分方程 不能用函數(shù)dde23()來(lái)求解,9/3/2020星期六, 2

29、008-9- 6, 13:10:16,7.4.5 切換微分方程的求解,切換系統(tǒng)是控制理論中的一個(gè)重要的研究領(lǐng)域，切換系統(tǒng)就是在某種規(guī)律下其模型在多個(gè)模型間切換的系統(tǒng) 切換微分方程的數(shù)學(xué)描述,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.27,假設(shè)已知系統(tǒng)模型 若 ，切換到系統(tǒng) ，而 切換到 初始狀態(tài)為 試求解該系統(tǒng)的微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,切換系統(tǒng)的M-函數(shù)表示 由匿名函數(shù)表示切換系統(tǒng),9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9-

30、 6, 13:10:16,7.4.6 隨機(jī)線性微分方程的求解,隨機(jī)線性連續(xù)狀態(tài)方程模型為 式中， 為兼容矩陣， 為確定性輸入向量， 為Gauss白噪聲向量，滿(mǎn)足,9/3/2020星期六, 2008-9- 6, 13:10:16,定義一個(gè)變量 ， 亦為Gauss白噪聲，滿(mǎn)足 其中， 為一個(gè)協(xié)方差矩陣 有下式成立 狀態(tài)變量的解析解可以寫(xiě)成：,9/3/2020星期六, 2008-9- 6, 13:10:16,離散化解析解，則有 式中，,9/3/2020星期六, 2008-9- 6, 13:10:16,且有 其中， 對(duì)V進(jìn)行Taylor冪級(jí)數(shù)展開(kāi),9/3/2020星期六, 2008-9- 6, 13:

31、10:16,其中， 與 可以由下式遞推地求出 遞推初值為 由Cholesky分解 ，且,9/3/2020星期六, 2008-9- 6, 13:10:16,其中， 為 向量，且有各個(gè)分量 系統(tǒng)的離散形式的遞推解可以寫(xiě)成,9/3/2020星期六, 2008-9- 6, 13:10:16,構(gòu)造隨機(jī)輸入下連續(xù)線性系統(tǒng)離散化函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,函數(shù)的調(diào)用格式為 其中，G為系統(tǒng)模型，s輸入信號(hào)協(xié)方差矩陣，Dt為采樣周期，(F,G,D,C)為離散化狀態(tài)方程的相應(yīng)矩陣,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.28,受控對(duì)象的

32、傳遞函數(shù)模型為 用白噪聲信號(hào)激勵(lì)該系統(tǒng)，試對(duì)其進(jìn)行仿真分析并得出輸出信號(hào)的統(tǒng)計(jì)規(guī)律,9/3/2020星期六, 2008-9- 6, 13:10:16,得出離散化的狀態(tài)方程模型 進(jìn)行仿真,9/3/2020星期六, 2008-9- 6, 13:10:16,分析統(tǒng)計(jì)規(guī)律，疊印理論概率密度函數(shù)曲線,9/3/2020星期六, 2008-9- 6, 13:10:16,7.5 邊值問(wèn)題的計(jì)算機(jī)求解,線性方程邊值問(wèn)題的打靶算法 非線性方程邊值問(wèn)題的打靶算法 一般邊值微分方程的求解方法,9/3/2020星期六, 2008-9- 6, 13:10:16,二階微分方程的邊值問(wèn)題的數(shù)學(xué)描述為： 假設(shè)想在區(qū)間a, b上

33、研究該方程的解，且已知在這兩個(gè)邊界點(diǎn)上滿(mǎn)足 上面的方程稱(chēng)為邊界a和b的邊界條件,9/3/2020星期六, 2008-9- 6, 13:10:16,7.5.1 線性方程邊值問(wèn)題的打靶算法,給定下列微分方程 其中， 均為給定函數(shù)，且可以得出 取變量 ，得一階微分方程組,9/3/2020星期六, 2008-9- 6, 13:10:16,運(yùn)用打靶算法求解上述方程 求出下面方程初值問(wèn)題的數(shù)值解 求出下面方程初值問(wèn)題的數(shù)值解 求出下面方程初值問(wèn)題的數(shù)值解,9/3/2020星期六, 2008-9- 6, 13:10:16,若 ， 則計(jì)算 計(jì)算下面初值問(wèn)題的數(shù)值解，則 即為原邊值問(wèn)題的數(shù)值解,9/3/2020

34、星期六, 2008-9- 6, 13:10:16,構(gòu)造MATLAB函數(shù)實(shí)現(xiàn)上述算法,9/3/2020星期六, 2008-9- 6, 13:10:16,該函數(shù)中，tspan為初始和終止仿真時(shí)間構(gòu)成的向量； 為邊界值 定義兩個(gè)輔助函數(shù)f1()和f2()來(lái)描述原模型 f2()描述 f1()描述原來(lái)方程的齊次,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.29,試求下述線性微分方程在0,1段方程的數(shù)值解 編輯出下面的兩個(gè)匿名函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句： 由微分方程理論得出原方程的解析解為： 結(jié)果檢驗(yàn)：,9/3

35、/2020星期六, 2008-9- 6, 13:10:16,7.5.2非線性方程邊值問(wèn)題的打靶算法,假定原始問(wèn)題可以轉(zhuǎn)換成下面的初值問(wèn)題 則問(wèn)題轉(zhuǎn)換成求解 Newton迭代法求取m,9/3/2020星期六, 2008-9- 6, 13:10:16,式中 且可以由下面的微分方程初值問(wèn)題來(lái)求解,9/3/2020星期六, 2008-9- 6, 13:10:16,其中，要求能顯式地求出下式： 構(gòu)造MATLAB函數(shù)實(shí)現(xiàn)上述算法：,9/3/2020星期六, 2008-9- 6, 13:10:16,試求解如下非線性微分方程邊值問(wèn)題 推導(dǎo)出偏導(dǎo)數(shù) 得出：,例 7.30,9/3/2020星期六, 2008-9-

36、 6, 13:10:16,構(gòu)造MATLAB匿名函數(shù)： MATALB求解語(yǔ)句： 解析解為： 驗(yàn)證：,9/3/2020星期六, 2008-9- 6, 13:10:16,7.5.3一般邊值微分方程的求解方法,bvp5c() 函數(shù)求解出一般邊值微分方程 參數(shù)初始化 微分方程和邊值問(wèn)題的MATLAB函數(shù)描述和邊值問(wèn)題的求解,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.31,試用bvp5c()數(shù)重新求解下述邊值問(wèn)題 選擇狀態(tài)變量 狀態(tài)方程變?yōu)?9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句： 注意最后一條命令,9/3/2020星期六,

37、2008-9- 6, 13:10:16,若邊值問(wèn)題修改成 則可以如下修改f2()，并求解該方程,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.32,已知某常微分方程模型為 且已知 試求出 并求解本微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,令 描述系統(tǒng)的動(dòng)態(tài)模型和邊值問(wèn)題,9/3/2020星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句： 初值向量 選擇若不當(dāng)，可能使得求解過(guò)程中的Jacobian矩陣奇異，若出現(xiàn)此現(xiàn)象，選擇其他的初值,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6 偏微

38、分方程求解入門(mén),偏微分方程組求解 二階偏微分方程的數(shù)學(xué)描述 偏微分方程的求解界面應(yīng)用舉例,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.1 偏微分方程組求解,MATLAB語(yǔ)言提供的pdepe()函數(shù)，可以直接求解以下類(lèi)型的偏微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,偏微分方程可以編寫(xiě)下面的函數(shù)描述，其入口為 其中，pdefun為函數(shù)名。這樣，由給定的輸入變量即可計(jì)算出 這3個(gè)函數(shù),9/3/2020星期六, 2008-9- 6, 13:10:16,邊界條件可以用下面的函數(shù)描述為 初始條件函數(shù) 簡(jiǎn)單函數(shù)描述初始條件，函數(shù)調(diào)用格式 求解偏

39、微分方程的函數(shù)調(diào)用格式,9/3/2020星期六, 2008-9- 6, 13:10:16,試求解下面的偏微分方程 其中，,例 7.33,9/3/2020星期六, 2008-9- 6, 13:10:16,初始條件： 邊界條件：,9/3/2020星期六, 2008-9- 6, 13:10:16,將原方程改寫(xiě)為如下的形式： 可見(jiàn)， ，且,9/3/2020星期六, 2008-9- 6, 13:10:16,描述偏微分方程的M-函數(shù)： 寫(xiě)出邊值方程： 左邊界： 右邊界：,9/3/2020星期六, 2008-9- 6, 13:10:16,描述邊界條件的函數(shù)： 描述初值條件的函數(shù)：,9/3/2020星期六,

40、2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.2 二階偏微分方程的數(shù)學(xué)描述,橢圓型偏微分方程 拋物線型偏微分方程 雙曲型偏微分方程 特征值型偏微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.2.1 橢圓型偏微分方程,橢圓型偏微分方程的一般表示形式為： 其中， 梯度： 散度：,9/3/2020星期六, 2008-9- 6, 13:10:16,若c為常數(shù)，上式可以化簡(jiǎn)為： 橢圓型偏微分方程可以更簡(jiǎn)單地寫(xiě)成：,9/3/2020星期六, 2008-9- 6, 13:10:16

41、,拋物型偏微分方程的一般形式為： 若c為常數(shù)，則該方程可以更簡(jiǎn)單地寫(xiě)成,7.6.2.2 拋物線型偏微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,雙曲型偏微分方程的一般形式為： 若c為常數(shù)，則該方程可以更簡(jiǎn)單地寫(xiě)成,7.6.2.3 雙曲型偏微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,特征值型偏微分方程的一般形式為： 若c為常數(shù)，則該方程可以更簡(jiǎn)單地寫(xiě)成,7.6.2.4 特征值型偏微分方程,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.3 偏微分方程的求解界面應(yīng)用舉例,偏微分方程求解程序概述 偏微分方程求解區(qū)

42、域繪制 偏微分方程邊界條件描述 偏微分方程求解舉例 時(shí)變解的動(dòng)畫(huà)顯示 函數(shù)參數(shù)的偏微分方程求解,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.3.1 偏微分方程求解程序概述,偏微分方程工具箱概述 啟動(dòng)偏微分方程工具箱求解界面GUI 在MATLAB下鍵入pdetool 該界面分為四個(gè)部分 菜單系統(tǒng) 工具欄 集合編輯 求解區(qū)域,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.3.2 偏微分方程求解區(qū)域繪制,繪制并且定義偏微分方程求解區(qū)域 繪制好求解區(qū)域后，單擊工具欄中 按鈕就可以定義求解區(qū)域 單擊D按鈕將求解區(qū)域用三角形劃分，單擊右側(cè)的按鈕加密

43、網(wǎng)格；網(wǎng)格越密，計(jì)算的結(jié)果越精確，但計(jì)算時(shí)間則越長(zhǎng),9/3/2020星期六, 2008-9- 6, 13:10:16,Dirichlet條件： 其中， 表示求解區(qū)域的邊界 Neumann條件： 其中， 為 向量法向的偏導(dǎo)數(shù),7.6.3.3 偏微分方程邊界條件描述,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.3.4 偏微分方程求解舉例,設(shè)置了求解區(qū)域和邊界條件，并選擇了合適的偏微分方程后，就可以單擊工具欄的等號(hào)按鈕 = 立即得出微分方程的解 微分方程的結(jié)果可以用其他很多方式顯示，如繪制等值線圖，引力線圖，網(wǎng)格型三維圖形,9/3/2020星期六, 2008-9- 6,

44、 13:10:16,例 7.34,試求解雙曲型偏微分方程：,9/3/2020星期六, 2008-9- 6, 13:10:16,設(shè)置時(shí)間向量 選擇其中的動(dòng)畫(huà)Animation選項(xiàng) 設(shè)置動(dòng)畫(huà)的播放速度Options 動(dòng)畫(huà)輸出到MATLAB工作空間 在MATLAB圖形窗口中播放得出的動(dòng)畫(huà),7.6.3.5 時(shí)變解的動(dòng)畫(huà)顯示,9/3/2020星期六, 2008-9- 6, 13:10:16,7.6.3.6 函數(shù)參數(shù)的偏微分方程求解,偏微分方程工具箱能處理含有非線性系數(shù)的橢圓偏微分方程問(wèn)題 使用x，y表示微分方程中 或 ux和uy表示 和,9/3/2020星期六, 2008-9- 6, 13:10:16,

45、假設(shè)偏微分方程為 它滿(mǎn)足橢圓型偏微分方程 其中邊界為u=0，試求解該偏微分方程,例 7.35,9/3/2020星期六, 2008-9- 6, 13:10:16,7.7 微分方程的框圖求解,Simulink簡(jiǎn)介 Simulink相關(guān)模塊 微分方程的Simulink建模與求解,9/3/2020星期六, 2008-9- 6, 13:10:16,7.7.1 Simulink簡(jiǎn)介,1990 年前后出現(xiàn)最早的 Simulink，當(dāng)時(shí)名為SimuLAB，1992 年改為 Simulink Simulink 的名字有兩重含義 仿真 (simu) 與模型/模塊 連接 (link) 使用odegroup命令可以打

46、開(kāi)自定義模塊集.,9/3/2020星期六, 2008-9- 6, 13:10:16,7.7.2 Simulink相關(guān)模塊,常用的模塊： 輸入輸出端口(In1, out1) 時(shí)鐘模塊：產(chǎn)生時(shí)間變量t 常用輸入模塊： Sine模塊產(chǎn)生正弦信號(hào) Step模塊產(chǎn)生階躍信號(hào) Constant模塊產(chǎn)生恒值信號(hào) 積分器模塊(Int)：描述一階導(dǎo)數(shù)，搭建積分的動(dòng)態(tài)系統(tǒng) 延遲模塊(Transport Delay),9/3/2020星期六, 2008-9- 6, 13:10:16,增益模塊： Gain：主要用于信號(hào)的放大 Sliding Gain：滾動(dòng)桿，實(shí)時(shí)地改變模塊的增益 Matrix Gain：矩陣增益模塊

47、 數(shù)學(xué)運(yùn)算模塊 +, -, *, 數(shù)學(xué)函數(shù)模塊 Sin, exp 信號(hào)向量化模塊 Mux：將若干路信號(hào)混成向量型的信號(hào) DeMux：將向量型信號(hào)解出單路的信號(hào),9/3/2020星期六, 2008-9- 6, 13:10:16,7.7.3 微分方程的Simulink建模與求解,建立起相應(yīng)的微分方程的 Simulink 模型 可以用sim()函數(shù)對(duì)其模型直接求解，該函數(shù)的調(diào)用格式 得出微分方程的數(shù)值解,9/3/2020星期六, 2008-9- 6, 13:10:16,例 7.36,求解Lorenz方程的問(wèn)題，其中，b=8/3, r=10, s=28，各個(gè)狀態(tài)變量的初值為 x1(0)= x2(0)= 0, x3(0)=10-10，試用Simulink搭建該模型并得出仿真結(jié)果.,9/3/2020星期六, 2008-9- 6, 13:10:16,搭建Simulink框圖