常微分方程數(shù)值解法簡介.ppt

1、科學計算與數(shù)學建模,中南大學數(shù)學科學與計算技術(shù)學院,第7章 常微分方程數(shù)值解法簡介,第七章常微分方程數(shù)值解法簡介,微分方程在科學和工程技術(shù)中有很廣泛的應(yīng)用。許多實際問題的數(shù)學模型都可以用微分方程來描述，歸結(jié)為常微分方程的定解問題。很多偏微分方程問題，也可以化為常微分方程問題來近似求解，但是求出所需的解絕非易事。實際上，除了極特殊情形外，人們不可能求出微分方程的解析解，只能用各種近似方法得到滿足一定精度的近似解。 在常微分方程中已經(jīng)熟悉了級數(shù)解法和Picard逐步逼近法，這些方法可以給出解的近似表達式，稱為近似解析方法。另一類方法只給出解在一些離散點上的值，稱為數(shù)值方法。數(shù)值方法應(yīng)用范圍更廣，特

2、別適合用計算機計算，本章主要介紹常用的常微分方程數(shù)值解法,函數(shù)是事物的內(nèi)部聯(lián)系在數(shù)量方面的反映,如何尋找變量之間的函數(shù)關(guān)系,在實際應(yīng)用中具有重要意義.在許多實際問題中，往往不能直接找出變量之間的函數(shù)關(guān)系,但是有時卻容易找出變量的改變量之間的關(guān)系，從而建立描述問題的微分方程模型,1 實際問題的微分方程模型,例7.1 將初始溫度 的一碗湯放置于環(huán)境溫度 保持在 的桌上，10 分鐘后測得湯的溫度為1000C。如果湯的溫度低于550C才可以喝，問再過20分鐘后這碗湯能喝了嗎,設(shè)物體在 時刻的溫度為 ，從 ，溫度從 ，注意到熱量總是從溫度高的物體向溫度低的物體傳導，因而 ，所以溫度差 恒正，又因物體將隨

3、時間而逐漸冷卻，則溫度的改變量為 兩邊除以 ，并令 得溫度變化速度為,解： 為了解決這一問題，需要了解有關(guān)熱力學的一些基本規(guī)律.熱量總是從溫度高的物體向溫度低的物體傳導的；在一定的溫度范圍內(nèi)，一個物體的溫度變化速度與這個物體的溫度和其所在介質(zhì)溫度的差值成正比,其中 是比例常數(shù).從而得出描述物體冷卻過程的微分方程模型為 容易求出這個一階微分方程初值問題的解為 根據(jù)問題所給的條件知 ,當時， ，得到 將 ， 代入，得,7.1.1,7.1.2,從而得到這碗湯的溫度隨時間變化的函數(shù)關(guān)系為 于是，將 代入計算得到再過20 min湯的溫度 ，這說明再過20 min后這碗湯能喝了. 不過，并不是所有的微分方

4、程模型都可求出解析解。例如，看似簡單的微分方程 ，自德國數(shù)學家Wilhelmvon Leibniz提出100多年后才被法國數(shù)學家Joseph Liouville證明它沒有解析解，只能借助于數(shù)值的方法求數(shù)值解,7.1.3,例7.2 某地區(qū)發(fā)現(xiàn)一種有免疫性的傳染病，為了控制疫情擴散對該地人 群進行隔離處理.為了分析受感染人數(shù)的變化規(guī)律，需要建立描述傳染 病傳播過程的數(shù)學模型. 解 設(shè)該地區(qū)的總?cè)藬?shù)為常數(shù)，任意時刻病人、健康人和病人治愈后 移出感染系統(tǒng)的移出者的比例分別為 ，病人的日接觸率 ， 日治愈率 ，則容易得出從 時刻，病人和健康人的改變量為,每個方程兩邊除以 ，并令 ，化簡后得,7.1.4,

