《高等應(yīng)用數(shù)學(xué)問(wèn)題的matlab求解-第二版》math-chap07

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

2、入門(mén)微分方程的框圖求解,3/19/2024星期六, 2008-9- 6, 13:10:16,7.1 常系數(shù)線性微分方程的解析解方法,線性常系數(shù)微分方程解析解的數(shù)學(xué)描述微分方程的解析解方法Laplace變換在線性微分方程求解中的應(yīng)用線性狀態(tài)空間方程的解析解特殊非線性微分方程的解析解,3/19/2024星期六, 2008-9- 6, 13:10:16,7.1.1 線性常系數(shù)微分方程解析解的數(shù)學(xué)描述,常系數(shù)線性微分方程

3、的一般描述方法為其中， 均為常數(shù),3/19/2024星期六, 2008-9- 6, 13:10:16,對(duì)零初值問(wèn)題：對(duì)應(yīng)得出Laplace變換為設(shè)代數(shù)方程的特征根si均相異，則解析解為 其中Ci為待定系數(shù)，并且g(t)是滿足u(t)輸入的一個(gè)特解,3/19/2024星期六, 2008-9- 6, 13:10:16,簡(jiǎn)單形式，自變量設(shè)為t指明自變量x,7.1.2 微分方程

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

5、/2024星期六, 2008-9- 6, 13:10:16,最終的近似結(jié)果為：,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.2,給定輸入信號(hào) ，求解 其中，,3/19/2024星期六, 2008-9- 6, 13:10:16,結(jié)果：,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.3,試求解線性微分方程組的解析解

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

7、tcos(2t+1)+5，假定所有的初始條件都為0，試求出下述微分方程的解析解,3/19/2024星期六, 2008-9- 6, 13:10:16,用Laplace變換法,3/19/2024星期六, 2008-9- 6, 13:10:16,在時(shí)域中輸出信號(hào)的解析解：,3/19/2024星期六, 2008-9- 6, 13:10:16,使用dsolve()函數(shù),3/19/2024星期六, 2008-9- 6, 13:10:

8、16,繪制解曲線：解析解相同,3/19/2024星期六, 2008-9- 6, 13:10:16,7.1.4 線性狀態(tài)空間方程的解析解,假設(shè)線性狀態(tài)空間模型的一般表示為 其中，A,B,C,D是常數(shù)矩陣，且已知狀態(tài)向量初值 ，該方程的解析解是：,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.5,給定輸入信號(hào)為u(t) = 2 + 2e-3t sin 2t，求出下面矩陣

9、描述的狀態(tài)空間方程的解析解,3/19/2024星期六, 2008-9- 6, 13:10:16,直接積分法求出結(jié)果,3/19/2024星期六, 2008-9- 6, 13:10:16,7.1.5 特殊非線性微分方程的解析解,只有少數(shù)非線性微分方程可以通過(guò)dsolve()函數(shù)得出解析解通過(guò)例子演示非線性方程的解析解求解 同時(shí)還將演示不能求解的例子,3/19/2024星期六, 2008-9- 6, 13:10:16

10、,例 7.6,求解非線性微分方程：改變?cè)⒎址匠痰男问?3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.7,試求出著名的Van der Pol方程的解析解嘗試如下的MATLAB命令,3/19/2024星期六, 2008-9- 6, 13:10:16,7.2 微分方程問(wèn)題的數(shù)值解法,微分方程問(wèn)題算法概述四階定步長(zhǎng)Runge-Kutta算法及MATLAB實(shí)現(xiàn)一階微分方程組的數(shù)值解

11、微分方程數(shù)值解的驗(yàn)證,3/19/2024星期六, 2008-9- 6, 13:10:16,7.2.1 微分方程問(wèn)題算法概述,一階顯式的微分方程組標(biāo)準(zhǔn)形式 其中， 狀態(tài)向量非線性函數(shù),3/19/2024星期六, 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ù)近似為：

