1、1 1、常微分方程與解、常微分方程與解為為n n階常微分方程階常微分方程。0 ), , ,()(nyyyyxF如果函數如果函數 在區間在區間a,ba,b內內n n階可導，稱方程階可導，稱方程)(xyy )(xyy 滿足方程的函數滿足方程的函數稱為微分方程的稱為微分方程的解解。則則如如為任意常數）為任意常數）xy2 CCxy(, 2一般稱為方程的一般稱為方程的通解通解。為 方 程 的 解為 方 程 的 解。12 xy如果如果則有則有10 )(y為方程滿足定解條件的解。為方程滿足定解條件的解。第第10章章 常微分方程的數值解常微分方程的數值解10.1引引 言言 科學研究和工程技術中的問題往往歸結為

2、求某個常微分科學研究和工程技術中的問題往往歸結為求某個常微分方程的定解問題方程的定解問題. 常微分方程的理論指出，除少數簡單情況能獲得初值問常微分方程的理論指出，除少數簡單情況能獲得初值問題的初等解（用初等函數表示的解）外，絕大多數情況下是題的初等解（用初等函數表示的解）外，絕大多數情況下是求不出初等解的求不出初等解的.有些初值問題即便有初等解，也往往由于形有些初值問題即便有初等解，也往往由于形式過于復雜而不便處理。式過于復雜而不便處理。 常微分方程的數值解法常用來求近似解，由于它提供的常微分方程的數值解法常用來求近似解，由于它提供的算法能通過計算機便捷地實現，因此近年來得到迅速的發展算法能通

3、過計算機便捷地實現，因此近年來得到迅速的發展和廣泛的應用。和廣泛的應用。10.2 初值問題解法的基本概念 科學技術中常常需要求解常微分方程的定解問科學技術中常常需要求解常微分方程的定解問題題. 這類問題最簡單的形式，是本章將著重考察的這類問題最簡單的形式，是本章將著重考察的一一階方程的初值問題階方程的初值問題00( , ),(2.1)().yf x yy xy 我們知道，只有我們知道，只有f(x, y)適當光滑適當光滑譬如關于譬如關于y滿足滿足利普希茨利普希茨(Lipschitz)條件條件理論上就可以保證初值問題的解理論上就可以保證初值問題的解yf(x)存在并且唯一存在并且唯一.我們以下的討論

4、，都在滿足上述條件下進行。我們以下的討論，都在滿足上述條件下進行。( , )( , ).(2.2)f x yf x yL yy 雖然求解常微分方程有各種各樣的解析方法，但雖然求解常微分方程有各種各樣的解析方法，但解析方法只能用來求解一些特殊類型的方程，實際問解析方法只能用來求解一些特殊類型的方程，實際問題中歸結出來的微分方程主要靠數值解法題中歸結出來的微分方程主要靠數值解法. 所謂所謂數值解法數值解法, 就是尋求解就是尋求解y(x)在一系列離散節點在一系列離散節點 121nnxxxx上的近似值上的近似值 y1,y2,yn,yn+1,. 相鄰兩個節點的間距相鄰兩個節點的間距hn=xn+1- -x

5、n稱為稱為步長步長. 今后如不特別說明，總是假定今后如不特別說明，總是假定 hi=h(i=1,2,)為為定數定數, 這時節點為這時節點為xn=x0+nh(i=0,1,2,) (等距節點等距節點). 初值問題的初值問題的數值解法數值解法有個有個基本特點基本特點，他們都采，他們都采取取“步進式步進式”，即求解過程順著節點排列的次序一步，即求解過程順著節點排列的次序一步一步地向前推進一步地向前推進. 描述這類算法，只要給出用已知信描述這類算法，只要給出用已知信息息yn,yn- -1,yn- -2,計算計算yn+1的遞推公式的遞推公式. 首先，要對微分方程離散化，建立求解數值解的首先，要對微分方程離散

6、化，建立求解數值解的遞推公式遞推公式. 一類是計算一類是計算yn+1時時只用到只用到xn+1，xn和和yn，即，即前一步的值。因此，有了初值以后就可以逐步往下計前一步的值。因此，有了初值以后就可以逐步往下計算，其代表是龍格庫塔法算，其代表是龍格庫塔法稱為稱為單步法單步法. 另一類是用另一類是用到到yn+1前面前面 k 點的值點的值yn,yn-1, yn-k+1，稱為，稱為多步法多步法. 其其次，要研究公式的次，要研究公式的局部截斷誤差局部截斷誤差和和階階，數值解，數值解yn與精與精確解確解y(xn)的的誤差估計誤差估計及及收斂性收斂性，還有遞推公式的，還有遞推公式的計算計算穩定性穩定性等問題等

7、問題.數值解的思想數值解的思想（1 1）將連續變量）將連續變量 離散為離散為,bax bxxxxank 10nkxyykk,)(21 （2 2）用代數的方法求出解函數）用代數的方法求出解函數 在在 點的近似值點的近似值)(xyy kx)(kxy* *ky)(xyy 數學界關注數學界關注工程師關注工程師關注如果找不到解函數如果找不到解函數數學界還關注：數學界還關注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振動性解的振動性解的周期性解的周期性解的穩定性解的穩定性解的混沌性解的混沌性第一步：連續變量離散化第一步：連續變量離散化,nkxxxxx10第二步：用直線步進第二步：用直

8、線步進1 1、EulerEuler格式格式10.3 簡單單步法10.3.1 歐拉歐拉(Euler)方法方法00(,)xy000()(,)y xf xy0000(,)()yyf xyxx1xx100010(,)()yyf xyxx11()y xy11( ,)x y111( )( ,)y xf x y1111( ,)()yyf x yxx2xx211121( ,)()yyf x yxx22()y xy(,)nnxy()(,)nnny xf xy(,)()nnnnyyf xyxx過過 做以做以 為切線斜率的方程為切線斜率的方程 當當時，得時，得，取取 當當時，得時，得，取，取 過過做以做以為切線斜率

9、的方程為切線斜率的方程一般地，過一般地，過做以做以為切線斜率的方程為切線斜率的方程kp0p1p1npnpkx0 x1x1nxnx),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式nxx11(,)()nnnnnnyyf xyxx()nny xy當當時，得時，得，取取 (3.1) 例例1 用歐拉公式求解初值問題用歐拉公式求解初值問題2(01),(0)1.xyyxyy 解解 取步長取步長h=0.1，歐拉公式的具體形式為，歐拉公式的具體形式為)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得

