第九章常微分方程數(shù)值解法

1、第第9 9章章 常微分方程數(shù)值解法常微分方程數(shù)值解法 /* Numerical Methods for Ordinary Differential Equations */ 考慮考慮一階一階常微分方程的常微分方程的初值問題初值問題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上連續(xù)，且關(guān)于上連續(xù)，且關(guān)于 y 滿足滿足 Lipschitz 條條件件，即存在與，即存在與 x, y 無關(guān)的常數(shù)無關(guān)的常數(shù) L 使使對任意定義在對任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則上

2、述都成立，則上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 理論上雖然已經(jīng)比較完善，但許多實(shí)際問題比如理論上雖然已經(jīng)比較完善，但許多實(shí)際問題比如, , ()200 11xdyexdxy x是無法用解析方法求解的，而本章將研究比較通用的微分方程是無法用解析方法求解的，而本章將研究比較通用的微分方程數(shù)值解法。數(shù)值解法。要計算出解函數(shù)要計算出解函數(shù) y(x) 在一系列節(jié)點(diǎn)在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值處的近似值),., 1()(nixyyii 節(jié)點(diǎn)間距節(jié)點(diǎn)間距 為步長，通常采用為步長，通常采用等距節(jié)點(diǎn)等距節(jié)點(diǎn)，即取即取 hi = h (常數(shù)常數(shù)

3、)。于是。于是) 1,., 0(1 nixxhiii所謂數(shù)值解法，就是求在一系列離散節(jié)點(diǎn)所謂數(shù)值解法，就是求在一系列離散節(jié)點(diǎn)( )yf x011nnxxxx上的近似值上的近似值,011nnyyyy,(),(,)()(, ,)01 2 nnnnnnnxxnh yy xff xyy xn初值問題數(shù)值解法的基本思想是初值問題數(shù)值解法的基本思想是”離散化離散化“，”步進(jìn)式步進(jìn)式“，即求，即求解過程依節(jié)點(diǎn)的排列次序一步一步向前推進(jìn)，一般有兩種做法：解過程依節(jié)點(diǎn)的排列次序一步一步向前推進(jìn)，一般有兩種做法：單步法：單步法：計算時只用到前一點(diǎn)的值。計算時只用到前一點(diǎn)的值。多步法多步法：用到前面多個點(diǎn)的值。：用

4、到前面多個點(diǎn)的值。 差商逼近法差商逼近法是用適當(dāng)?shù)牟钌瘫平鼘?dǎo)數(shù)值使方程離散化的一種方法。如取是用適當(dāng)?shù)牟钌瘫平鼘?dǎo)數(shù)值使方程離散化的一種方法。如取()()()()(),()111nnnnnny xy xy xy xy xy xhh引入?yún)?shù)由上式得引入?yún)?shù)由上式得 , 0 1 ()()()()(),111nnnny xy xy xy xh()()()()()111nnnny xy xhy xy x()111nnnnyyhff取取1 1nnnyyhf顯式顯式Euler法法取取 0 11nnnyyhf隱式隱式Euler法法取取12 ()112nnnnhyyff梯形法梯形法 數(shù)值積分法數(shù)值積分法基本思想

5、是將原微分方程轉(zhuǎn)化為積分方程基本思想是將原微分方程轉(zhuǎn)化為積分方程()()( , ( ),(),00mnxmnxy xy xf x y x dx y xy mn然后用數(shù)值積分法離散化，得到一個離散的差分格式，如利用矩然后用數(shù)值積分法離散化，得到一個離散的差分格式，如利用矩形公式，同時取形公式，同時取m=n+1則得則得(),(), , 111110 1nnnnnnyyhfxxyy取取1 1nnnyyhf取取 0 11nnnyyhf1 歐拉方法歐拉方法 /* Eulers Method */ 歐拉公式：歐拉公式：x0 x1向前差商近似導(dǎo)數(shù)向前差商近似導(dǎo)數(shù)hxyxyxy)()()(010 ),()()