5、其中 （對任意的 t,式（7.1.4）就是描述病人和健康人的比例 和 隨時間變化的微分方程模型，這是一個微分方程組的初值問題.但是，這一初值問題的解析解是無法求出的，因此不能直接利用 和 的解析式來分析和解決問題。,在數(shù)學建模課程中學到的大量數(shù)學模型都是用微分方程形式給出的，各類微分方程本身和它們的解所具有的特性在常微分方程及數(shù)學物理方程中有所解釋.雖然，求解微分方程有許多解析方法，但解析方法只能夠求解一些特殊類型的方程，在實際應(yīng)用中人們更關(guān)心的是某些特定的自變量在某一個定義范圍內(nèi)的一系列離散點上的近似值.這樣一組近似解稱為微分方程在該范圍內(nèi)的數(shù)值解，尋找微分方程數(shù)值解的過程稱為微分方程的數(shù)值

6、解法,2 簡單的數(shù)值方法與基本概念,設(shè) 在區(qū)域 上連續(xù)，求 滿足 其中 是已知常數(shù),這就是一階常微分方程的初值問題. 為使問題（7.2.1）的解存在、唯一且連續(xù)依賴初值 ，即初值問題（7.2.1）適定，還必須對右端項 加以適當限制，通常要求 關(guān)于 是已知函數(shù),且滿足Lipschitz條件，即存在常數(shù)L，使,7.2.1 常微分方程初值問題,7.2.1,7.2.2,對所有 及 成立。 本章總假定滿足條件（7.2.2）。 1. Euler方法的導出與幾何意義 最簡單的數(shù)值解法是Euler法。 將區(qū)間 作N等分，小區(qū)間的長度 稱為步長，點列 稱為節(jié)點， 。 由已知初值 ，可算出 在 的導數(shù),7.2.2

7、 Euler法及改進的Euler法,其中 ，并略去二階小量 ，得,下面用3種方法導出Euler法.本章用 表示函數(shù) 在 點 的精確值， 表示 的近似值,就是 的近似值。利用 可算出 ，如此下去可算 出 在所有節(jié)點上的值，它的一般遞推公式為,7.2.3,1） 冪級數(shù)展開法 利用Taylor展式,7.2.4,這就是Euler法,實際上，初值問題（7.2.1）的解是xy平面上過點 的一條積分曲線。按Euler法，過初始點 作經(jīng)過此點的積分曲線的切線（斜率為 ），沿切線取點 （ 按式（7.2.4）計算）作為 的近似.然后，過 作經(jīng)過此點的積分曲線的切線（斜率為 ），沿切線取點 （ 按式（7.2.4）計

8、算）作為 的近似.如此下去，即得一條以 為頂點的折線，這就是用Euler法得到的近似積分曲線，如圖7-1所示.從幾何圖形上看， 越小，此折線逼近積分曲線越好，因此也稱Euler法為Euler折線法,Euler法有明顯的幾何意義,2 ) 數(shù)值微分法,利用向前差商近似導數(shù),從而得出Euler法的一般遞推公式,3 ) 數(shù)值積分法,將初值問題（7.2.1）寫成等價的積分形式,取 ，得,用左矩形公式作為右端積分的近似，并用 替代 ，即得,從而也可得出Euler法的一般遞推公式為,7.2.5,2. 改進的Euler方法,由Euler方法的數(shù)值積分導出法可知只要給出式（7.2.5）右端定積分的一種近似計算方

9、法，就可得出初值問題（7.2.1）的一種數(shù)值求解方法,如果用右矩形公式作為式（7.2.5）右端積分的近似，則可得,從而也可得出一般遞推公式為,稱式（7.2.6）為后退Euler法,7.2.6,顯然，改進的Euler法比Euler法精度更高。后退Euler法和改進的Euler法，由于未知數(shù) 同時出現(xiàn)在等式的兩邊，故稱為隱式算法；如果未知數(shù) 由已知量直接計算（即不出現(xiàn)在等式右端），則稱為顯式算法。對于隱式算法，每步計算需要解關(guān)于 的方程，而這樣的方程往往是非線性的，通常將初值取為 ，用迭代法求解，一般只需迭代幾步即可收斂。一般先用顯式公式計算一個初值，再用隱式公式迭代求解,如果用梯形公式近似替代式

