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

1、第六章第六章 常微分方程數(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 無(wú)關(guān)的常數(shù)無(wú)關(guān)的常數(shù) L 使使對(duì)任意定義在對(duì)任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則上述都成立

2、，則上述IVP存存在唯一解在唯一解。| ),(),(|2121yyLyxfyxf 要計(jì)算出解函數(shù)要計(jì)算出解函數(shù) y(x) 在一系列節(jié)點(diǎn)在一系列節(jié)點(diǎn) a = x0 x1 xn= b 處的近似值處的近似值),., 1()(nixyyii 節(jié)點(diǎn)間距節(jié)點(diǎn)間距 為步長(zhǎng)，通常采用為步長(zhǎng)，通常采用等距節(jié)點(diǎn)等距節(jié)點(diǎn)，即取即取 hi = h (常數(shù)常數(shù))。) 1,., 0(1 nixxhiii1 歐拉方法歐拉方法 /* Eulers Method */ 歐拉公式：歐拉公式：x0 x1向前差商近似導(dǎo)數(shù)向前差商近似導(dǎo)數(shù)hxyxyxy)()()(010 ),()()()(000001yxfhyxyhxyxy 1y記為

3、記為)1,., 0(),(1 niyxfhyyiiii定義定義在假設(shè)在假設(shè) yi = y(xi)，即第，即第 i 步計(jì)算是精確的前提下，考步計(jì)算是精確的前提下，考慮的截?cái)嗾`差慮的截?cái)嗾`差 Ri = y(xi+1) yi+1 稱為稱為局部截?cái)嗾`差局部截?cái)嗾`差 /* local truncation error */。定義定義若某算法的局部截?cái)嗾`差為若某算法的局部截?cái)嗾`差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。 歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：),()()()()()(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR )()(322hOxyih

4、 歐拉法具有歐拉法具有 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 同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故同時(shí)出現(xiàn)在等式的兩邊，不能直接得到，故稱為稱為隱式隱式

5、 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。一般先用顯式計(jì)算一個(gè)初值，再一般先用顯式計(jì)算一個(gè)初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截?cái)嗾`差：歐拉法的局部截?cái)嗾`差：11)(iiiyxyR)()(322hOxyih 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。 Hey! Isnt the leading term of the local truncation error of Eulers method ? Seems that we can make a good use of i

6、t )(22ihxy 梯形公式梯形公式 /* trapezoid formula */ 顯、隱式兩種算法的顯、隱式兩種算法的平均平均)1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：的確有局部截?cái)嗾`差的確有局部截?cái)嗾`差 ， 即梯形公式具有即梯形公式具有2 階精度，比歐拉方法有了進(jìn)步。階精度，比歐拉方法有了進(jìn)步。但注意到該公式是但注意到該公式是隱式隱式公式，計(jì)算時(shí)不得不用到公式，計(jì)算時(shí)不得不用到迭代法，其迭代收斂性與歐拉公式相似。迭代法，其迭代收斂性與歐拉公式相似。)()(311hOyxyRiii 中點(diǎn)歐拉公式中點(diǎn)歐拉公式 /* midpoint formula *

7、/中心差商近似導(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個(gè)初值個(gè)初值 y0和和 y1來(lái)啟動(dòng)遞推來(lái)啟動(dòng)遞推過程，這樣的算法稱為過程，這樣的算法稱為雙步法雙步法 /* double-step method */，而前面的三種算法都是，而前面的三種算法都是單步法單步法 /* single-step method */

8、。方方 法法 顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點(diǎn)公式中點(diǎn)公式簡(jiǎn)單簡(jiǎn)單精度低精度低穩(wěn)定性最好穩(wěn)定性最好精度低精度低, 計(jì)算量大計(jì)算量大精度提高精度提高計(jì)算量大計(jì)算量大精度提高精度提高, 顯式顯式多一個(gè)初值多一個(gè)初值, 可能影響精度可能影響精度 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 make it possible. 改進(jìn)歐拉法改進(jìn)

9、歐拉法 /* modified Eulers method */Step 1: 先用先用顯式顯式歐拉公式作歐拉公式作預(yù)測(cè)預(yù)測(cè)，算出，算出),(1iiiiyxfhyy Step 2: 再將再將 代入代入隱式隱式梯形公式的右邊作梯形公式的右邊作校正校正，得到，得到1 iy),(),(2111 iiiiiiyxfyxfhyy注：注：此法亦稱為此法亦稱為預(yù)測(cè)預(yù)測(cè)-校正法校正法 /* predictor-corrector method */?？梢宰C明該算法具有可以證明該算法具有 2 階精度，同時(shí)可以看到它是個(gè)階精度，同時(shí)可以看到它是個(gè)單單步步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代

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

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

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

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

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

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

16、低階算法低階算法而將步長(zhǎng)而將步長(zhǎng)h 取小取小。3 收斂性與穩(wěn)定性收斂性與穩(wěn)定性 /* Convergency and Stability */ 收斂性收斂性 /* Convergency */定義定義 若某算法對(duì)于任意固定的若某算法對(duì)于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同時(shí)同時(shí) i ) 時(shí)有時(shí)有 yi y( xi )，則稱該算法是，則稱該算法是收斂收斂的。的。 例：例：就初值問題就初值問題 考察歐拉顯式格式的收斂性?？疾鞖W拉顯式格式的收斂性。 0)0(yyyy 解：解：該問題的精確解為該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy

17、)1 (1 0)1 (yhyii 對(duì)任意固定的對(duì)任意固定的 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)的歐拉格式計(jì)算數(shù)值解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計(jì)算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進(jìn)歐拉法改進(jìn)歐拉法 歐歐拉隱式拉隱式歐拉顯式歐拉顯式 節(jié)點(diǎn)節(jié)點(diǎn) xixey30 1.0000 2.0000