6、()(000001yxfhyxyhxyxy 1y記為記為)1,., 0(),(1 niyxfhyyiiii定義定義在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計算是精確的前提下，考步計算是精確的前提下，考慮的截斷誤差慮的截斷誤差 Ri = y(xi+1) yi+1 稱為稱為局部截斷誤差局部截斷誤差 /* local truncation error */。定義定義若某算法的局部截斷誤差為若某算法的局部截斷誤差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。 歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：),()()()()()(32112iiiihiiiiiyxhfyh

7、OxyxyhxyyxyR )()(322hOxyih 歐拉法具有歐拉法具有 1 階精度。階精度。Ri 的的主項(xiàng)主項(xiàng)/* leading term */亦稱為亦稱為歐拉折線法歐拉折線法 /* Eulers polygonal arc method*/ 歐拉公式的改進(jìn)：歐拉公式的改進(jìn)： 隱式歐拉法隱式歐拉法 /* implicit Euler method */向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy )1,., 0(),(111 niyxfhyyiiii由于未知數(shù)由于未知數(shù) yi+1 同時出現(xiàn)在等式的兩邊，不能直接得到，故同

8、時出現(xiàn)在等式的兩邊，不能直接得到，故稱為稱為隱式隱式 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。一般先用顯式計算一個初值，再一般先用顯式計算一個初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11)(iiiyxyR)()(322hOxyih 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。 Hey! Isnt the leading term of the local truncation error of Eulers method ? Seems that w

9、e can make a good use of it )(22ihxy 梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的顯、隱式兩種算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的確有局部截斷誤差的確有局部截斷誤差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。階精度，比歐拉方法有了進(jìn)步。但注意到該公式是但注意到該公式是隱式隱式公式，計算時不得不用到公式，計算時不得不用到迭代法，其迭代收斂性與歐拉公式相似。迭代法，其迭代收斂性與歐拉公式相似。)()(311hOyxyRiii 中點(diǎn)歐拉公式中點(diǎn)

10、歐拉公式 /* midpoint formula */中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù)hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假設(shè)假設(shè) ，則可以導(dǎo)出，則可以導(dǎo)出即中點(diǎn)公式具有即中點(diǎn)公式具有 2 階精度。階精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 需要需要2個初值個初值 y0和和 y1來啟動遞推來啟動遞推過程，這樣的算法稱為過程，這樣的算法稱為雙步法雙步法 /* double-step method */，而前面的三種算法都是，而前面的三種算法都是單步法單步

11、法 /* single-step method */。方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點(diǎn)公式中點(diǎn)公式簡單簡單精度低精度低穩(wěn)定性最好穩(wěn)定性最好精度低精度低, 計算量大計算量大精度提高精度提高計算量大計算量大精度提高精度提高, 顯式顯式多一個初值多一個初值, 可能影響精度可能影響精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, lets

