MATLAB-ch014數(shù)值計(jì)算-常微分方程2ppt課件

1、第第1414講講 數(shù)值計(jì)算數(shù)值計(jì)算常微分方程常微分方程張建瓴張建瓴內(nèi)容提要內(nèi)容提要 在科學(xué)研究和工程教學(xué)中會(huì)經(jīng)常遇到常微分方程。在科學(xué)研究和工程教學(xué)中會(huì)經(jīng)常遇到常微分方程。只含有一個(gè)自變量的微分方程稱為常微分方程只含有一個(gè)自變量的微分方程稱為常微分方程(Ordinary Differential Equations，ODE)。它或者沒(méi)有解析解，或。它或者沒(méi)有解析解，或者求取解析解的代價(jià)無(wú)法忍受，或者只有數(shù)值解等。者求取解析解的代價(jià)無(wú)法忍受，或者只有數(shù)值解等。在在MATLAB中，時(shí)常微分方程的解法一般有兩種：數(shù)值中，時(shí)常微分方程的解法一般有兩種：數(shù)值解和符號(hào)解解和符號(hào)解(解析解解析解)。常微分方

2、程的求解問(wèn)題可分為：初值問(wèn)題常微分方程的求解問(wèn)題可分為：初值問(wèn)題(Initial Value Problem，IVP)和邊值問(wèn)題和邊值問(wèn)題(Boundary Vaule Problem，BVP)。14.1 14.1 常微分方程的數(shù)值解常微分方程的數(shù)值解一般微分方程式描述系統(tǒng)內(nèi)部變量的變化率如何受系一般微分方程式描述系統(tǒng)內(nèi)部變量的變化率如何受系統(tǒng)內(nèi)部變量和外部激勵(lì)，如輸入的影響。統(tǒng)內(nèi)部變量和外部激勵(lì)，如輸入的影響。當(dāng)常微分方程式能夠解析求解時(shí)，可用當(dāng)常微分方程式能夠解析求解時(shí)，可用MATLAB的的符號(hào)工具箱中的功能找到精確解。符號(hào)工具箱中的功能找到精確解。在微分方程難以獲得解析解的情況下，可以方

3、便地在在微分方程難以獲得解析解的情況下，可以方便地在數(shù)值上求解。數(shù)值上求解。微分方程數(shù)值解是數(shù)值計(jì)算的基本內(nèi)容，由于微分方微分方程數(shù)值解是數(shù)值計(jì)算的基本內(nèi)容，由于微分方程的多樣性，有不同的解法。程的多樣性，有不同的解法。MATLAB給出的給出的7種解種解法分別由法分別由7個(gè)不同的函數(shù)來(lái)完成。個(gè)不同的函數(shù)來(lái)完成。 一、剛性問(wèn)題一、剛性問(wèn)題(Stiff Problem)1、剛性方程、剛性方程對(duì)于一個(gè)常微分方程組，如果其對(duì)于一個(gè)常微分方程組，如果其Jacobian矩陣的特征矩陣的特征值相差十分懸殊，那么這個(gè)方程組就稱為剛性方程組。值相差十分懸殊，那么這個(gè)方程組就稱為剛性方程組。對(duì)于剛性方程組，為保持

4、解法的穩(wěn)定，步長(zhǎng)選取根對(duì)于剛性方程組，為保持解法的穩(wěn)定，步長(zhǎng)選取根困難。有些解法不能用來(lái)解剛性方程組，而有些解困難。有些解法不能用來(lái)解剛性方程組，而有些解法對(duì)穩(wěn)定性的要求不嚴(yán)格，可以用來(lái)解決剛性問(wèn)題。法對(duì)穩(wěn)定性的要求不嚴(yán)格，可以用來(lái)解決剛性問(wèn)題。二、求解微分方程初值問(wèn)題的常用命令二、求解微分方程初值問(wèn)題的常用命令 MATLAB為解常微分為解常微分方程初值問(wèn)題提供的方程初值問(wèn)題提供的命令包括：微分方程命令包括：微分方程結(jié)算命令、被解算命結(jié)算命令、被解算命令、調(diào)用的常微分方令、調(diào)用的常微分方程文件格式命令、積程文件格式命令、積分算法參數(shù)選項(xiàng)分算法參數(shù)選項(xiàng)(Options)處理命令以處理命令以及輸出