12、 時(shí)刻微分方程的近似解為:,3/19/2024星期六, 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),3/19/2024星期六, 2008-9- 6, 13:10:16,不能無(wú)限制地減小h 的值的兩條原因：減慢計(jì)算速度增加累

13、積誤差在對(duì)微分方程求解過(guò)程中應(yīng)采取的三個(gè)措施：選擇適當(dāng)?shù)牟介L(zhǎng)改進(jìn)近似算法精度采用變步長(zhǎng)方法,3/19/2024星期六, 2008-9- 6, 13:10:16,7.2.2 四階定步長(zhǎng)Runge-Kutta算法及MATLAB實(shí)現(xiàn),四階定步長(zhǎng)Runge-Kutta算法的數(shù)學(xué)描述：,3/19/2024星期六, 2008-9- 6, 13:10:16,其中h為計(jì)算步長(zhǎng)下一個(gè)步長(zhǎng)的狀態(tài)變量值為：MATLAB調(diào)用格式

14、函數(shù)中使用了循環(huán)結(jié)構(gòu),3/19/2024星期六, 2008-9- 6, 13:10:16,構(gòu)造MATLAB函數(shù),3/19/2024星期六, 2008-9- 6, 13:10:16,7.2.3 一階微分方程組的數(shù)值解,四階五級(jí)Runge-Kutta-Felhberg算法基于MATLAB的微分方程求解函數(shù)MATLAB下帶有附加參數(shù)的微分方程求解微分方程數(shù)值解的驗(yàn)證,3/19/2024星期六, 2008-9- 6

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

16、分方程的MATLAB函數(shù)調(diào)用格式,3/19/2024星期六, 2008-9- 6, 13:10:16,描述需要求解的微分方程組：修改控制變量:,3/19/2024星期六, 2008-9- 6, 13:10:16,微分方程的求解步驟,寫(xiě)出微分方程的數(shù)學(xué)形式編輯方程的MATLAB代碼M-函數(shù)匿名函數(shù)Inline函數(shù)，不推薦使用求解方程將繪制結(jié)果證明結(jié)果的正確性,3/19/2024星期六, 2008-9-

17、6, 13:10:16,例 7.8,求解下列Lorenz模型 式中參數(shù)為初始條件為,3/19/2024星期六, 2008-9- 6, 13:10:16,方程的描述Inline函數(shù)（不推薦使用）,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解命令用comet3()繪制相空間軌跡,3/19/2024星期六, 2008-9- 6, 13:10:16,例

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

19、 MATLAB下帶有附加參數(shù)的微分方程求解,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.10,試編寫(xiě)帶有附加參數(shù)的MATLAB函數(shù)來(lái)描述Lorenz方程，求解 和一組新參數(shù) 下方程的數(shù)值解編寫(xiě)帶附加參數(shù)的方程，并求解：,3/19/2024星期六, 20

20、08-9- 6, 13:10:16,接上頁(yè),3/19/2024星期六, 2008-9- 6, 13:10:16,選擇新參數(shù)新的調(diào)用命令,3/19/2024星期六, 2008-9- 6, 13:10:16,7.2.4 微分方程數(shù)值解的驗(yàn)證,仿真算法和控制參數(shù)選擇不當(dāng)，如相對(duì)誤差限，則可能得出不可信的結(jié)果，甚至是錯(cuò)誤的結(jié)果修改仿真控制參數(shù)，如可以接受的誤差限，再判斷是否和上次得出的結(jié)果一致選擇不同的微分方程求解算法,

21、3/19/2024星期六, 2008-9- 6, 13:10:16,7.3 微分方程轉(zhuǎn)換,單個(gè)高階常微分方程處理方法高階常微分方程組的變換方法矩陣微分方程的變換與求解方法,3/19/2024星期六, 2008-9- 6, 13:10:16,7.3.1 單個(gè)高階常微分方程處理方法,一個(gè)高階常微分方程的一般形式為：輸出變量的各階導(dǎo)數(shù)初始值為：選擇一組狀態(tài)變量：,3/19/2024星期六, 2008-9-