12、 make it possible. 改進(jìn)歐拉法改進(jìn)歐拉法 /* modified Eulers method */Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)測預(yù)測，算出，算出),(1iiiiyxfhyy Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：注：此法亦稱為此法亦稱為預(yù)測預(yù)測-校正法校正法 /* predictor-corrector method */。可以證明該算法具有可以證明該算法具有 2 階精度，同時可以看到它是個階精度，同時可以看到它是個單單步步遞推格式

13、，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后面將。后面將看到，它的看到，它的穩(wěn)定性高穩(wěn)定性高于顯式歐拉法。于顯式歐拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii2 龍格龍格 - 庫塔法庫塔法 /* Runge-Kutta Method */建立高精度的單步遞推格式。建立高精度的單步遞推格式。單步遞推法的單步遞推法的基本思想基本思想是從是從 ( xi , yi ) 點(diǎn)出發(fā)，以點(diǎn)出發(fā)，以某一斜某一斜率率沿直線達(dá)到沿直線達(dá)到 ( xi+1 , yi+1 ) 點(diǎn)。歐拉法及其各種變形所點(diǎn)。歐拉法及其各種變形所能達(dá)到的最高精度為能達(dá)

14、到的最高精度為2階階。 考察改進(jìn)的歐拉法，可以將其改寫為：考察改進(jìn)的歐拉法，可以將其改寫為：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的的平均值平均值嗎？嗎？步長一定是一個步長一定是一個h 嗎？嗎？首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點(diǎn)作點(diǎn)作 Taylor 展開展開)(),(),(),(),(2112h

15、OyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將改進(jìn)歐拉法推廣為：將改進(jìn)歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點(diǎn)的點(diǎn)的泰勒泰勒展開作比較展

16、開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hOyxyRiii21,1221 p 這里有這里有 個未知個未知數(shù)，數(shù)， 個方程。個方程。32存在存在無窮多個解無窮多個解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式塔格式。21, 121 p注意到，注意到， 就是改進(jìn)的歐拉法。就是改進(jìn)的歐拉法。 Q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中其中 i ( i = 1, , m )， i (

17、i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為待定均為待定系數(shù)，確定這些系數(shù)的系數(shù)，確定這些系數(shù)的步驟與前面相似。步驟與前面相似。 ).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用為四級最常用為四級4階階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法 /* Classical Runge-Kutta Method */ ：),(),(),(),()22(34222312221432161hKyhxfKKy

18、xfKKyxfKyxfKKKKKyyiihihihihiiihii 注：注： 龍格龍格-庫塔法庫塔法的主要運(yùn)算在于計算的主要運(yùn)算在于計算 Ki 的值，即計算的值，即計算 f 的的值。值。Butcher 于于1965年給出了計算量與可達(dá)到的最高精年給出了計算量與可達(dá)到的最高精度階數(shù)的關(guān)系：度階數(shù)的關(guān)系：753可達(dá)到的最高精度可達(dá)到的最高精度642每步須算每步須算Ki 的個數(shù)的個數(shù))(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好解

19、函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用采用低階算法低階算法而將步長而將步長h 取小取小。HW: p.202 #1, 23 收斂性與穩(wěn)定性收斂性與穩(wěn)定性 /* Convergency and Stability */ 收斂性收斂性 /* Convergency */定義定義 若某算法對于任意固定的若某算法對于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同時同時 i ) 時有時有 yi y( xi )，則稱該算法是，則稱該算法是收斂收斂的。的。 例：例：就初值問題就初值問題 考察歐拉顯式格式的收斂性。考察歐拉顯式格式的收斂性。 0)0(yyyy 解：解：該問題的精確

20、解為該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對任意固定的對任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 穩(wěn)定性穩(wěn)定性 /* Stability */例：例：考察初值問題考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進(jìn)歐拉法改進(jìn)歐拉法 歐歐

21、拉隱式拉隱式歐拉顯式歐拉顯式 節(jié)點(diǎn)節(jié)點(diǎn) xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theore

22、ms are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!定義定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都算中都逐步衰減逐步衰減，則稱該算法是，則稱該算法是絕對穩(wěn)定的絕對穩(wěn)定的 /*absolutely stable */。一般分析時為簡單起見，只考慮一般分析時為簡單起見，只考慮試驗(yàn)方程試驗(yàn)方程 /* test equation */yy 常數(shù)，可以常數(shù)，可以是復(fù)數(shù)是復(fù)數(shù)當(dāng)步長取為當(dāng)步長取為 h

23、 時，將某算法應(yīng)用于上式，并假設(shè)只在初值時，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該，則若此誤差以后逐步衰減，就稱該算法相對于算法相對于 絕對穩(wěn)定絕對穩(wěn)定， 的全體構(gòu)成的全體構(gòu)成絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域。我們稱我們稱算法算法A 比算法比算法B 穩(wěn)定穩(wěn)定，就是指，就是指 A 的絕對穩(wěn)定區(qū)域比的絕對穩(wěn)定區(qū)域比 B 的的大大。000yy h h h例：例：考察顯式歐拉法考察顯式歐拉法011)1(yhyhyyiiii 000yy 011)1(yhyii 01111)1( iiiihyy由此可見，要保證初始誤差由此可見，要保證初始誤差 0 以后逐步衰減，以后逐

24、步衰減，必須滿足：必須滿足：hh 1|1| h0-1-2ReImg例：例：考察隱式歐拉法考察隱式歐拉法11 iiiyhyy iiyhy 11101111 iih可見絕對穩(wěn)定區(qū)域?yàn)椋嚎梢娊^對穩(wěn)定區(qū)域?yàn)椋?|1| h210ReImg注：注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。法的好。例：例：隱式龍格隱式龍格-庫塔法庫塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而顯式顯式 1 4 階方法的絕對穩(wěn)定階方法的絕對穩(wěn)定區(qū)域?yàn)閰^(qū)域?yàn)?)2,2(1111KhyhxfKhKyyiiii其中其中2階

25、方法階方法 的絕對穩(wěn)定區(qū)域?yàn)榈慕^對穩(wěn)定區(qū)域?yàn)?ReImgk=1k=2k=3k=4-1-2-3-123ReImg無條件穩(wěn)定無條件穩(wěn)定HW: p.202 #64 線性多步法線性多步法 /* Multistep Method */用用若干若干節(jié)點(diǎn)處的節(jié)點(diǎn)處的 y 及及 y 值的值的線性組合線性組合來近似來近似y(xi+1)。).(.110111101kikiiikikiiiffffhyyyy 其通式可寫為：其通式可寫為：),(jjjyxff 當(dāng)當(dāng) 1 0 時，為時，為隱式公隱式公式式; 1=0 則為則為顯式公式顯式公式。 基于數(shù)值積分的構(gòu)造法基于數(shù)值積分的構(gòu)造法將將 在在 上積分，得到上積分，得到)