5、處理命令等。及輸出處理命令等。這些命令的含義特點(diǎn)這些命令的含義特點(diǎn)和使用范圍如表所示。和使用范圍如表所示。 求解微分方程初值問(wèn)題的常用命令求解微分方程初值問(wèn)題的常用命令 （1表中單步法表中單步法是指只需要前一步是指只需要前一步的解即可計(jì)算出當(dāng)?shù)慕饧纯捎?jì)算出當(dāng)前的解，不需要附前的解，不需要附加初始值。在計(jì)算加初始值。在計(jì)算過(guò)程中可隨便改變過(guò)程中可隨便改變步長(zhǎng)而不會(huì)增加任步長(zhǎng)而不會(huì)增加任何附加的計(jì)算量。何附加的計(jì)算量。（2表中多步法表中多步法是指需要前幾次是指需要前幾次的解來(lái)計(jì)算當(dāng)前的解來(lái)計(jì)算當(dāng)前的解。的解。三、求解常微分方程的一般步驟三、求解常微分方程的一般步驟對(duì)于對(duì)于n階微分方程初值問(wèn)履，由于

6、函數(shù)及其直至階微分方程初值問(wèn)履，由于函數(shù)及其直至(n-1)階導(dǎo)階導(dǎo)數(shù)在某自變量點(diǎn)的值已知，所以由泰勒級(jí)數(shù)展開(kāi)，可算數(shù)在某自變量點(diǎn)的值已知，所以由泰勒級(jí)數(shù)展開(kāi)，可算出新的函數(shù)及導(dǎo)數(shù)值。出新的函數(shù)及導(dǎo)數(shù)值。在在MATLAB中，具體利用其命令來(lái)解初值問(wèn)題。常微分中，具體利用其命令來(lái)解初值問(wèn)題。常微分方程的步驟是：方程的步驟是：（1根據(jù)在工程實(shí)際中各學(xué)科的規(guī)律、定理和公式列出微根據(jù)在工程實(shí)際中各學(xué)科的規(guī)律、定理和公式列出微分方程和相應(yīng)的初始條件；分方程和相應(yīng)的初始條件；（2運(yùn)用變量替換，把一個(gè)高階方程寫成一階微分方運(yùn)用變量替換，把一個(gè)高階方程寫成一階微分方程組，初始條件也要做相應(yīng)地替換；程組，初始條件

7、也要做相應(yīng)地替換；求解常微分方程的一般步驟求解常微分方程的一般步驟（3根據(jù)變換后的一階微分方程組，編寫計(jì)算導(dǎo)數(shù)的根據(jù)變換后的一階微分方程組，編寫計(jì)算導(dǎo)數(shù)的M文件文件(在在MATLAB中稱為中稱為ODE文件文件)；（4使編寫好的使編寫好的ODE函數(shù)文件和變換后的初值供函數(shù)文件和變換后的初值供Y0微微分方程解算命令調(diào)用，運(yùn)行后即可得到分方程解算命令調(diào)用，運(yùn)行后即可得到Y(jié)(包含包含y及其導(dǎo)及其導(dǎo)數(shù)數(shù))在指定的時(shí)間區(qū)恥的數(shù)值解。在指定的時(shí)間區(qū)恥的數(shù)值解。四、求解常微分方程的命令格式四、求解常微分方程的命令格式MATLAB提供了多種解算常微分方程的命令提供了多種解算常微分方程的命令(如表如表13-1中所

8、列中所列)，這些命令函數(shù)的用法完全相同，以函數(shù)，這些命令函數(shù)的用法完全相同，以函數(shù)ode45為例，其最完整的調(diào)用格式為：為例，其最完整的調(diào)用格式為： T,Y=ode45(f,tspan,y0,options,p1,p2,) 式中：式中：f f是定義常微分方程的文件名字符串，這些函數(shù)是針是定義常微分方程的文件名字符串，這些函數(shù)是針對(duì)一階常微分方程組設(shè)計(jì)的，對(duì)于高階微分方程，必須先對(duì)一階常微分方程組設(shè)計(jì)的，對(duì)于高階微分方程，必須先化為形如化為形如 的一階微分方程組。在這里要注意：的一階微分方程組。在這里要注意：定義常微分方程組的函數(shù)必須以定義常微分方程組的函數(shù)必須以t t、x x為輸入?yún)?shù)為輸入?yún)?/p>