22、 6, 13:10:16,原高階常微分方程模型變換為下述一階微分方程組：初值為：,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.11,Van der Pol 方程初值為 選擇初值為新微分方程組MATLAB代碼：,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解：一個(gè)反例,3/19/2024星期六, 2008-9

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

24、8-9- 6, 13:10:16,例 7.12,Apollo衛(wèi)星的運(yùn)動(dòng)軌跡(x, y)滿足下面的方程其中， 并且初始狀態(tài),3/19/2024星期六, 2008-9- 6, 13:10:16,選擇一組狀態(tài)變量：得出一階常微分方程組：其中，并且,3/19/2024星期六, 2008-9- 6, 13:10:16,M

25、ATLAB描述,3/19/2024星期六, 2008-9- 6, 13:10:16,求解：改變精度：,3/19/2024星期六, 2008-9- 6, 13:10:16,繪制計(jì)算步長(zhǎng)的曲線，求解全程所采用的最小步長(zhǎng)：評(píng)注不應(yīng)當(dāng)太依賴(lài)于默認(rèn)值最好設(shè)置一下控制精度到步長(zhǎng)值必須達(dá)到0.0001對(duì)結(jié)果應(yīng)該進(jìn)行驗(yàn)證,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.13,用定步長(zhǎng)的四

26、階Runge-Kutta算法求下式其中，,3/19/2024星期六, 2008-9- 6, 13:10:16,選擇步長(zhǎng)為0.01：選擇步長(zhǎng)為0.001：,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.14,轉(zhuǎn)換成一階微分方程組選擇狀態(tài)變量： 并且,3/19/2024星期六, 2008-9- 6, 13:10:16,帶入第二個(gè)等式得出得出一階微分方程組：,3

27、/19/2024星期六, 2008-9- 6, 13:10:16,使用solve()函數(shù)和前面完全一致的結(jié)果,3/19/2024星期六, 2008-9- 6, 13:10:16,7.3.3 矩陣微分方程的變換與求解方法,Lagrange方程其中， 均為 矩陣 均為 矩陣選擇狀態(tài)變量,3/19/2024星期六, 200

28、8-9- 6, 13:10:16,寫(xiě)出系統(tǒng)的狀態(tài)方程模型可以重新寫(xiě)為其中，因此，矩陣方程可以通過(guò)函數(shù)ode45()或 ode15s()來(lái)求解,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.15,考慮二級(jí)倒立擺系統(tǒng)的數(shù)學(xué)模型其中，q = [a,q1,q2]T，和a為小車(chē)位置，q1，q2為下擺桿、上擺桿與垂直方向夾角，則矩陣 為：,3/19/2024星期六, 200

29、8-9- 6, 13:10:16,參數(shù)為,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB描述語(yǔ)句,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解代碼：原系統(tǒng)并不穩(wěn)定,3/19/2024星期六, 2008-9- 6, 13:10:16,若矩陣 與 無(wú)關(guān)，則該微分方程為線性微分方程，相應(yīng)的

30、線性狀態(tài)方程模型為,3/19/2024星期六, 2008-9- 6, 13:10:16,Riccati微分方程,Riccati微分方程的數(shù)學(xué)描述求解這樣的方程同樣需要轉(zhuǎn)換成向量型一階顯式微分方程組，然后進(jìn)行求解,3/19/2024星期六, 2008-9- 6, 13:10:16,描述微分Riccati方程的MATLAB函數(shù)MATLAB求解語(yǔ)句調(diào)用格式允許終止時(shí)間小于起始時(shí)間,3/19/2024星期六,

31、 2008-9- 6, 13:10:16,例 7.16,若已知某微分Riccati方程中矩陣及初值如下，試求解該方程,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,3/19/2024星期六, 2008-9- 6, 13:10:16,7.4 特殊微分方程的數(shù)值解,剛性微分方程的求解隱式微分方程求解微分代數(shù)方程的求解延遲微分方程求解切換微分方程的求解隨機(jī)線性微分方程的求解