10、（7.2.5）右端的定積分，則可得,從而得出一般遞推公式為,稱式（7.2.7）為改進的Euler法,7.2.7,如果先用顯式Euler公式作預測，算出 再將 代入隱式梯形公式的右邊作校正，得到,從而可得,這種方法稱為預估校正法?？梢钥吹剿秋@示格式，迭代求解過程比隱式公式的簡單。后面將看到，它的穩(wěn)定性高于顯式Euler法,如果在區(qū)間,上對初值問題（7.2.1）的方程兩邊積分，則有,并用中矩形求積公式近似替代右端的定積分，則得出一般遞推公式為,稱式（7.2.8）為Euler中點公式,7.2.8,和,這樣的方法稱為雙步法（或,如果計算,的近似值,時只用到前一節(jié)點的值,則從初值,樣的方法稱為單步法；

11、而Euler中點公式計算,到前兩個節(jié)點的值,多步法.多步法附加初值才能逐一計算出以后各節(jié)點的值,出發(fā)可逐一計算出以后各節(jié)點的值，這,時需要用,二步法）；計算,時需要用到前面多個節(jié)點值的方法稱為,7.2.3 截斷誤差與算法精度的階,從Euler法的幾何意義得知，由Euler法所得的折線明顯偏離了積分曲線，可見此方法非常粗糙（即誤差太大）?，F(xiàn)在分析一下求解初值問題（7.2.1）的數(shù)值方法誤差的來源。為使問題簡化，不考慮因計算機字長限制引起的攝入誤差,在假設(shè)第 步計算是精確的前提下 （即 ），第 步計算 的截斷誤差 稱為局部截斷誤差。若某算法的局部截斷誤差為 ，即為 的同階無窮小，則稱該算法有 階精

12、度。若局部截斷誤差,則稱 為誤差主項， 為誤差主項系數(shù),1. Euler法的局部截斷誤差,由Euler法的一般遞推公式和 的Taylor展式，得,所以，Euler法的局部截斷誤差為 ，即Euler法為1階精度算法，其誤差主項為,2. 后退Euler法的局部截斷誤差,同理，由后退Euler法的一般遞推公式， 和 的Taylor 展式，得,所以，后退Euler法的局部截斷誤差也是 ，即后退Euler法也是1階精度算法，其誤差主項為,3. 改進Euler法的局部截斷誤差,由改進Euler法的一般遞推公式可得其局部截斷誤差為,所以，改進Euler法的局部截斷誤差為 ，即改進Euler法是2階精度算法，

13、其誤差主項為,4. 整體截斷誤差,當然人們更關(guān)心的是近似解 的誤差，即 ，稱為整體截斷誤差.由 的Taylor展式與Euler法的一般遞推公式相減，得,記 ，因 關(guān)于 滿足Lipschitz條件，所以存在常數(shù)L，使得,以此遞推，得,注意到 ， ，于是,上式右端依賴于初始誤差 和局部截斷誤差的上界R 。對于Euler法，可取 （C 是與n 無關(guān)的常數(shù)），若 ，即 ，則,所以 ，比局部截斷誤差低一階。用同樣的方法可以證明改進Euler法的整體截斷誤差的階為 ，也比局部截斷誤差低一階,5. Euler算法的穩(wěn)定性,在實際計算中，由于測量誤差、舍入誤差等因素的影響，初值 往往不能精確給出，其誤差將依次

14、傳遞下去。如果傳遞誤差能夠被控制，精確地說，傳遞誤差連續(xù)依賴于初始誤差，則稱算法穩(wěn)定；否則就說算法不穩(wěn)定。顯然，不穩(wěn)定的算法是不能用的。下面僅考察Euler法的穩(wěn)定性,設(shè)從初值 和 算出的節(jié)點值分別為 和 ，則滿足,兩式相減，并令 ，則,從而,因,這說明 連續(xù)依賴初始誤差 ，即Euler法穩(wěn)定,同樣可證改進的Euler法也穩(wěn)定,例7.3 取步長 ,分別用Euler法、后退Euler法和中點法 求解初值問題,解 步長 ，所以各節(jié)點,1) 因為 利用Euler法的計算公式,可得,2) 利用后退Euler法的計算公式,可解得的顯示表達,于是，可得,3) 由中點法的計算公式,可知需要兩個初值。在此，我