9、數(shù)( (順序不能順序不能變變) )，以，以 為輸出參數(shù)，并且輸出為輸出參數(shù)，并且輸出 必須是列向量。必須是列向量。 ),(txfx x x 參數(shù)說(shuō)明參數(shù)說(shuō)明tspantspan是一個(gè)向量，當(dāng)被賦予二元向量是一個(gè)向量，當(dāng)被賦予二元向量t0 tfinalt0 tfinal時(shí)，時(shí)，指定微分方程指定微分方程y=F(t,y)y=F(t,y)在初始條件在初始條件y0y0下，從下，從t0t0到到tfinaltfinal進(jìn)行積分；進(jìn)行積分；當(dāng)當(dāng)tspantspan被賦予多元向量被賦予多元向量t0,t1,tfinalt0,t1,tfinal時(shí)，在時(shí)，在tspantspan指定的時(shí)刻序列上求數(shù)值解。此時(shí)指定的時(shí)刻

10、序列上求數(shù)值解。此時(shí)tspantspan中的元素中的元素必須按單調(diào)次序排列。必須按單調(diào)次序排列。 y0是初始狀態(tài)列向量。是初始狀態(tài)列向量。 )1(0000nyyyy參數(shù)說(shuō)明參數(shù)說(shuō)明 options是一些可選的算法綜合參數(shù)，由函數(shù)是一些可選的算法綜合參數(shù)，由函數(shù)odeset進(jìn)行設(shè)置。進(jìn)行設(shè)置。當(dāng)輸入?yún)?shù)只有三個(gè)時(shí)，或第四個(gè)輸入?yún)?shù)位置上用空矩陣當(dāng)輸入?yún)?shù)只有三個(gè)時(shí)，或第四個(gè)輸入?yún)?shù)位置上用空矩陣“ ”賦值時(shí)，算法將使用賦值時(shí)，算法將使用options的缺省設(shè)置。的缺省設(shè)置。p1和和p2是傳遞給是傳遞給F的參數(shù)，從第的參數(shù)，從第5個(gè)輸入?yún)?shù)起，以后的各輸個(gè)輸入?yún)?shù)起，以后的各輸入?yún)?shù)都是向入?yún)?shù)都