18、 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 theorems are so unstable that a small perturbat

19、ion on the conditions will cause a crash on the conclusions!定義定義若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)若某算法在計(jì)算過程中任一步產(chǎn)生的誤差在以后的計(jì)算中都算中都逐步衰減逐步衰減，則稱該算法是，則稱該算法是絕對(duì)穩(wěn)定的絕對(duì)穩(wěn)定的 /*absolutely stable */。一般分析時(shí)為簡(jiǎn)單起見，只考慮一般分析時(shí)為簡(jiǎn)單起見，只考慮試驗(yàn)方程試驗(yàn)方程 /* test equation */yy 常數(shù)，可以常數(shù)，可以是復(fù)數(shù)是復(fù)數(shù)當(dāng)步長(zhǎng)取為當(dāng)步長(zhǎng)取為 h 時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值時(shí)，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生

20、誤差產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該，則若此誤差以后逐步衰減，就稱該算法相對(duì)于算法相對(duì)于 絕對(duì)穩(wěn)定絕對(duì)穩(wěn)定， 的全體構(gòu)成的全體構(gòu)成絕對(duì)穩(wěn)定區(qū)域絕對(duì)穩(wěn)定區(qū)域。我們稱我們稱算法算法A 比算法比算法B 穩(wěn)定穩(wěn)定，就是指，就是指 A 的絕對(duì)穩(wěn)定區(qū)域比的絕對(duì)穩(wěn)定區(qū)域比 B 的的大大。000yy h h h例：例：考察顯式歐拉法考察顯式歐拉法011)1(yhyhyyiiii 000yy 011)1(yhyii 01111)1( iiiihyy由此可見，要保證初始誤差由此可見，要保證初始誤差 0 以后逐步衰減，以后逐步衰減，必須滿足：必須滿足：hh 1|1| h0-1-2ReImg例：例：考察隱式歐拉法考察隱式歐拉法11 iiiyhyy iiyhy 11101111 iih可見絕對(duì)穩(wěn)定區(qū)域?yàn)椋嚎梢娊^對(duì)穩(wěn)定區(qū)域?yàn)椋?|1| h210ReImg注：注：一般來(lái)說，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式一般來(lái)說，隱式歐拉法的絕對(duì)穩(wěn)定性比同階的顯式法的好。法的好。例：例：隱式龍格隱式龍格-庫(kù)塔法庫(kù)塔法 ),., 1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii 而而顯式顯式 1 4 階方法的絕對(duì)穩(wěn)定階方法的絕對(duì)穩(wěn)