10、由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy依次計算下去，依次計算下去，部分計算結果部分計算結果見下表見下表. xy21 與準確解與準確解 相比，可看出歐拉公式的計算結相比，可看出歐拉公式的計算結果精度很差果精度很差. xn 歐拉公式數值解歐拉公式數值解yn準確解準確解y(xn) 誤差誤差 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1

11、.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719 歐拉公式具有明顯的幾何意義歐拉公式具有明顯的幾何意義, , 就是就是用折線近似用折線近似代替方程的解曲線代替方程的解曲線，因而常稱公式，因而常稱公式(3.1)為為歐拉折線法歐拉折線法. .( )yy xxynx1nxnp1np1np x 還可以通過幾何直觀來考察歐拉方法的精度還可以通過幾何直觀來考察歐拉方法的精度. .假假設設yn=y(xn), ,即頂點即頂點Pn落在積分曲線落在積分曲線y=y(x)上，那么，上，那么，按歐拉方法做出的折線按歐拉方法做出的折線PnPn+1便

12、是便是y=y(x)過點過點Pn的切線的切線. .從圖形上看從圖形上看, ,這這樣定出的頂點樣定出的頂點Pn+1顯著顯著地偏離了原來的積分曲地偏離了原來的積分曲線，可見歐拉方法是線，可見歐拉方法是相相當粗糙當粗糙的的. .l 1818世紀最杰出的數學家之一，世紀最杰出的數學家之一，1313歲歲時入讀巴塞爾大學，時入讀巴塞爾大學，1515歲大學畢業，歲大學畢業，1616歲獲得碩士學位。歲獲得碩士學位。l 1727 1727年年-1741-1741年（年（2020歲歲-34-34歲）在彼歲）在彼得堡科學院從事研究工作，在分析學得堡科學院從事研究工作，在分析學、數論、力學方面均有出色成就，并、數論、力

13、學方面均有出色成就，并應俄國政府要求，解決了不少地圖學應俄國政府要求，解決了不少地圖學、造船業等實際問題。、造船業等實際問題。l 24 24歲晉升物理學教授。歲晉升物理學教授。l 1735 1735年（年（2828歲）右眼失明。歲）右眼失明。l 17411741年年 - 1766- 1766（3434歲歲-59-59歲）任德國科學院物理數學所所歲）任德國科學院物理數學所所長，任職長，任職2525年。在行星運動、剛體運動、熱力學、彈道學、人年。在行星運動、剛體運動、熱力學、彈道學、人口學、微分方程、曲面微分幾何等研究領域均有開創性的工作口學、微分方程、曲面微分幾何等研究領域均有開創性的工作。l