11、是向ODE文件傳遞的參數(shù)。它們的使用可避免采用文件傳遞的參數(shù)。它們的使用可避免采用參數(shù)的全局變量方式傳遞。參數(shù)的全局變量方式傳遞。T是輸出參數(shù)，為所求數(shù)值解的時(shí)間列向量，是輸出參數(shù)，為所求數(shù)值解的時(shí)間列向量，Y是數(shù)值解的矩是數(shù)值解的矩陣，它的每一行對(duì)應(yīng)著列向量陣，它的每一行對(duì)應(yīng)著列向量T中的一個(gè)時(shí)間值。中的一個(gè)時(shí)間值。例例14-1example14_1.m考慮描述振蕩器的經(jīng)典的范得波考慮描述振蕩器的經(jīng)典的范得波Var der Pol微分方微分方程的初值問(wèn)題。程的初值問(wèn)題。已知初始條件為已知初始條件為x0=1，x0=0。0)1 (222xdtdxxdtxd例例14-2 example14_2.m

12、解微分方程組解微分方程組(剛性問(wèn)題剛性問(wèn)題)：已知二階微分方程：已知二階微分方程y=1000(1-y2)y-y初始條件：初始條件：y(0)=0，y(0)=1，求該微分方程的解。，求該微分方程的解。 14.2 ODE14.2 ODE文件文件ODE文件是專門被微分方程解算命令文件是專門被微分方程解算命令(如如ode45)調(diào)用的調(diào)用的M函數(shù)文件。函數(shù)文件。ODE文件的功能很多，其中計(jì)算每一步的導(dǎo)文件的功能很多，其中計(jì)算每一步的導(dǎo)數(shù)是其最基本、必不可少的功能。另外，還有一些常用的數(shù)是其最基本、必不可少的功能。另外，還有一些常用的功能，如解算指令中功能，如解算指令中options的參數(shù)設(shè)置、事件設(shè)置等。

13、的參數(shù)設(shè)置、事件設(shè)置等。MATLAB給給ODE文件制定了比較嚴(yán)整的格式和規(guī)定。而文件制定了比較嚴(yán)整的格式和規(guī)定。而對(duì)于初學(xué)者來(lái)說(shuō)，這樣的對(duì)于初學(xué)者來(lái)說(shuō)，這樣的ODE文件就顯得比較復(fù)雜。為文件就顯得比較復(fù)雜。為幫助用戶理解和編寫幫助用戶理解和編寫ODE文件，文件，MATLAB提供了一個(gè)提供了一個(gè)“模板文件模板文件odefile.m。用戶在自己編寫。用戶在自己編寫ODE文件時(shí)，可文件時(shí)，可以借助此模板，在模板中填入適當(dāng)?shù)膬?nèi)容，刪去不需要以借助此模板，在模板中填入適當(dāng)?shù)膬?nèi)容，刪去不需要的內(nèi)容。的內(nèi)容。一、一、ODEODE文件模板文件模板一個(gè)一個(gè)ODEODE文件就是被微分方程解算指令文件就是被微分方程

14、解算指令( (如如od45)od45)調(diào)用的調(diào)用的M M函數(shù)文件，用戶可根據(jù)需要，使用任何其他文件名。函數(shù)文件，用戶可根據(jù)需要，使用任何其他文件名。 在缺省情況下，在缺省情況下，ODEODE解算指令處理形如解算指令處理形如dy/dt=F(t,y)dy/dt=F(t,y)的微分的微分方程問(wèn)題，方程問(wèn)題，t t是自變量，是自變量，y y是一階微分方程組的列向量形式是一階微分方程組的列向量形式函數(shù)。解算時(shí)，解算指令重復(fù)調(diào)用函數(shù)。解算時(shí)，解算指令重復(fù)調(diào)用F=odefile(T,Y)F=odefile(T,Y)，此處，此處T T是一個(gè)時(shí)間序列，是一個(gè)時(shí)間序列，Y Y是一個(gè)列向量，輸出參數(shù)是一個(gè)列向量，輸

15、出參數(shù)F F是一個(gè)同長(zhǎng)是一個(gè)同長(zhǎng)度的列向量度的列向量ODEODE文件的經(jīng)典最簡(jiǎn)形式為：文件的經(jīng)典最簡(jiǎn)形式為：function F=odefile(t,y)function F=odefile(t,y)F=;F=;ODE文件模板格式文件模板格式在在ODEODE文件中，最完整的格式為：文件中，最完整的格式為：function varargout=odefile(t,y,flag,p1,p2,)function varargout=odefile(t,y,flag,p1,p2,)其中：其中：t t和和y y是積分變量，是積分變量，t t必須作為第一輸入?yún)?shù)，而必須作為第一輸入?yún)?shù)，而y y必須作為第