15、們利用后退歐拉法計算的結(jié)果 再依次計算得,而該初值問題的解析解是 用它計算各節(jié)點的函數(shù)值，可得,將上述3種方法計算的結(jié)果同精確值對照,可以看出它們確實都是精確值的近似值,只是誤差不一樣。Euler法的誤差較小,后退Euler法誤差偏大,中點法誤差最小,3 線性多步法,用Euler法計算節(jié)點 的近似值 時只用到前一節(jié)點的值 ，是線性的單步法.為了提高解的精度，需要構(gòu)造線性多步法，其一般形式為,其中 ， 和 是常數(shù)，且 ， 和 不同時為0 。按公式（7.3.1）計算 時，要用到前面 個節(jié)點的值 ，因此式（7.3.1）稱為多步法（或 步法）。又因為方程（7.3.1）關(guān)于 是線性的，所以稱為線性多步法

16、。顯然，若 ，則線性多步法（7.3.1）是顯式的；若 ，則線性多步法（7.3.1）是隱式的.用線性多步法進行計算 時，除需要給定 外，還要附加初值 ，這可以用其他方法計算。由于多步法每計算一步用到的信息更多，因此希望由此構(gòu)造出精度更高的算法,7.3.1 數(shù)值積分法,將微分方程 在 上積分，得,7.3.2,適當選取 個節(jié)點，作被積函數(shù) 的 次Lagrange插值多項式 并用 近似代替 ，就可得到形如式（7.3.1）的線性多步法。插值節(jié)點的不同取法就可得出不同的線性多步法,1. Adams外插法,取 為節(jié)點，構(gòu)造 的Lagrange插值多項式 ，則,其中是插值余項.將式（7.3.3）代入式（7.3

17、.2），得,舍去余項,并用 代替 ，即得,由此可見，Adams法的局部截斷誤差為,7.3.3,7.3.4,7.3.5,7.3.6,下面給出Adams法（7.3.6）的具體形式。假設(shè)前 個節(jié)點處的函數(shù)值已知，即 的近似值 已算出，從而函數(shù)值 也已知。這樣，就可以利用 個數(shù)據(jù)點 , ， 構(gòu)造被積函數(shù) 的 次插值多項式,其中 是Lagrange插值基函數(shù).從而可得,7.3.6,記 , 則有,這就是求解初值問題的Adams顯式公式。 是關(guān)于 和 的線性表達式，所以它是線性 步法,在上述Adams顯式公式的推導中，選用了 , 作為 插值節(jié)點，但構(gòu)造的 次插值多項式 是代替區(qū) 間 上的未知函數(shù) ，因此屬于

18、“外插”，稱為 Adams外插法，也稱為Adams-Bashorth法.顯然，當 時， Adams外插法就是Euler法,2. Adams內(nèi)插法,如果將Adams外插法推導過程中的節(jié)點改為 , 公式（7.3.7）就相應(yīng)地變?yōu)?由于式（7.3.8）右端含有 (可能是非線性表達式)，所以式（7.3.8）屬于隱式公式，稱為Adams隱式公式，且插值區(qū)間包含了積分區(qū)間 ，因此屬于“內(nèi)插”，稱為Adams內(nèi)插法，也稱為Adams-Moulton法。顯然，當 時，Adams內(nèi)插法就是改進的Euler法,表7-1和表7-2分別給出了當 時Adams顯式公式和隱式公式的系數(shù)值,7.3.8,表7-1Adams外

19、插法系數(shù)值,表7-2Adams內(nèi)插法系數(shù)值,利用插值多項式的余項可以求出Adams法的局部截斷誤差，對指定的 ，表7-3 列出了它們的局部截斷誤差主項,表7-3Adams法（步）的局部截斷誤差主項,應(yīng)該指出，用數(shù)值積分法只能構(gòu)造一類特殊的多步法，其系數(shù)滿足 ， ， （當,下面介紹更為一般的待定系數(shù)法,7.3.2待定系數(shù)法,為了分析一般線性多步法的局部截斷誤差，令,設(shè) 是初值問題的解，將 和 在點 用Taylor公式展開，代入式（7.3.9）按 的同次冪合并同類項，得,其中,若 有 次連續(xù)微商，則可選取足夠大的 和 ， ， 使 ,而,即選取 , 滿足,此時，有,令 ，則,7.3.10,于是，由滿