32、,3/19/2024星期六, 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()完全一致,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.17,試求解m=1000時(shí)的Van der P

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

34、之降低,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.19,對(duì)下式用數(shù)值解法求解：初值條件MATLAB求解語(yǔ)句：,3/19/2024星期六, 2008-9- 6, 13:10:16,變步長(zhǎng)解法：用函數(shù)ode15s()替代函數(shù)ode45()：,3/19/2024星期六, 2008-9- 6, 13:10:16,7.4.2 隱式微分方程求解,所謂隱式微分方程，就是那些不能轉(zhuǎn)換成式一

35、階顯式常微分方程組的微分方程可以借用顯式微分方程求解的函數(shù)來(lái)求解這類(lèi)問(wèn)題或者使用函數(shù)ode15i(),3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.20,給定x1(0)=x2(0)=0 ，求下式的數(shù)值解令x=[x1,x2]T，可以將原微分方程改寫(xiě)成矩陣形式：,3/19/2024星期六, 2008-9- 6, 13:10:16,其中MATLAB求解語(yǔ)句：,3/19/2024星

36、期六, 2008-9- 6, 13:10:16,例 7.21,給定x(0)=[1,0,0,1]T，求下式的數(shù)值解 定義狀態(tài)變量： 改寫(xiě)原隱式方程：,3/19/2024星期六, 2008-9- 6, 13:10:16,用MATLAB描述原方程求解方程：,3/19/2024星期六, 2008-9- 6, 13:10:16,使用函數(shù)ode15i(),MATLAB 7.0版本的新函數(shù)ode15i()

37、可以直接用于隱式微分方程的求解.隱式微分方程的數(shù)學(xué)描述為：,3/19/2024星期六, 2008-9- 6, 13:10:16,編寫(xiě)函數(shù)fun()描述該隱式微分方程，再decic()函數(shù)求出未完全定義的初值條件函數(shù)decic()的調(diào)用格式求解隱式微分方程的函數(shù)調(diào)用格式,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.22,給定初始狀態(tài) x(0)=[1,0,0,1]T，用隱式微分方程求解的方

38、法解出：選擇狀態(tài)變量：,3/19/2024星期六, 2008-9- 6, 13:10:16,原方程可以變換成：考慮下述形式,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,3/19/2024星期六, 2008-9- 6, 13:10:16,7.4.3 微分代數(shù)方程的求解,所謂微分代數(shù)方程(DAE)，是指在微分方程中，某些變量間滿足某些代數(shù)方程的約束微分方程

39、的更一般形式可以寫(xiě)成：微分代數(shù)方程中， 矩陣為奇異矩陣微分代數(shù)方程的求解中不能使用,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.23,給定初值條件，x1(0)=0.8, x2(0)=x3(0)=0.1，求數(shù)值解：矩陣形式表示該微分代數(shù)方程：,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：解此微分代數(shù)方程

40、,3/19/2024星期六, 2008-9- 6, 13:10:16,可以將這個(gè)微分代數(shù)方程可以轉(zhuǎn)換成常微分方程求解，從約束式子中求出 代入其他兩個(gè)微分方程式子可以寫(xiě)出匿名函數(shù)描述微分方程,3/19/2024星期六, 2008-9- 6, 13:10:16,再次求解,3/19/2024星期六, 2008-9- 6, 13:10:16,通過(guò)的隱式微分方程求解函數(shù)ode15i()求解描述隱式微分方程：,3/19

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

42、求解微分代數(shù)方程,3/19/2024星期六, 2008-9- 6, 13:10:16,7.4.4 延遲微分方程求解,延遲微分方程組的一般形式為：其中， 為狀態(tài)變量x(t)的延遲常數(shù)隱式Runge-Kutta算法dde23(),3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.25,用數(shù)值法求解下述延遲微分方程組初始狀態(tài)：選擇狀態(tài)變量：,3/19/2024星期六,

