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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

16、。顯然，若 ，則線性多步法（7.3.1）是顯式的；若 ，則線性多步法（7.3.1）是隱式的.用線性多步法進行計算 時，除需要給定 外，還要附加初值 ，這可以用其他方法計算。由于多步法每計算一步用到的信息更多，因此希望由此構造出精度更高的算法,7.3.1 數值積分法,將微分方程 在 上積分，得,7.3.2,適當選取 個節點，作被積函數 的 次Lagrange插值多項式 并用 近似代替 ，就可得到形如式（7.3.1）的線性多步法。插值節點的不同取法就可得出不同的線性多步法,1. Adams外插法,取 為節點，構造 的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）的具體形式。假設前 個節點處的函數值已知，即 的近似值 已算出，從而函數值 也已知。這樣，就可以利用 個數據點 , ， 構造被積函數 的 次插值多項式,其中 是Lagrange插值基函數.從而可得,7.3.6,記 , 則有,這就是求解初值問題的Adams顯式公式。 是關于 和 的線性表達式，所以它是線性 步法,在上述Adams顯式公式的推導中，選用了 , 作為 插值節點，但構造的 次插值多項式 是代替區 間 上的未知函數 ，因此屬于

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

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

20、足線性方程組（7.3.10）的 ， 得到的線性多步法（7.3.1）的局部截斷誤差為 ， 可以證明此線性多步法的整體截斷誤差的階是 ，所以此線性多步法為 階 步法。顯然， 的大小和 有關,因為線性多步法（7.3.1）可以相差一個非零乘數，所以不妨設 。當 時， 可用 直接表示，稱為顯式法;反之，當 時，求 需要解一個方程（一般用迭代法），稱為隱式法.用待定系數法構造線性多步法的一個基本要求是選取 ， 使局部截斷誤差的階盡可能高,1. Milne方法,作為待定系數法的一個應用，下面討論一般的2步 法。此時 ， ，其余5個系數 由 確定，即,其中 為任意常數,所以，一般的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內插法；若取 ，則2步法（7.3.11）為顯式法,它的余項為,2. Hmaming方法,用待定系數法容易求出，4階3步的公式,這就是著名的4階Hmaming公式，它的余項為,更多常用的線性多步法和多步法計算中的問題等可以閱讀相關參考文獻，在此不再贅述.關于線性多步法的穩定性、收斂性和誤差估計，以及絕對穩定性和相對穩定性

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

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

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

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

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

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

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

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