20、足線性方程組（7.3.10）的 ， 得到的線性多步法（7.3.1）的局部截斷誤差為 ， 可以證明此線性多步法的整體截斷誤差的階是 ，所以此線性多步法為 階 步法。顯然， 的大小和 有關(guān),因為線性多步法（7.3.1）可以相差一個非零乘數(shù)，所以不妨設(shè) 。當 時， 可用 直接表示，稱為顯式法;反之，當 時，求 需要解一個方程（一般用迭代法），稱為隱式法.用待定系數(shù)法構(gòu)造線性多步法的一個基本要求是選取 ， 使局部截斷誤差的階盡可能高,1. Milne方法,作為待定系數(shù)法的一個應(yīng)用，下面討論一般的2步 法。此時 ， ，其余5個系數(shù) 由 確定，即,其中 為任意常數(shù),所以，一般的2步法為,7.3.11,由,

21、得,這是4階2步法，是具有最高階的2步法，稱之為Milne方法,7.3.12,所以，當 時， ，此時的2步法（7.3.11）是 3階2步法；當 時， ，但 ，此時方法（7.3.11）化為,這一方法也可用Simpson公式導出.此外，若取 ，則2步法（7.3.11）為2步Adams內(nèi)插法；若取 ，則2步法（7.3.11）為顯式法,它的余項為,2. Hmaming方法,用待定系數(shù)法容易求出，4階3步的公式,這就是著名的4階Hmaming公式，它的余項為,更多常用的線性多步法和多步法計算中的問題等可以閱讀相關(guān)參考文獻，在此不再贅述.關(guān)于線性多步法的穩(wěn)定性、收斂性和誤差估計，以及絕對穩(wěn)定性和相對穩(wěn)定性

22、等基本理論問題也請參閱相關(guān)的參考文獻,7.3.13,4 非線性高階單步法Runge-Kutte法,Euler法是最簡單的單步法，單步法不需要附加初值，所需的存儲量小，改變步長靈活，但是線性單步法的階最多是2.本節(jié)介紹非線性（關(guān)于 ）高階單步法，主要介紹Runge-Kutta法,7.4.1 泰勒展開法,設(shè)初值問題（7.2.1）的解充分光滑，將在點用Taylor公式展開,其中 是介于 與 之間的常數(shù)，是未知函數(shù)在點的階導數(shù)，但它的值可以利用微分方程本身來計算,舍去式（7.4.1）中的含有的項，則得到求解初值問題的非線性單步法,7.4.3,其中 按式（7.4.2）來計算.這一方法的局部截斷誤差為 ，

23、因此，它是非 線性階單步法,由于需要計算 的高階偏導數(shù)，計算量太大，所以一般不用式 （7.4.3）作數(shù)值計算，但可用它計算附加初值,7.4.2 Runge-Kutta法,將初值問題在區(qū)間上寫成積分形式,同樣，利用Euler法又可以算出,其中,7.4.7,7.4.8,下面推導一些常用的計算方案.將,展開 到,的3次冪,于是，由R-K法得,1. 一級R-K法,可見一級R-K法就是Euler法,2. 二級R-K法,令,此時,于是,與,比較,的系數(shù)，得,它有無窮多組解，從而有無窮多個二級二階R-K算法.兩個常見的方法分別如下,稱此為中點法，這是一種修正的Euler法,3. 三級R-K法,4個方程不能完

24、全確定2個系數(shù)，含有2個自由參數(shù)，因此有無窮多個三級三階R-K算法,僅列舉兩個常見的方法,7.4.9,稱此為Heun三階方法,7.4.10,4. 四級R-K法,7.4.11,和,7.4.12,這是最常用的四階R-K法，稱為標準（或經(jīng)典）R-K法,但是，它需要計算,的次數(shù)比預估-校正法多.這里介紹的R-K法都是,仿此，還可以構(gòu)造出更多的算法，容易,階R-K法整體截斷誤差,R-K法和預估-校正法是求解常微分方程初值問題,證明的階為,的兩大類有效數(shù)值方法.R-K法的絕對穩(wěn)定域一般比預估-校正法的大,顯式的，為了進一步改善穩(wěn)定域，可采用隱式R-K法（請參考相關(guān) 文獻,解 為了比較3種方法的計算精度，將