43、 2008-9- 6, 13:10:16,新的標(biāo)準(zhǔn)形式M-函數(shù)描述方程組：,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.26,中性延遲微分方程不能用函數(shù)dde23()來(lái)求解,3/19/2024星期六, 2008-9- 6, 13:10:16,7.4.5 切換微分方程的求解

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

45、16,切換系統(tǒng)的M-函數(shù)表示由匿名函數(shù)表示切換系統(tǒng),3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,3/19/2024星期六, 2008-9- 6, 13:10:16,7.4.6 隨機(jī)線性微分方程的求解,隨機(jī)線性連續(xù)狀態(tài)方程模型為式中， 為兼容矩陣， 為確定性輸入向量， 為Gauss白噪聲向量，滿足,3/19/2024星期

46、六, 2008-9- 6, 13:10:16,定義一個(gè)變量 ， 亦為Gauss白噪聲，滿足其中， 為一個(gè)協(xié)方差矩陣有下式成立狀態(tài)變量的解析解可以寫(xiě)成：,3/19/2024星期六, 2008-9- 6, 13:10:16,離散化解析解，則有式中，,3/19/2024星期六, 2008-9- 6, 13

47、:10:16,且有其中，對(duì)V進(jìn)行Taylor冪級(jí)數(shù)展開(kāi),3/19/2024星期六, 2008-9- 6, 13:10:16,其中， 與 可以由下式遞推地求出遞推初值為由Cholesky分解 ，且,3/19/2024星期六, 2008-9- 6, 13:10:16,其中，

48、 為 向量，且有各個(gè)分量系統(tǒng)的離散形式的遞推解可以寫(xiě)成,3/19/2024星期六, 2008-9- 6, 13:10:16,構(gòu)造隨機(jī)輸入下連續(xù)線性系統(tǒng)離散化函數(shù),3/19/2024星期六, 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)矩

49、陣,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.28,受控對(duì)象的傳遞函數(shù)模型為用白噪聲信號(hào)激勵(lì)該系統(tǒng)，試對(duì)其進(jìn)行仿真分析并得出輸出信號(hào)的統(tǒng)計(jì)規(guī)律,3/19/2024星期六, 2008-9- 6, 13:10:16,得出離散化的狀態(tài)方程模型進(jìn)行仿真,3/19/2024星期六, 2008-9- 6, 13:10:16,分析統(tǒng)計(jì)規(guī)律，疊印理論概率密度函數(shù)曲線,3/19/2024星期

50、六, 2008-9- 6, 13:10:16,7.5 邊值問(wèn)題的計(jì)算機(jī)求解,線性方程邊值問(wèn)題的打靶算法非線性方程邊值問(wèn)題的打靶算法一般邊值微分方程的求解方法,3/19/2024星期六, 2008-9- 6, 13:10:16,二階微分方程的邊值問(wèn)題的數(shù)學(xué)描述為： 假設(shè)想在區(qū)間[a, b]上研究該方程的解，且已知在這兩個(gè)邊界點(diǎn)上滿足 上面的方程稱(chēng)為邊界a和b的邊界條件,3/19/2024星期六,

51、2008-9- 6, 13:10:16,7.5.1 線性方程邊值問(wèn)題的打靶算法,給定下列微分方程其中， 均為給定函數(shù)，且可以得出取變量 ，得一階微分方程組,3/19/2024星期六, 2008-9- 6, 13:10:16,運(yùn)用打靶算法求解上述方程求出下面方程初值問(wèn)題的數(shù)值解求出下面方程初值問(wèn)題的數(shù)值解

52、求出下面方程初值問(wèn)題的數(shù)值解,3/19/2024星期六, 2008-9- 6, 13:10:16,若 ， 則計(jì)算計(jì)算下面初值問(wèn)題的數(shù)值解，則 即為原邊值問(wèn)題的數(shù)值解,3/19/2024星期六, 2008-9- 6, 13:10:16,構(gòu)造MATLAB函數(shù)實(shí)現(xiàn)上述算法,3/19/2024星期六, 2008-9- 6, 13:10:16,該函數(shù)中，tspan為初始