14、1766 1766年應沙皇禮聘重回彼得堡，在年應沙皇禮聘重回彼得堡，在17711771年（年（6464歲）左眼失歲）左眼失明。明。l Euler Euler是數學史上最多產的數學家，平均以每年是數學史上最多產的數學家，平均以每年800800頁的速頁的速度寫出創造性論文。他去世后，人們用度寫出創造性論文。他去世后，人們用3535年整理出他的研究成年整理出他的研究成果果7474卷。卷。 14 方法一化導數為差商的方法10()()()()()limnnnnnhy xhy xy xy xy xhh 由于在逐步求解的過程中，由于在逐步求解的過程中，y(xn) 的準確值無法求解的準確值無法求解出來，因此用

15、其近似值代替。出來，因此用其近似值代替。為避免混淆，以下學習簡記：為避免混淆，以下學習簡記：y(xn)：待求函數：待求函數 y(x) 在在 xn 處的精確函數值處的精確函數值yn ：待求函數：待求函數 y(x) 在在 xn 處的近似函數值處的近似函數值歐拉歐拉(Euler)方法（幾種推導法）方法（幾種推導法）15 代入初值問題表達式可得：根據根據 y0 可以一步步計算出函數可以一步步計算出函數 y y(x) 在在 x1, x2, x3 x4, 上的近似值上的近似值 y1, y2, y3, y4 , 常微分方程數值解是一組離散的函數值數據，它的常微分方程數值解是一組離散的函數值數據，它的精確表達