25、Euler法的步長取為0.025，改進的Euler法的步長取為0.05，R-K方法的步長取為0.1，則3種方法的變量每增加0.1時，都需要計算4個函數(shù)值.現(xiàn)將計算結(jié)果列于表7-4中,從計算結(jié)果可以看出,標準R-K方法比另外兩種方法的精度好很多。在,處，3種方法的誤差分別是,表7-4 例7.4的3種計算結(jié)果比較,解,利用標準的R-K公式,依次計算： (1) 當,時,接下來，利用修正的Milne-Hamming組成的預估校正法公式,依次計算,5 一階方程組和高階方程的初值問題,前面研究了單個方程 初值問題的數(shù)值解法，只要把 和 理解為向量，那么，所提供的各種計算公式即可應(yīng)用到一階方程組的情形,含個

26、方程的一階方程組初值問題的一般形式為,7.5.1,如果實際問題不是一階方程組而是高階方程，也可以把它化為一階方程組,例如，m階微分方程,7.5.2,只要引進新變量組,式（7.5.2）就可化為一階方程組,7.5.3,這種轉(zhuǎn)化不僅是理論上的需要，在計算上也可能更為方便,引進向量記號,則式（7.5.1）可寫為向量形式,7.5.4,若關(guān)于滿足Lipschitz條件，則初值問題（7.5.1）有唯一解,7.5.5,其中 是標量， 是向量值算子,例如，解方程組初值問題的線性多步法為,前面介紹的線性多法、預估-校正法和R-K法都可以直接推廣到一階方程組，這只需用向量代替相應(yīng)的標量。所有關(guān)于階、相容性、穩(wěn)定性和

27、收斂性的定義和結(jié)論都可推廣到方程組，而將絕對值符號 換成m維歐氏空間向量的模,6 常微分方程邊值問題的數(shù)值解法,常微分方程邊值問題的一般形式為：求函數(shù),使之滿足,7.6.1,常微分方程邊值問題的基本數(shù)值解法分為兩類，一類是將它轉(zhuǎn)化成初值問題來求解，第二類方法是利用數(shù)值微商的方法將它轉(zhuǎn)化成線性或非線性方程組求解,7.6.1試射法,將邊值問題轉(zhuǎn)化成如下形式初值問題,7.6.2,即依據(jù)邊值條件尋求與它等價的初始條件,然后，令,從而使問題（7.6.2）轉(zhuǎn)化為一階微分方程組,最后，令,其中,這樣方程組（7.6.3）就具有如下形式,7.6.4,顯然，方程組（7.6.4）與標準的常微分方程初值問題（7.2.

28、1）具有相同的形式，因而，所有阿求解初值問題（7.2.1）的方法都可以用來求解方程組（7.6.4）。所不同的是，只需要把原公式中等分別換成向量。下面以Euler法為例說明試射法的基本方法，并用例子說明標準R-K方法的試射法求解過程,或用分量形式表示成,因此，利用試射法求解邊值問題的關(guān)鍵是如何把邊值條件轉(zhuǎn)化為等價的初值條件，即確定 。具體方法如下,1)憑經(jīng)驗提供 的兩個預測值,并按這兩個斜率值“試射,所謂試射，就是按上述試射法的基本步驟分別求解對應(yīng)的初值問 題。假設(shè)它們的解為,計算,以得到,的兩個近似值,2)如果,均不滿足預定的精度，就用線性插值方法校正,即選擇新的斜率值,將Euler法寫成相應(yīng)

29、的向量形式為,即利用 再選擇新的斜率值。 重復上述計算過程，直到找到合適的斜率值的近似值，顯然在該初值條件下得到的初值問題的解也是原問題（7.6.1）的解。,3) 用 試射，又會得到對應(yīng)的初值問題的解,和,的近似值,如果,不滿足預定的精度，回到計算過程(2,和,解 假設(shè),則對應(yīng)一階微分方程組初值問題為,例7.6以標準R-K方法為基礎(chǔ)，用試射法求解如下問題,選取,的近似值,用標準R-K方法求解上述方,再選取,的近似值,用標準R-K方法求解上述方程組，求得,由,作線性插值，計算,的新的近似值,程組，求得,并由此得,由,作線性插值，計算,的新的近似值,并由此得,由,作線性插值，計算,的新的近似值,計算得,算法終止。這時所得到的原問題邊值問題的數(shù)值解如下,1. 差分法的計算步驟 為求解邊值問題（7.6.1），差分法的思想是用向前差商,代替,或向后差商,或中