16、二輸入?yún)?shù)。必須作為第二輸入?yún)?shù)。plpl和和p2p2是被傳遞的參數(shù)，這里是被傳遞的參數(shù)，這里僅列出兩個(gè)參數(shù)，實(shí)際上，傳遞參數(shù)的數(shù)目不受限制。僅列出兩個(gè)參數(shù)，實(shí)際上，傳遞參數(shù)的數(shù)目不受限制。flagflag是切換變量，它必須處在第三個(gè)輸入?yún)?shù)位置上是切換變量，它必須處在第三個(gè)輸入?yún)?shù)位置上二、二、ODE模板的使用方法模板的使用方法1 1、ODEODE模板使用步驟模板使用步驟如果用戶需要編寫比較復(fù)雜的如果用戶需要編寫比較復(fù)雜的ODE文件，該模板將十分有文件，該模板將十分有用。其具體步驟如下：用。其具體步驟如下：（1在在MATLAB命令窗口中運(yùn)行命令窗口中運(yùn)行help odefile，就會(huì)出現(xiàn)，就

17、會(huì)出現(xiàn)ODE模板；模板；（2將此模板文件全部?jī)?nèi)容復(fù)制到文件編輯器中；將此模板文件全部?jī)?nèi)容復(fù)制到文件編輯器中；（3用戶根據(jù)需要，將模板文件中的那些不被使用的情況用戶根據(jù)需要，將模板文件中的那些不被使用的情況(Case)連同相應(yīng)的子函數(shù)刪除；連同相應(yīng)的子函數(shù)刪除；如無(wú)特別原因，盡量少改動(dòng)模板中所有的子函數(shù)名、輸如無(wú)特別原因，盡量少改動(dòng)模板中所有的子函數(shù)名、輸入和輸出參數(shù)名。用戶根據(jù)具體問(wèn)題編寫相應(yīng)的程序。入和輸出參數(shù)名。用戶根據(jù)具體問(wèn)題編寫相應(yīng)的程序。1 1、ODEODE模板使用步驟模板使用步驟（4 4把編寫好的文件存入用戶自己的工作目錄，以備調(diào)用。把編寫好的文件存入用戶自己的工作目錄，以備調(diào)用。

18、編寫好的編寫好的ODEODE文件用戶可以給一個(gè)新名字，但要注意文件用戶可以給一個(gè)新名字，但要注意以下三者間的一致：存放在工作目錄上的名稱；函數(shù)以下三者間的一致：存放在工作目錄上的名稱；函數(shù)文件第一行中等式右邊的函數(shù)名稱；在微分方程解算文件第一行中等式右邊的函數(shù)名稱；在微分方程解算指令調(diào)用該指令調(diào)用該ODEODE文件時(shí)的字符串名稱。文件時(shí)的字符串名稱。2 2、MATLABMATLAB提供的微分方程幫助資源提供的微分方程幫助資源MATLABMATLAB為微分方程指令的使用提供了豐富的幫助資源。它們?yōu)槲⒎址匠讨噶畹氖褂锰峁┝素S富的幫助資源。它們的分布形式和位置為：的分布形式和位置為：（1 1PDEP

19、DE文件文件matlabRl2helppdf_docmatlabusing_ml.pdfmatlabRl2helppdf_docmatlabusing_ml.pdfmatlabR12helppdf_docmatlabrefrefbook.pdfmatlabR12helppdf_docmatlabrefrefbook.pdf（2 2在線幫助在線幫助使用使用help funfunhelp funfun以及對(duì)具體命令進(jìn)一步使用以及對(duì)具體命令進(jìn)一步使用helphelp，或使，或使用用helpwinhelpwin查閱。查閱。M A T L A BM A T L A B 實(shí) 現(xiàn)實(shí) 現(xiàn) O D EO D E

20、 演 示 的演 示 的 M M 文 件 位 于 ：文 件 位 于 ：matlabRl2toolboxmatlabdemos matlabRl2toolboxmatlabdemos 14.3 常微分方程的符號(hào)解常微分方程的符號(hào)解從數(shù)值計(jì)算角度看，與初值問(wèn)題求解相比，微分方程從數(shù)值計(jì)算角度看，與初值問(wèn)題求解相比，微分方程邊值問(wèn)題的求解顯得更為復(fù)雜和困難。而這一問(wèn)題對(duì)邊值問(wèn)題的求解顯得更為復(fù)雜和困難。而這一問(wèn)題對(duì)于符號(hào)計(jì)算來(lái)說(shuō)，不論是初值問(wèn)題，還是邊值問(wèn)題，于符號(hào)計(jì)算來(lái)說(shuō)，不論是初值問(wèn)題，還是邊值問(wèn)題，其求解微分方程的命令形式相同，并且相當(dāng)簡(jiǎn)單。其求解微分方程的命令形式相同，并且相當(dāng)簡(jiǎn)單。但符號(hào)計(jì)算