16、式很難求解得到，但可以進行插值計算后精確表達式很難求解得到，但可以進行插值計算后用插值函數逼近用插值函數逼近 y(x)10()()()()()limnnnnnhy xhy xy xy xy xhh 100(,)0,1,2,()nnnnyyhf xynyy x 1(,)nnnnyyf xyh 為了分析計算公式的精度，通?？捎锰├照归_為了分析計算公式的精度，通常可用泰勒展開將將y(xn+1)在在xn處展開，則有處展開，則有).,()(2),()()(2)()()()(1221 nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy 在在yn=y(xn)的前提下，的前提下，f(xn,y

17、n )=f(xn,y(xn)=y ( (xn n) ). .于是于是可得歐拉法可得歐拉法(3.1)的的公式誤差公式誤差為為2211()( )( ),(3.2)22nnnnhhy xyyy x稱為此方法的稱為此方法的局部截斷誤差局部截斷誤差. . 方法二泰勒級數展開法泰勒級數展開法17 方法三數值積分法數值積分法1111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x 同樣以近似值同樣以近似值 yn 代替精確值代替精確值 y(xn) 可得可得：100(,)()nnnnyyhf xyyy x 將微分方程將微

18、分方程 y f (x, y) 在區間在區間 xn, xn+1 上積分：上積分：182.隱式歐拉法隱式歐拉法(后退后退) 在數值積分法推導中，積分的近似值取為積分區間寬在數值積分法推導中，積分的近似值取為積分區間寬度與右端點處的函數值乘積，即：度與右端點處的函數值乘積，即：這樣便得到了隱式歐拉法：這樣便得到了隱式歐拉法：11100(,)()nnnnyyh f xyyy x 含有未知含有未知的函數值的函數值11111111d()() , ( )d() , (), ()nnnnxxnnxxnnnnnnyxy xy xf x y xxxxf xy xhf xy x （3.3） 隱式歐拉公式與歐拉公式有

19、著本質的區別隱式歐拉公式與歐拉公式有著本質的區別, 后者后者是關于是關于yn+1的一個直接計算公式，這類公式稱作是的一個直接計算公式，這類公式稱作是顯顯式的式的；前者公式的；前者公式的右端含有未知的右端含有未知的yn+1，它實際上是，它實際上是關于關于yn+1的一個函數方程的一個函數方程, ,這類方程稱作是這類方程稱作是隱式的隱式的. . 顯式顯式與與隱式隱式兩類方法各有特點，考慮到數值穩兩類方法各有特點，考慮到數值穩定性等其他因素，人們有時需要選用定性等其他因素，人們有時需要選用隱式隱式方法，但方法，但使用使用顯式顯式算法遠比算法遠比隱式隱式方便方便. . 隱式方程通常用迭代法求解，而迭代過

20、程的實隱式方程通常用迭代法求解，而迭代過程的實質是質是逐步逐步顯式化顯式化. . 設用歐拉公式設用歐拉公式),() 0(1nnnnyxhfyy 給出迭代初值給出迭代初值 ，用它代入，用它代入(3.1)式的式的右端，使之轉右端，使之轉化為顯式，直接計算得化為顯式，直接計算得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(3.1)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反復進行，得如此反復進行，得(1)( )111(,) (0,1, ).(3.4)kknnnnyyhf xyk由于由于f(x, y)對

21、對y滿足滿足Lipschitz條件條件(2.1). 由由(3.4)減減(3.3)得得),(),(11)(111) 1(1 nnknnnknyxfyxfhyy.1)(1 nknyyhL由此可知，只要由此可知，只要hL n) 上產生的擾動為上產生的擾動為 ，如果：，如果：md d0nd d (1,2,)mnmnndddd定義定義：設在節點：設在節點 xn 處用數值算法得到的理想數值解為處用數值算法得到的理想數值解為 yn，而實際計算得到的近似解為，而實際計算得到的近似解為 ，稱差值：，稱差值：ny nnnyyd d 為第為第 n 步的數值解的步的數值解的擾動擾動。則稱該數值方法是則稱該數值方法是穩

22、定穩定的。的。 下面以歐拉法為例考察計算穩定性下面以歐拉法為例考察計算穩定性. 例例4 用歐拉公式求解初值問題用歐拉公式求解初值問題 . 1)0(,100yyy 解解 用歐拉法解方程用歐拉法解方程y=- -100y 得得其準確解其準確解 是一個按指數曲線衰減很快的是一個按指數曲線衰減很快的函數函數.xexy100)( .)1001(1nnyhy 若若取步長取步長h=0.025，則歐拉公式的具體形式為，則歐拉公式的具體形式為.5 . 11nnyy 節點節點xn歐拉方法歐拉方法yn后退歐拉方法后退歐拉方法yn0.0250.0500.0750.100- -1.5 2.25 - -3.375 5.06

23、250.28570.08160.02330.0067計算結果見表計算結果見表, , 明顯計算過程不穩定明顯計算過程不穩定, , 但取但取h=0.005, , yn+1=- -1.5yn, , 則計算過程穩定則計算過程穩定. .對后退的歐拉公式，取對后退的歐拉公式，取h=0.025時，則計算公式時，則計算公式為為yn+1=- -(1/3.5)yn . .計算結果見表計算結果見表, , 這時計算過程是這時計算過程是穩定的穩定的. . 例題表明穩定性不但例題表明穩定性不但與方法有關與方法有關，也，也與步長與步長h有有關，當然關，當然與方程中的與方程中的f(x, y)有關有關. 為了只考察數值方為了只

24、考察數值方法本身，通常只檢驗數值方法用于解法本身，通常只檢驗數值方法用于解模型方程模型方程的穩的穩定性，定性，模型方程模型方程為為,(4.4)yy 其中其中為一直實數或復數為一直實數或復數(Re()0) ，這個方程分析，這個方程分析較簡單，對一般方程可以通過局部線性化化為這種較簡單，對一般方程可以通過局部線性化化為這種形式。形式。53 定義定義6 6 單步法（單步法（4.24.2）用于解模型方程（）用于解模型方程（4.44.4），若得），若得到的解到的解 ，滿足，滿足 ，則稱方法，則稱方法（4.14.1）是）是絕對穩定絕對穩定的的. . 1()nnyE hy1)(hE 在在 的平面上，使的平面

25、上，使 的變量圍成的的變量圍成的區域，稱為區域，稱為絕對穩定域絕對穩定域，h1)(hE它與實軸的交稱為它與實軸的交稱為絕對穩定區間絕對穩定區間. . 54歐拉法：歐拉法：1(,)nnnnyyhf xy 考察模型方程：考察模型方程：(0)yy 1(1)nnyhy 即：即：假設在節點值假設在節點值 yn 上有擾動上有擾動 d dn，在節點值，在節點值 yn 1 上有上有擾動擾動 d dn 1，且，且 d dn 1 僅由僅由 d dn 引起（即：計算過程中引起（即：計算過程中不再引起新的誤差）不再引起新的誤差）55111(1)(1)(1)()(1)nnnnnnnnyyhyhyhyyhd d d d

26、歐拉法穩定歐拉法穩定1nndddd 11h 即：即：111h 歐拉法穩定的條件：歐拉法穩定的條件：20h 0 1(1)nnyhy 針對模型方程：針對模型方程：的顯式歐拉法：的顯式歐拉法：1(,)nnnnyyhf xy 化簡得：化簡得：20h 56隱式歐拉法：隱式歐拉法：111(,)nnnnyyhf xy111111nnnnnnyyyyhhhd dd d 考察模型方程：考察模型方程：(0)yy 即：即：11nnnyyh y 化簡為：化簡為：11nnyyh 假設假設 yn 上有擾動上有擾動 ，則，則 yn 1 的擾動為：的擾動為：nnnyyd d 隱式歐拉法穩定隱式歐拉法穩定1nndddd 111

27、h 0 0h ，上式均成立，所以：上式均成立，所以：隱式歐拉法穩定是恒穩定的隱式歐拉法穩定是恒穩定的57,21211nnyhhy故故 .2/12/1)(hhhE對對 有有 ，0)Re( h 12/12/1)(hhhE故絕對穩定域為故絕對穩定域為 的左半平面，的左半平面，hyy 梯形法的穩定性梯形法的穩定性絕對穩定區間為絕對穩定區間為 ，即，即 時梯形法均時梯形法均是穩定的是穩定的. . 0h h0 隱式歐拉法與梯形方法的絕對穩定域均為隱式歐拉法與梯形方法的絕對穩定域均為 在具體計算中步長在具體計算中步長 的選擇只需考慮計算精度及迭代收斂的選擇只需考慮計算精度及迭代收斂性要求而不必考慮穩定性，具

28、有這種特點的方法特別重要性要求而不必考慮穩定性，具有這種特點的方法特別重要. . ,0)Re(hhh8.3 龍格龍格-庫塔庫塔(Runge-Kutta)法法歐拉方法是顯式的一步法歐拉方法是顯式的一步法,使用方便使用方便,但精度較低但精度較低.本節將構造出本節將構造出高精度的顯式一步法高精度的顯式一步法：龍格：龍格-庫塔法庫塔法,簡稱簡稱R-K法法.8.3.1 二階二階R-K法法歐拉法的公式為：歐拉法的公式為：yi+1=yi+h f (xi ,yi) i=0,1,2,n-1決定其精度的是函數決定其精度的是函數f (xi ,yi).如能改進這個函數如能改進這個函數,就可能提高公就可能提高公式的精度

29、式的精度.為此把公式改寫成：為此把公式改寫成：yi+1=yi+h (xi ,yi,h) i=0,1,2,n-1 (8.10)選擇函數選擇函數 (xi ,yi,h),一種方法是用若干個點的函數值的線性組合一種方法是用若干個點的函數值的線性組合代替代替 (xi ,yi,h),如：如：1111, 3 , 2),(),(),(jlljlijijiipjjjiipjKbhyhaxfKyxfKKchyx其中其中cj,aj,bjl是待定參數是待定參數, aj和和bjl滿足滿足pjbajljlj,3,211以上方法稱為以上方法稱為p級級R-K法法,選擇選擇cj,aj和和bjl,可能使以上方法為可能使以上方法為

30、p階階方法方法.顯然歐拉法就是一階顯然歐拉法就是一階R-K法法.二級二級R-K法的形式是法的形式是:),(),(1222122111hKayhaxfKyxfKKcKchyyiiiiii此時此時)11. 8 ()()(,()()(,()(,()(,(),(22212221iiiiiiiiiiiixyhaxyhaxfcxycxyxhfaxyhaxfcxyxfchyx由二元函數的泰勒展開：由二元函數的泰勒展開：)()()(,()()(,(22222hOfxyhafhaxyxfxyhaxyhaxfyixiiiii其中所有的偏導數都是它們在點其中所有的偏導數都是它們在點(xi,y(xi)的值的值,下同下

31、同又由于：又由于：)()()(,()(xyffxyxyxfxyyx 所以所以)()()()()(,(2222hOxyhaxyxyhaxyhaxfiiiii 代入代入(8.11)()()(),(2221hOxyhaxycchyxiiii 代入代入(8.10)()()()(3222211hOxyhcaxyhccxyyiiii 而而Taylor展開式展開式)()(2)()()(321hOxyhxyhxyxyiiii 二式相減二式相減,得局部截斷誤差得局部截斷誤差)()(21)(1)(322221111hOxyhcaxyhccyxyRiiiii :021,012221得令cacc2112221cacc

32、只要只要c1,c2,a2滿足以上方程滿足以上方程,就得到一個二階的就得到一個二階的R-K法法.這是一個不定方程這是一個不定方程,有無窮多解有無窮多解.比如：比如：63結束結束(1)取取c1=c2=1/2,a2=1得得)12. 8 (),(),(2121211hKyhxfKyxfKKKhyyiiiiii這實際上是這實際上是(8.9)公式公式,即梯形公式的預估即梯形公式的預估-校正公式只迭代一次校正公式只迭代一次的形式的形式,通常稱為通常稱為改進的歐拉法改進的歐拉法.(2)取取c1=0,c2=1,a2=1/2 得得:)13. 8 (2,2),(12121KhyhxfKyxfKhKyyiiiiii這

33、公式又稱這公式又稱中點公式中點公式.我們還可以構造其他的二階我們還可以構造其他的二階R-K法法.8.3.2 四階四階 R-K 法法用類似的方法可以確定三級和四級用類似的方法可以確定三級和四級R-K法的參數法的參數,構造出三階和構造出三階和四階的四階的R-K法法.但最常用的是四階但最常用的是四階R-K法法,四階四階R-K法也不只一個法也不只一個,下面給出的是最常用的四階經典的下面給出的是最常用的四階經典的R-K公式：公式：)14. 8 (1, 2 , 1 , 0,2,22,2),(226342312143211nihKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyyiiiiiiiii

34、i65結束結束例例3 用經典的四階用經典的四階R-K法計算例法計算例2題目題目,取步長為取步長為0.2,且與準確值比且與準確值比較較.計算結果列入表計算結果列入表8-3：可見即使用可見即使用h=0.2計算計算,也比一階和二階方法精度好得多也比一階和二階方法精度好得多66結束結束8.4 線性多步法線性多步法單步法只利用前一步的結果單步法只利用前一步的結果,只要給出初值只要給出初值,就能開始計算就能開始計算但也因但也因為它只利用前一步的值為它只利用前一步的值,為了提高精度就要計算一些非結點處的函為了提高精度就要計算一些非結點處的函數值數值,增加了計算量增加了計算量.R-K法就是通過這一途徑提高精度

35、的法就是通過這一途徑提高精度的.下面介下面介紹的線性多步法紹的線性多步法,在求在求yi+1時時,不僅用到不僅用到yi的值的值, 還用到前若干步的還用到前若干步的yi-1,yi-k的值的值,這些值都是已知的這些值都是已知的,因此可在計算量增加不多的情況下因此可在計算量增加不多的情況下提高精度提高精度.8.4.1 用待定系數法構造線性多步法用待定系數法構造線性多步法線性多步法的一般形式是：線性多步法的一般形式是：kikiikikiifffhyyy110110kjjijkjjijknifhy00)15. 8(., 2 , 1 , 0,或寫為：或寫為：其中其中j,j(j=0,1,k)都是實常數都是實常

36、數,且且k0, 0+00, fi+j= f (xi+j,yi+j), j=0,1,k,由由(8.15)可看出要計算可看出要計算yi+k,要利用它前面的要利用它前面的k個值個值yi,yi+1,yi+k-1,又又因為因為(8.15)關于關于yi+j和和 f i+j都是線性組合都是線性組合,所以這一類方法都稱為線所以這一類方法都稱為線性性k步法步法.歐拉法歐拉法,隱式歐拉法和梯形法都是線性一步法隱式歐拉法和梯形法都是線性一步法,歐拉中點歐拉中點公式是線性二步法公式是線性二步法.,)()(),(0kjjjjhxyhjhxyhxyL令把把y(x+jh)和和y(x+jh)作作Taylor展開：展開： ,)

37、(!)(! 2)()()()(2 xypjhxyjhxyjhxyjhxypp ,)()!1()()()()(1 xypjhxyjhxyjhxypp68結束結束,)()()()()!1(!)()(),(,)(100)(1xyhCxyhCxyCxyhpjpjxyhjxyhxyLpppkjppjpjpjjj于是其中其中kC100 kkkC102112kkkkC21221212222 kkkkC2212162321613222)16. 8 ()!1/(2!/2121121pkpkCkppkppp69結束結束若選擇若選擇j,j,使使C0=C1=Cp=0,Cp+10,則則)()()()(),(2)1(11

38、0ppppkjjjhOxyhCjhxyhjhxyhxyL將將x=xi代入上式代入上式,設設k=1,并注意到并注意到y (xi+jh)=f(xi+jh,y(xi+jh),可推可推出出)()()(,()()(2)1(11010ppppkjiijkjijkihOxyhCjhxyjhxfhjhxyxy即即)17.8()()()(,()()(2)1(11010ppppkjjijijkjjijkihOxyhCxyxfhxyxy70結束結束設設yi=y(xi),yi+1=y(xi+1),yi+k-1=y(xi+k-1),記左端為記左端為yi+k,并舍去最后兩項并舍去最后兩項,(8.17)變為：變為：)18.

39、8(),(010kjjijijkjjijkiyxfhyy就是一種就是一種p階的線性階的線性k步方法步方法,Cp+1hp+1y(p+1)(xi)稱為局部截斷誤差稱為局部截斷誤差的主項的主項.當當k=0時時,是顯式方法是顯式方法,當當k0時是隱式方法時是隱式方法.下面構造幾個實用的線性多步法公式下面構造幾個實用的線性多步法公式例例4 形如：形如：yi+4= -0yi+h(1fi+1+2fi+2+3fi+3)的線性的線性4步法公式步法公式,試確定試確定0,1,2,3并求其局部截斷誤差主項并求其局部截斷誤差主項.71結束結束解解: 由由(8.15)知知1=2=3=0,4=1,0=4=0,因為有因為有4

40、個待定系數個待定系數,由由(8.16)寫出前寫出前4個方程：個方程：0!232! 340)32(!240)(401322213332122321100CCCC解之解之,得得0=-1,1=3=8/3,2=-4/3故所求的公式為故所求的公式為)19.8()22(343214iiiiifffhyy72結束結束將將0,1,2,3代入代入C40!332!443323144C再求再求C54514!432!543424155C公式公式(8.19)稱為稱為米爾尼米爾尼(Milne)公式公式,它的局部截斷誤差為：它的局部截斷誤差為：)()(45146)5(5hOxyhRi它也可以寫成它也可以寫成)21.8()(

41、4514)5(5yhR )20. 8 (1, 3 , 2)22(341231nifffhyyiiiii73結束結束例例5 試確定下列公式的系數和局部截斷誤差試確定下列公式的系數和局部截斷誤差yi+1=a yi+b yi-2+h (c fi+1+d fi+e fi-1)解解 按按(8.15)和和(8.16)可知可知k=3,0= - b,1=0,2= - a,3=1,0=0,1=e,2=d,3=c,有五個未知參數有五個未知參數,寫出前五個方程：寫出前五個方程：0!332!4230!232!3230)32(!2230)(3201334442233322210cdeaCcdeaCcdeaCcdeaCb

42、aC74結束結束解之得：解之得：a=9/8, b= - 1/8, c=3/8, d=6/8, e= - 3/8代入代入C5401! 432! 52344555cdeaC此公式稱為此公式稱為哈明哈明(Hamming)公式公式,寫為：寫為：)22. 8 (1, , 3 , 2)2(83)9 (811121nifffhyyyiiiiii局部截斷誤差記為：局部截斷誤差記為：)23.8()(401)5(5yhR它是一個四階隱式三步法它是一個四階隱式三步法.在在(8.15)中若中若0=1=k-2=0,k=1,則稱此類方法為阿達姆斯則稱此類方法為阿達姆斯(Adams)類方法類方法.當當k=0時時,稱顯式稱顯

43、式,否則為隱式否則為隱式.這類方法也可以這類方法也可以用以上待定系數法求出用以上待定系數法求出這類方法的通式為：這類方法的通式為：75結束結束)24.8(011kjjkijiifhyy例例6 求求(8.24)中的隱式三步法中的隱式三步法解解:k=3,0=1=0,2=-1,3=10632243202326320)32(2320)(32033231342443222133233321322223210321320CCCCC76結束結束46527831994253213213213213210解之解之,0=1/24,1=-5/24,2=19/24,3=9/24.)25. 8 (1, , 3 , 2)

44、9195(241121niffffhyyiiiiii公式為：公式為：72019)32(241)32(12013424135255C77結束結束)26.8()(72019)5(5yhR局部截斷誤差為局部截斷誤差為它是一個線性三步四階隱式公式它是一個線性三步四階隱式公式,稱為稱為四階阿達姆斯內插公式四階阿達姆斯內插公式,應應用十分普遍用十分普遍.此外常用的還有線性四步此外常用的還有線性四步四階的顯式阿達姆斯公式四階的顯式阿達姆斯公式,常稱為阿達常稱為阿達姆斯外推公式：姆斯外推公式：)28.8()(720251)5(5yhR )27. 8 (1, , 3 , 2)9375955(243211niff

45、ffhyyiiiiii78結束結束8.4.2 用數值積分法構造線性多步法公式用數值積分法構造線性多步法公式在在(8.15)式中式中,若令若令k=1,0,1,k-1中只有一個不為零中只有一個不為零,這樣的這樣的線性多步法公式也可用數值積分方法構造線性多步法公式也可用數值積分方法構造.如阿達姆斯外推公式如阿達姆斯外推公式.)2.8()()1.8()(,(00yyfyxxyx對對(8.1)兩端在兩端在xi,xi+1上積分上積分)29.8()(,()()(11dxxyxxyxyiixxiif選擇節點選擇節點xi-3,xi-2,xi-1,xi作插值節點作函數作插值節點作函數f (x,y(x)的三次插值多

46、的三次插值多項式：項式：79結束結束3212311323213)3)(2()()()(2()()()3)()()()3)(2()()()(iiiiiiiiiiiiiiiifhhhxxxxxxfhhhxxxxxxfhhhxxxxxxfhhhxxxxxxxP且插值余項且插值余項)()()(!4)(321)4(3iiiixxxxxxxxfR把把f (x,y (x)=P3(x)+R3(x)代入代入(9.29)得：得：dxxRdxxPxyxyiiiixxxxii11)()()()(331經計算積分后得經計算積分后得3213937595524)(1iiiixxffffhdxxPii80結束結束)31.9(

47、9375955243211iiiiiiffffhyy所以得所以得就是阿達姆斯外推公式就是阿達姆斯外推公式1)()()(!4)(321)4(iixxiiiidxxxxxxxxxfR其局部截斷誤差為其局部截斷誤差為用積分中值定理依然可以推出用積分中值定理依然可以推出:)(720251)(720251)()()(!4)()5(5)4(5321)4(1yhfhdxxxxxxxxxfRiixxiiii81結束結束)32.8(1, 3 ,2431111nifffhyyiiiii用類似的方法還可以構造二步四階的辛普生用類似的方法還可以構造二步四階的辛普生(Simpson)公式：公式：)33.8()(901)

48、5(5yhR它的局部截斷誤差為它的局部截斷誤差為線性多步法要用到二個以上的初始值才能開始計算線性多步法要用到二個以上的初始值才能開始計算,不能不能“自自啟動啟動”.一般應選擇同階或高階的單步法算出前幾個值一般應選擇同階或高階的單步法算出前幾個值,再使用再使用線性多步法線性多步法.由于線性多步法每一步實際上只計算一次函數值由于線性多步法每一步實際上只計算一次函數值,又能獲得較高的精度又能獲得較高的精度,所以應用比較廣泛所以應用比較廣泛.82結束結束8.5 預估校正法預估校正法在在8.2中介紹過梯形公式的預估校正法中介紹過梯形公式的預估校正法,本節中介紹幾種高階的本節中介紹幾種高階的實用的預估校正

49、法實用的預估校正法.介紹幾個英文縮寫：介紹幾個英文縮寫：P(Predictor)預估預估C(Corrector)校正校正E(Evaluation)求函數值求函數值M(Modifier)改進改進83結束結束8.5.1 阿達姆斯公式的阿達姆斯公式的PEC模式模式)34. 8 (519924),(93759552421)0(11)0(11)0(1321)0(1iiiiiiiiiiiiiiiffffhyyyxfffffhyy它是用阿達姆斯外推公式作預測它是用阿達姆斯外推公式作預測,用阿達姆斯內插公式作一次校正用阿達姆斯內插公式作一次校正例例7 用用Adams外推公式和外推公式和AdamsPEC方法計算

50、方法計算1)0(2yyxyy取步長取步長h=0.1,前三步用四階前三步用四階R-K法求值法求值,計算值保留七位有效數計算值保留七位有效數字字.84結束結束解解: 計算結果列表如下計算結果列表如下( 略略)8.5.2 阿達姆斯公式的阿達姆斯公式的PMECME模式模式令令Pi和和Ci分別為分別為(8.34)中的第中的第i步的預測值和校正值：步的預測值和校正值：)(72019)()(720251)()5(511)5(511yhCxyyhPxyiiii二式相減二式相減,設設y(5)(x) 變化不大變化不大,都約等于都約等于y(5)()(720270)(720251)(72019)5(5)5(5)5(5

51、11yhyhyhPCii11)5(5270720)(iiPCyh有85結束結束)(27019)(27072072019)()(270251)(270720720251)(1111iiiiiiiiiiiiPCPCCxyPCPCPxy由此得：由此得：于是在于是在(8.34)式中加上兩個改進步式中加上兩個改進步(M),最后再求一次函數值最后再求一次函數值(E)供下一次預測用供下一次預測用,就形成阿達姆斯法的就形成阿達姆斯法的PMECME模式模式)35. 8 (1, 4 , 3),(:)(27019:519924:),(:)(270251:937595524:11111112111111113211n

52、iyxffPCCyfffMhyCMxfMPCPMffffhyPiiiiiiiiiiiiiiiiiiiiiiiiii計算改進校正計算改進預測86結束結束開始求開始求M4時時,無無C3和和P3,認為認為C3-P3=0 8.5.3 哈明哈明(Hamming)法法PMECME模式模式利用米爾尼公式利用米爾尼公式(8.20)進行預測進行預測,再用哈明公式再用哈明公式(8.22)進行校正進行校正,中間中間利用兩公式局部截斷誤差的主項系數的比例進行改進兩次利用兩公式局部截斷誤差的主項系數的比例進行改進兩次,就得到就得到另一種另一種PMECME模式模式,此模式稱此模式稱哈明法哈明法,哈明法為：哈明法為：)36. 8(1, 4 , 3),(:)(1219:283981:),(:)(121112:2234:11111111121111112131niyxffCPCyffMhyyCMxfMCPPMfffhyPiiiiiiiiiiiiiiiiiiiiiiiii計算改進校正計算改進預測開始求開始求M4時時,無無C3和和P3,認為認為P3-C3=087結束結束8.6 一階常微分方程組和高階方程一階常微分方程組和高階方程本章本章8.1.2已提到已提到,一