26、,(yxfy ,1iixx 1)(,()()(1iixxiidxxyxfxyxy只要只要近似地算出右邊的積分近似地算出右邊的積分 ，則可通，則可通過過 近似近似y(xi+1) 。而。而選用不同近似式選用不同近似式 Ik，可得到不，可得到不同的計算公式同的計算公式。 1)(,(iixxkdxxyxfIkiiIyy 1 亞當(dāng)姆斯顯式公式亞當(dāng)姆斯顯式公式 /* Adams explicit formulae */利用利用k+1 個節(jié)點(diǎn)上的被積函數(shù)值個節(jié)點(diǎn)上的被積函數(shù)值 構(gòu)造構(gòu)造 k 階牛頓階牛頓后插后插多項(xiàng)式多項(xiàng)式 , 有有kiiifff ,.,11, 0, )( thtxNik 1010)()()

27、(,(1dthhtxRdthhtxNdxxyxfikxxikiiNewton插值余項(xiàng)插值余項(xiàng)dthtxNhyyikii)(101 /* 顯式計算公式顯式計算公式 */局部截斷誤差為：局部截斷誤差為： 1011)()(dthtxRhyxyRikiii例：例：k=1 時有時有)()(11 iiiiiifftfftfhtxN 10111)3(2)(iiiiiiiiffhydtfftfhyydththtdxyfdhRxxi)1(!21)(,(1022 )(1253iyh 注：注：一般有一般有 ，其中，其中Bk 與與yi+1 計算公式計算公式中中 fi , , fi k 各項(xiàng)的各項(xiàng)的系數(shù)系數(shù)均可查表得到

28、均可查表得到 。 )()2(2ikkkiyhBR 10123k2123122324552112162459125243724912583720251fifi 1fi 2fi 3BkMisprint on p.204常用的是常用的是 k = 3 的的4階亞當(dāng)姆斯顯式公式階亞當(dāng)姆斯顯式公式)9375955(243211 iiiiiiffffhyy 亞當(dāng)姆斯隱式公式亞當(dāng)姆斯隱式公式 /* Adams implicit formulae */利用利用k+1 個節(jié)點(diǎn)上的被積函數(shù)值個節(jié)點(diǎn)上的被積函數(shù)值 fi+1 , fi , , fi k+1 構(gòu)造構(gòu)造 k 階階牛頓牛頓前插前插多項(xiàng)式。與顯式多項(xiàng)式完全類似

29、地可得到一系列多項(xiàng)式。與顯式多項(xiàng)式完全類似地可得到一系列隱式公式隱式公式，并有，并有 ，其中，其中 與與 fi+1 , fi , , fi k+1 的系數(shù)亦可查表得到。的系數(shù)亦可查表得到。)()2(2ikkkiyhBR kB10123k21 21125249211282419121 245 241121 241 72019 fi+1fifi 1fi 2Bk常用的是常用的是 k = 3 的的4階亞當(dāng)姆斯隱式公式階亞當(dāng)姆斯隱式公式)5199(242111 iiiiiiffffhyy小于小于Bk較同階顯較同階顯式式穩(wěn)定穩(wěn)定 亞當(dāng)姆斯預(yù)測亞當(dāng)姆斯預(yù)測-校正系統(tǒng)校正系統(tǒng) /* Adams predict

30、or-corrector system */Step 1: 用用Runge-Kutta 法法計算前計算前 k 個初值；個初值；Step 2: 用用Adams 顯式顯式計算計算預(yù)測預(yù)測值；值；Step 3: 用同階用同階Adams 隱式隱式計算計算校正校正值。值。注意：注意：三步所用公三步所用公式的精度必須相同。式的精度必須相同。通常用通常用經(jīng)典經(jīng)典Runge-Kutta 法法配合配合4階階Adams 公式公式。 Hey! Look at the local truncation error of the explicit and implicit Adams methods: and Dont

31、 you think theres something you can do?)(720251)5(5iyh )(72019)5(5iyh 4階階Adams隱式公式的截斷誤差為隱式公式的截斷誤差為)(72019)()5(511iiiyhyxy 4階階Adams顯式公式的截斷誤差為顯式公式的截斷誤差為)(720251)()5(511iiiyhyxy 當(dāng)當(dāng) h 充分小充分小時，可近似認(rèn)為時，可近似認(rèn)為 i i ，則：，則：19251)()(1111 iiiiyxyyxy)(270251)(1111 iiiiyyyxy)(27019)(1111 iiiiyyyxyPredicted value pi

32、+1Modified value mi+1Corrected value ci+1Modified final value yi+1外推技術(shù)外推技術(shù) /* extrapolation */ 基于泰勒展開的構(gòu)造法基于泰勒展開的構(gòu)造法).(.110111101kikiiikikiiiffffhyyyy 將通式中的右端各項(xiàng)將通式中的右端各項(xiàng) yi 1, , yi k ; fi+1, fi 1, , fi k 分別在分別在 xi 點(diǎn)作點(diǎn)作泰勒展開泰勒展開，與精確解，與精確解 y(xi+1) 在在 xi 點(diǎn)的泰勒展開作點(diǎn)的泰勒展開作比較比較。通過令。通過令同同類項(xiàng)系數(shù)相等類項(xiàng)系數(shù)相等，得到足以確定待定系

33、數(shù)，得到足以確定待定系數(shù) 0, , k ; 1, 0, , k 的等式，則可構(gòu)造出線性多步法的的等式，則可構(gòu)造出線性多步法的公式。公式。例：例：設(shè)設(shè))(3322110221101 iiiiiiiiyyyyhyyyy 確定式中待定系數(shù)確定式中待定系數(shù) 0, 1, 2, 0, 1, 2, 3, 使得公式具有使得公式具有4階階精度。精度。解：解：)(5)4(42413612211hOyhyhyhyhyyiiiiii )(225)4(43233422hOyhyhyhyhyyiiiiii )(4)4(3612211hOyhyhyhyyiiiii )(224)4(33422hOyhyhyhyyiiiii

34、)(34)4(3292293hOyhyhyhyyiiiii )()(5)4(42413612211hOyhyhyhyhyxyiiiiii /* y(xi) = yi */1210 hh )2(21021 22132121212)322(hh 36132921212341613)2(hh 424132923416123212414)(hh 個未知數(shù)個未知數(shù)個方程個方程75 令令 1 = 2 = 0Adams 顯式顯式公式公式 以以 y i+1 取代取代 y i 1，并取，并取 1 = 2 = 0Adams 隱式隱式公式公式 以以 yi 3 取代取代 y i 3 ，則可導(dǎo)出另一組，則可導(dǎo)出另一組4

35、 階顯式算法，其中階顯式算法，其中包含了著名的包含了著名的米爾尼米爾尼 /* Milne */ 公式公式)22(342131 iiiiiyyyhyy其局部截斷誤差為其局部截斷誤差為),(,)(45141)5(5 iiiiixxyhR 注：注：上式也可通過上式也可通過數(shù)值積分?jǐn)?shù)值積分導(dǎo)出，即將導(dǎo)出，即將 在區(qū)間在區(qū)間 上積分，得到上積分，得到 再過再過 做做 f 的插值多項(xiàng)式即可。的插值多項(xiàng)式即可。),(yxfy ,13iixx 13,)(,()()(31iixxiidxxyxfxyxy21, iiifff取取 1 = 1, 2 = 0得到得到辛甫生辛甫生 /* Simpson */ 公式公式與

36、與Milne 公式匹配使用公式匹配使用 辛甫生辛甫生 /* Simpson */ 公式公式)4(31111 iiiiiyyyhyy在區(qū)間在區(qū)間xi 1, xi+1上積分，并用上積分，并用Simpson數(shù)值積分?jǐn)?shù)值積分公式來近似積公式來近似積分項(xiàng)，亦可得此分項(xiàng)，亦可得此Simpson公式。公式。 Milne-Simpson 系統(tǒng)的缺點(diǎn)是系統(tǒng)的缺點(diǎn)是穩(wěn)定性差穩(wěn)定性差，為改善穩(wěn)定性，為改善穩(wěn)定性，考慮另一種隱式校正公式：考慮另一種隱式校正公式：)(11011221101iiiiiiiyyyhyyyy 要求公式具有要求公式具有4 階精度。通過泰勒展開，可得到階精度。通過泰勒展開，可得到 個等式，個等式，從中解出從中解出 個未知數(shù)，則有個未知數(shù)，則有 個自由度。個自由度。561取取 1 = 1 得得Simpson 公式公式哈明哈明 /* Hamming */ 用用 1 的不同數(shù)值進(jìn)行試驗(yàn)，發(fā)現(xiàn)當(dāng)?shù)牟煌瑪?shù)值進(jìn)行試驗(yàn)