21、可能花費(fèi)較多的時(shí)間，也可能得不到簡(jiǎn)單的但符號(hào)計(jì)算可能花費(fèi)較多的時(shí)間，也可能得不到簡(jiǎn)單的解析解以及得不到封閉式的解，甚至也可能無(wú)法求解。解析解以及得不到封閉式的解，甚至也可能無(wú)法求解。不論怎樣，微分方程的符號(hào)解法和數(shù)值解法有很好的互不論怎樣，微分方程的符號(hào)解法和數(shù)值解法有很好的互補(bǔ)作用。補(bǔ)作用。一、微分方程符號(hào)解的一般命令一、微分方程符號(hào)解的一般命令在在MATLABMATLAB中，函數(shù)中，函數(shù)dsolvedsolve可用來(lái)求常微分方程的符號(hào)解。在符號(hào)方程中，可用來(lái)求常微分方程的符號(hào)解。在符號(hào)方程中，用符號(hào)表達(dá)式中包含的字母用符號(hào)表達(dá)式中包含的字母“D“D來(lái)代替微分運(yùn)算，符號(hào)來(lái)代替微分運(yùn)算，符號(hào)D

22、2D2，D3D3，DNDN分別對(duì)應(yīng)于第二、第三、分別對(duì)應(yīng)于第二、第三、第、第N N導(dǎo)數(shù)。這樣，導(dǎo)數(shù)。這樣，D2yD2y就是就是d2y/dt2d2y/dt2在符號(hào)數(shù)在符號(hào)數(shù)學(xué)工具箱中對(duì)應(yīng)的形式。因變量是位于學(xué)工具箱中對(duì)應(yīng)的形式。因變量是位于D D后面的變量，缺省的自變量為后面的變量，缺省的自變量為t t。還可以在命令中加入包含自變量的參數(shù)，將自變量變?yōu)槠渌?hào)變量。還可以在命令中加入包含自變量的參數(shù)，將自變量變?yōu)槠渌?hào)變量。函數(shù)函數(shù)dsolvedsolve的調(diào)用格式為：的調(diào)用格式為：S=dsolve(eqnl,eqn2,)1、輸入?yún)?shù)的內(nèi)容、輸入?yún)?shù)的內(nèi)容輸入?yún)?shù)包含三部分內(nèi)容：微分方程、初始

23、條件和指定獨(dú)立輸入?yún)?shù)包含三部分內(nèi)容：微分方程、初始條件和指定獨(dú)立變量。其中，微分方程是必不可少的輸入內(nèi)容，其余內(nèi)容視變量。其中，微分方程是必不可少的輸入內(nèi)容，其余內(nèi)容視需要而定，可有可無(wú)。輸入的參數(shù)必須以字符形式編寫。需要而定，可有可無(wú)。輸入的參數(shù)必須以字符形式編寫。2 2、獨(dú)立變量的規(guī)定、獨(dú)立變量的規(guī)定關(guān)于指定獨(dú)立變量的規(guī)定：若要指定獨(dú)立變量，則總是關(guān)于指定獨(dú)立變量的規(guī)定：若要指定獨(dú)立變量，則總是由全部輸入?yún)?shù)由全部輸入?yún)?shù)eqnleqnl、eqn2eqn2、中的量后一個(gè)參數(shù)中的量后一個(gè)參數(shù)eqnneqnn定義。若不對(duì)獨(dú)立變量加以專門定義，則該指令定義。若不對(duì)獨(dú)立變量加以專門定義，則該指令將默認(rèn)小寫英文字母將默認(rèn)小寫