1、 8 常微分方程的數值解法常微分方程的數值解法 （ Numerical Solution of Ordinary Differential Equations ） n 本章主要內容本章主要內容n 8.1 基本概念與基本求解途徑基本概念與基本求解途徑 n 8.2 歐拉方法與局部截斷誤差歐拉方法與局部截斷誤差 n 8.3 龍格龍格庫塔法庫塔法n 8.4 單步法的收斂性與穩定性單步法的收斂性與穩定性n 8.5 線性多步法線性多步法 n 8.6 一階常微分方程組與高階常微分方程的數值解法一階常微分方程組與高階常微分方程的數值解法n 8.7 邊值問題的差分方法簡介邊值問題的差分方法簡介 n 重點：重點：

2、歐拉法與龍格歐拉法與龍格庫塔法、預測庫塔法、預測校正技術校正技術n 難點：難點：各種方法的收斂性與穩定性各種方法的收斂性與穩定性8.1 基本概念與基本求解途徑基本概念與基本求解途徑 常常微分方程的應用背景微分方程的應用背景 00)()(qtqRCtqUdtdq 橋梁擾度問題設一根長為L的兩端固定的矩形截面梁，q是均勻負荷強度，E是彈性模量，S是斷點作用力，I(x)是慣性矩。則描述梁的橈度y(x) 的數學模型是 電容充電問題 設q是電容器上的帶電量，C為電容，R為電阻，U為電源的電動勢。描述電容器充電過程的數學模型是 牛頓力學定律設s是物體的位移，f是物體的受力，M是物體的質量，則是物體的速度，

3、a是物體的加速度，牛頓定律可描述為000022)(,)(vtsstsMfadtsd一階常微分方程一階常微分方程的初值問題的初值問題二階常微分方程二階常微分方程的初值問題的初值問題二階常微分方程的二階常微分方程的兩點邊值問題兩點邊值問題0)()0()()(2)()(22LyyLxxEIqxxyxEISdxyd8.1 基本概念與基本求解途徑基本概念與基本求解途徑 常常微分方程的數學模型微分方程的數學模型 一階常微分方程初值問題的數學模型的一般形式00)(),()(yxyyxfxy二階常微分方程初值問題的數學模型的一般形式兩點邊值問題（第一類邊界條件）的數學模型的一般形式 0000)(,)()()(

4、)()()(yxyyxyxfxqxyxpxy )(,)(),()(byayyyxfxy 數學方法及存在問題數學方法及存在問題 只有少數簡單的微分方程能夠用初等方法求得它們的解，多數情形只能利用近似方法求解。在常微分方程教材中的級數解法級數解法、逐步逼近法逐步逼近法等都是近似方法，這些方法可以給出解的近似表達式，統稱為近似解析近似解析法法。但實際中遇到的常微分方程問題往往不能給出解析解，或者解析解不變于使用。如下面兩個問題及解析解分別為0)0(21yxydxdy1)0(yydxdyxxtxeydteey022解析解要用插值法或要用插值法或數值積分法求值數值積分法求值要查指數表要查指數表8.1 基

5、本概念與基本求解途徑基本概念與基本求解途徑8.1 基本概念與基本求解途徑基本概念與基本求解途徑 微分方程的數值求解途徑的數學模型微分方程的數值求解途徑的數學模型 數值方法的基本思想是離散化，即把求解區間a,b分成n等分（區間離散化），同時對微分算子離散化，求函數y(x)在n+1個離散的等距節點上的近似值，而不必求y(x)的解析表達式。 初初值問題的數值解法，就是尋求解函數值問題的數值解法，就是尋求解函數y(x)在在a,b的一系列等距的的一系列等距的離散節點離散節點 a=x0 x1n)的絕對誤差均不超的絕對誤差均不超過過 ，則稱該單步法是穩定的。，則稱該單步法是穩定的。 例例8-7 考察顯歐拉方