53、和終止仿真時(shí)間構(gòu)成的向量； 為邊界值定義兩個(gè)輔助函數(shù)f1()和f2()來(lái)描述原模型f2()描述f1()描述原來(lái)方程的齊次,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.29,試求下述線性微分方程在[0,1]段方程的數(shù)值解編輯出下面的兩個(gè)匿名函數(shù),3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：由微分

54、方程理論得出原方程的解析解為：結(jié)果檢驗(yàn)：,3/19/2024星期六, 2008-9- 6, 13:10:16,7.5.2非線性方程邊值問(wèn)題的打靶算法,假定原始問(wèn)題可以轉(zhuǎn)換成下面的初值問(wèn)題則問(wèn)題轉(zhuǎn)換成求解Newton迭代法求取m,3/19/2024星期六, 2008-9- 6, 13:10:16,式中且可以由下面的微分方程初值問(wèn)題來(lái)求解,3/19/2024星期六, 2008-9- 6, 13:10:16,

55、其中，要求能顯式地求出下式：構(gòu)造MATLAB函數(shù)實(shí)現(xiàn)上述算法：,3/19/2024星期六, 2008-9- 6, 13:10:16,試求解如下非線性微分方程邊值問(wèn)題推導(dǎo)出偏導(dǎo)數(shù)得出：,例 7.30,3/19/2024星期六, 2008-9- 6, 13:10:16,構(gòu)造MATLAB匿名函數(shù)：MATALB求解語(yǔ)句：解析解為：驗(yàn)證：,3/19/2024星期六, 2008-9- 6, 13:

56、10:16,7.5.3一般邊值微分方程的求解方法,bvp5c() 函數(shù)求解出一般邊值微分方程參數(shù)初始化微分方程和邊值問(wèn)題的MATLAB函數(shù)描述和邊值問(wèn)題的求解,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.31,試用bvp5c()數(shù)重新求解下述邊值問(wèn)題選擇狀態(tài)變量狀態(tài)方程變?yōu)?3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：

57、注意最后一條命令,3/19/2024星期六, 2008-9- 6, 13:10:16,若邊值問(wèn)題修改成則可以如下修改f2()，并求解該方程,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.32,已知某常微分方程模型為且已知試求出 并求解本微分方程,3/19/2024星期六, 2008-9- 6, 13:10:16,令 描述系統(tǒng)的動(dòng)態(tài)模型和邊值問(wèn)題,3/19/

58、2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：初值向量 選擇若不當(dāng)，可能使得求解過(guò)程中的Jacobian矩陣奇異，若出現(xiàn)此現(xiàn)象，選擇其他的初值,3/19/2024星期六, 2008-9- 6, 13:10:16,7.6 偏微分方程求解入門(mén),偏微分方程組求解二階偏微分方程的數(shù)學(xué)描述偏微分方程的求解界面應(yīng)用舉例,3/19/2024星期六, 2008

59、-9- 6, 13:10:16,7.6.1 偏微分方程組求解,MATLAB語(yǔ)言提供的pdepe()函數(shù)，可以直接求解以下類(lèi)型的偏微分方程,3/19/2024星期六, 2008-9- 6, 13:10:16,偏微分方程可以編寫(xiě)下面的函數(shù)描述，其入口為其中，pdefun為函數(shù)名。這樣，由給定的輸入變量即可計(jì)算出 這3個(gè)函數(shù),3/19/2024星期六, 2008-9- 6, 13:10:16,邊界條件可以

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

61、可見(jiàn)， ，且,3/19/2024星期六, 2008-9- 6, 13:10:16,描述偏微分方程的M-函數(shù)：寫(xiě)出邊值方程：左邊界：右邊界：,3/19/2024星期六, 2008-9- 6, 13:10:16,描述邊界條件的函數(shù)：描述初值條件的函數(shù)：,3/19/2024星期六, 2008-9- 6, 13:10:16,MATLAB求解語(yǔ)句：,3/19/2024星期六, 20