6、法和隱歐拉方法的絕對穩定性。考察顯歐拉方法和隱歐拉方法的絕對穩定性。 例例8-8 考查二階龍格考查二階龍格-庫塔公式的絕對穩定性庫塔公式的絕對穩定性 。 例例8-9 用四階經典用四階經典R-K公式求解初值問題公式求解初值問題(h分別取分別取0.1和和0.2 )1)0( 1 , 020yxydxdy 表表8-3 常用單步法的絕對穩定區間常用單步法的絕對穩定區間在計算在計算yn+1時，選擇前面已算出的時，選擇前面已算出的yn-r,yn-r+1, ,yn-1,yn多個值的線性組多個值的線性組合作為初值，相應近似斜率合作為初值，相應近似斜率yn=f(xn,yn)的線性組合構造增量表達式，的線性組合構造

7、增量表達式，以獲得較高的精度。以獲得較高的精度。8.5 線性多步法線性多步法 線性多步法的基本思想線性多步法的基本思想 線性線性r+1步法計算公式的一般形式步法計算公式的一般形式8.5.1 線性多步法公式的構造線性多步法公式的構造101(8.14)rrnknkknkkkya yhb yb-1=0時，顯式公式時，顯式公式b-10時，隱式公式時，隱式公式 基于泰勒級數展開的構造方法基于泰勒級數展開的構造方法01111(8.19)()()1(1,2, )rkkrrjjkkkkakajkbjp)()()() 1()(1 )!1(2) 1(11111pnprkrkkpkppnhOxybkpakphR p

8、+1個方程個方程,2r+1個未知數的線個未知數的線性方程組，其解解即可構造出性方程組，其解解即可構造出步步p階精度的線性多步公式階精度的線性多步公式 局部截斷誤差局部截斷誤差8.5 線性多步法線性多步法 8.5.1 線性多步法公式的構造線性多步法公式的構造 基于泰勒級數展開的構造方法基于泰勒級數展開的構造方法例例8-10 構造形如構造形如 的線性兩步隱式公式的線性兩步隱式公式 )(110111101nnnnnnfbfbfbhyayay例例8-11 構造形如構造形如 的線性兩步隱式公式的線性兩步隱式公式 )()(11011nnnnnfbfbhyyay8.5 線性多步法線性多步法 8.5.1 線性

9、多步法公式的構造線性多步法公式的構造 基于數值積分的構造方法基于數值積分的構造方法將初值問題（將初值問題（8-1）式的微分方程兩端同時從）式的微分方程兩端同時從xn-r到到xn+1積分，得積分，得1)(,()()(1nrnxxrnndxxyxfxyxy 對式中的被積函數對式中的被積函數f(x,y(x)用用Lagrang插值函數近似，積分后可得插值函數近似，積分后可得到不同的線性到不同的線性r+1步公式。在區間步公式。在區間xn-r,xn+1選取的插值節點越多，選取的插值節點越多，公式的精度越高，插值節點若包括公式的精度越高，插值節點若包括xn+1則為則為隱公式隱公式，否則為，否則為顯公式顯公式

10、。課后作業：課后作業：使用泰勒級數法使用泰勒級數法構造構造Milne公式公式8.5 線性多步法線性多步法 8.5.2 亞當姆斯亞當姆斯(Adams)公式公式 Adams顯式公式顯式公式（外插公式）（外插公式）r步步r階階Adams顯式公式的一般形式為顯式公式的一般形式為二階二階Adams顯公式顯公式三階三階Adams顯公式顯公式四階四階Adams顯公式顯公式101rkknknnybhyy Adams顯式公式既可以通過數值積分來構造，也可以用顯式公式既可以通過數值積分來構造，也可以用Taylor展開展開式來構造?，F通過數值積分來構造可得式來構造?，F通過數值積分來構造可得10101),(rjjnj

11、nnrjjnjnjnnnfAyyxfAyy1)(nnxxjnjndxxlA其中)3(211nnnnffhyy)51623(12211nnnnnfffhyy)9375955(243211nnnnnnffffhyy表表8-5顯式顯式Adams公式系數表公式系數表8.5 線性多步法線性多步法 8.5.2 亞當姆斯亞當姆斯(Adams)公式公式 Adams隱式公式隱式公式（內插公式）（內插公式）r步步r+1階階Adams隱式公式的一般形式為隱式公式的一般形式為211rkknknnybhyy Adams隱式公式通常隱式公式通常Taylor展開式來構造。展開式來構造。3步步4階階Adams隱式公式隱式公式

12、表表8-6隱式隱式Adams公式系數表公式系數表)2()()()()()()(2101111hxybhxybxybhxybhxyhxyyxyRnnnnnnnnn)()()3462241()()22261()()221()()1 (54)4(211321122112101hOhxybbbhxybbbhxybbbhxybbbbnnnn 為達到四階精度，令為達到四階精度，令h的的14次冪的系數為零，解之得次冪的系數為零，解之得241,245,2419,832101bbbb)5199(242111nnnnnnffffhyy從而有從而有3步步4階階Adams隱式公式隱式公式8.5 線性多步法線性多步法

13、8.5.3 線性多步法預測線性多步法預測校正公式校正公式 不論是單步法還是多步法，隱公式相比顯公式而言其穩定性好，但在不論是單步法還是多步法，隱公式相比顯公式而言其穩定性好，但在計算時會遇到兩個問題：計算時會遇到兩個問題：n 隱公式如何能方便的計算；隱公式如何能方便的計算；n 實際計算步長取多大。實際計算步長取多大。 一般地，不難驗證，如果預測公式（顯公式）是一般地，不難驗證，如果預測公式（顯公式）是p階或階或p+1階精度，校階精度，校正公式（隱公式）是正公式（隱公式）是p+1階精度，則用預測公式提供初值，校正公式迭代階精度，則用預測公式提供初值，校正公式迭代一次的效果也能達到一次的效果也能達

14、到p+1階精度，再迭代下去的話，效果就不明顯了。階精度，再迭代下去的話，效果就不明顯了。預預測測校正技術校正技術既保證了計算精度，又使隱式計算顯式化，克服了隱式公式既保證了計算精度，又使隱式計算顯式化，克服了隱式公式需要反復迭代計算的困難。通過對預測需要反復迭代計算的困難。通過對預測校正公式使用校正公式使用外推原理，外推原理，得到誤得到誤差估計式，用來調整計算步長，使其達到精度要求。差估計式，用來調整計算步長，使其達到精度要求。 4階階Adams顯公式顯公式預測預測，4階階Adams隱公式隱公式校正校正),(9),(37),(59),(55(243322111nnnnnnnnnnyxfyxfy

15、xfyxfhyy),(),(5),(19),(9 (242211111nnnnnnnnnnyxfyxfyxfyxfhyy8.5 線性多步法線性多步法 8.5.3 線性多步法預測線性多步法預測校正公式校正公式 Milne(4階階)公式公式預測預測，Hamming(4階階)公式公式校正校正),(2),(),(2(34221131nnnnnnnnyxfyxfyxfhyy),(),(2),(83)9(81111121nnnnnnnnnyxfyxfyxfhyyy 中點公式公式中點公式公式預測預測，梯形公式，梯形公式校正校正112(,)nnnnyyhf xy111 (,)(,)2nnnnnnhyyf xy

16、f xy由局部截斷誤差可知由局部截斷誤差可知 進而得到進而得到3階公式階公式)(51),(),(2),(2111111111nnnnnnnnnnnnnnyyyyyxfyxfhyyyxhfyy改進：校正：預測：注意：注意：隱式公式的計算需要隱式公式的計算需要用單步法（如龍格用單步法（如龍格-庫塔法）庫塔法）若干個初值才能啟動計算若干個初值才能啟動計算。)()(12)()()(3)(4311143111hOxyhyxyThOxyhyxyTnnnnnnnn )(414)(4111hOTTxynnn 改進公式是一個改進公式是一個事后誤事后誤差估計式差估計式。通常用下式。通常用下式來選擇合適的步長。來選

17、擇合適的步長。)(5111nnyy8.6 一階常微分方程組和高階常微分方程簡介一階常微分方程組和高階常微分方程簡介8.6.1 一階常微分方程組一階常微分方程組 一階常微分方程組初值問題的一般形式一階常微分方程組初值問題的一般形式 一階常微分方程組初值問題的解法一階常微分方程組初值問題的解法bxacaycaycayyyyxfxyyyyxfxyyyyxfxynnnnnnn)(,)(,)(),()(),()(),()(22112121222111TnTnTnccccyxfyxfyxfyxfxyxyxyxy),(),(,),(),(),()(,),(),()(212121cayyxfxy)(),()(

18、定義定義向量向量得到得到 類似于一階常微分方程初值問題的數值方法，但需注意求解公式也是用向量形式表示。例例8-12 二階亞當姆斯預測二階亞當姆斯預測校正公式求解初值問題方程組校正公式求解初值問題方程組8.6.2 高階常微分方程高階常微分方程 高階常微分方程初值問題的一般形式高階常微分方程初值問題的一般形式bxacaycaycayyyyxfxymmmm1) 1(10) 1()()(,)(,)(),()()1(21,mmyyyyyymjcayyyxfyyyyyyyjjmmmm,2, 1)(),(1113221引入變量：轉化為一階常微分方程組的初值問題 高階常微分方程初值問題的解法高階常微分方程初值

19、問題的解法求解一階常微分方程組的初值問題例例 用歐拉法求解2階常微分方程的初值問題。) 10(1)0(, 0)0(2 xyyyyxy8.7 邊值問題的數值解簡介邊值問題的數值解簡介8.7.1 2 階常微分方程邊值問題的階常微分方程邊值問題的打靶法打靶法 2 階常微分方程邊值問題的一般形式階常微分方程邊值問題的一般形式 打靶法的基本思想打靶法的基本思想 )(,)(),(),(byaybaxyyxfy先猜測一個初始斜率 y (a) = s，通過解初值問題 sayaayyyxfy)()(),(找出s*使得(s*) = ，即把問題轉化為求方程 (s) = 0 的根。yx0aby x( ) 斜率斜率 = s0 ()s0斜率斜率 = s1 ()s1每每計算一個計算一個 ( (s s) ) 都必須解一個初值問題都必須解一個初值問題y(b) = (s)8.7.2 2 階常微分方程邊值問題的階常微分方程邊值問題的有限差分法有限差分法 有限差分法的基本思想有限差分法的基本思想將求解區間a,b 等分為N 份，取節點 xi = a +ih (i = 0,N )，在每一個節點處將 y 和 y 離散化。 有限差分法公式推導有限差分法公式