62、08-9- 6, 13:10:16,7.6.2 二階偏微分方程的數(shù)學(xué)描述,橢圓型偏微分方程拋物線型偏微分方程雙曲型偏微分方程特征值型偏微分方程,3/19/2024星期六, 2008-9- 6, 13:10:16,7.6.2.1 橢圓型偏微分方程,橢圓型偏微分方程的一般表示形式為：其中，梯度：散度：,3/19/2024星期六, 2008-9- 6, 13:10:16,若c為常數(shù)，上式可以化簡(jiǎn)為

63、：橢圓型偏微分方程可以更簡(jiǎn)單地寫(xiě)成：,3/19/2024星期六, 2008-9- 6, 13:10:16,拋物型偏微分方程的一般形式為：若c為常數(shù)，則該方程可以更簡(jiǎn)單地寫(xiě)成,7.6.2.2 拋物線型偏微分方程,3/19/2024星期六, 2008-9- 6, 13:10:16,雙曲型偏微分方程的一般形式為：若c為常數(shù)，則該方程可以更簡(jiǎn)單地寫(xiě)成,7.6.2.3 雙曲型偏微分方程,3/19/2024星期六

64、, 2008-9- 6, 13:10:16,特征值型偏微分方程的一般形式為：若c為常數(shù)，則該方程可以更簡(jiǎn)單地寫(xiě)成,7.6.2.4 特征值型偏微分方程,3/19/2024星期六, 2008-9- 6, 13:10:16,7.6.3 偏微分方程的求解界面應(yīng)用舉例,偏微分方程求解程序概述偏微分方程求解區(qū)域繪制偏微分方程邊界條件描述偏微分方程求解舉例時(shí)變解的動(dòng)畫(huà)顯示函數(shù)參數(shù)的偏微分方程求解,3/19/2024星期六

65、, 2008-9- 6, 13:10:16,7.6.3.1 偏微分方程求解程序概述,偏微分方程工具箱——概述啟動(dòng)偏微分方程工具箱求解界面GUI在MATLAB下鍵入pdetool該界面分為四個(gè)部分菜單系統(tǒng)工具欄集合編輯求解區(qū)域,3/19/2024星期六, 2008-9- 6, 13:10:16,7.6.3.2 偏微分方程求解區(qū)域繪制,繪制并且定義偏微分方程求解區(qū)域繪制好求解區(qū)域后，單擊工具欄中 按

66、鈕就可以定義求解區(qū)域單擊D按鈕將求解區(qū)域用三角形劃分，單擊右側(cè)的按鈕加密網(wǎng)格；網(wǎng)格越密，計(jì)算的結(jié)果越精確，但計(jì)算時(shí)間則越長(zhǎng),3/19/2024星期六, 2008-9- 6, 13:10:16,Dirichlet條件：其中， 表示求解區(qū)域的邊界Neumann條件：其中， 為 向量法向的偏導(dǎo)數(shù),7.6.3.3 偏微分方程邊界條件描述,3/19/2024星期六, 20

67、08-9- 6, 13:10:16,7.6.3.4 偏微分方程求解舉例,設(shè)置了求解區(qū)域和邊界條件，并選擇了合適的偏微分方程后，就可以單擊工具欄的等號(hào)按鈕 = 立即得出微分方程的解微分方程的結(jié)果可以用其他很多方式顯示，如繪制等值線圖，引力線圖，網(wǎng)格型三維圖形,3/19/2024星期六, 2008-9- 6, 13:10:16,例 7.34,試求解雙曲型偏微分方程：,3/19/2024星期六, 2008-9- 6, 13:1

68、0: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à)顯示,3/19/2024星期六, 2008-9- 6, 13:10:16,7.6.3.6 函數(shù)參數(shù)的偏微分方程求解,偏微分方程工具箱能處理含有非線性系數(shù)的橢圓偏微分方程問(wèn)題使用x，y表示微分方程中 或