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

1、 工程技術和自然科學中的許多實際問題，在數(shù)學上往往可以歸結為求解微分線性方程的形式。高數(shù)中我們學過的解析方法主要是用在一些簡單和特殊的微分方程求解中，而對于大量的一般形式的微分方程，不能用解析方法求出其精確解，而只能用數(shù)值方法求近似解。 在具體求解微分方程時，要有某種定解條件，微分方程與定解條件在一起稱為定解問題。定解條件有兩種：一種是給出積分曲線在初始點的狀態(tài)（初始條件），相應問題稱作初值問題；另一種是給出積分曲線首尾兩端的狀態(tài)（邊界條件），相應的問題稱做邊值問題。 本章主要研究一階常微分方程初值問題的數(shù)值解法。00)(),(yxybxayxfy2121),(),(yyLyxfyxf一階常微

2、分方程的初值問題是指形如的問題。當函數(shù) f (x，y)滿足時初值問題（8-1）的解y = y(x)存在且唯一，此條件稱為李普希茲 (Lipschitz)條件。L為正常數(shù)，稱為Lipschitz常數(shù)。（8-1）我們討論的數(shù)值解法基本思想如下： bxxxxan210 常微分方程的數(shù)值解法期待建立求解y(xi)的近似值yi的遞推格式，通過選取適當?shù)某跏贾担蟮脃(x)在各離散點上的近似值差分方程。在區(qū)間a，b上插入一系列離散點 歐拉法是解初值問題最常用方法，我們用它來說明求解微分方程的基本的技巧和概念。 歐拉法的計算格式是), 2 , 1 , 0()(),(001ixyyyxhfyyiiii(8-2

3、) h=xi+1-xi稱為步長；一般取等步長 h = (b-a)/n。（1）利用y (x)在x = xi處的一階泰勒公 式，得)(! 2)()()(21iiiiyhxyhxyxy 舍去二階導數(shù)以后的余項，并以yi和yi+1分別表示y(xi)和y(xi+1)的近似值，就得到 ),(1iiiiyxhfyy 這個遞推公式可以在已知 yi 時推出yi+1，下面給出公式的幾種解釋：（2）對于充分小的h，將微分方程中的dxdy用差商 iiiixxxyxy11)()(代替，),()()()(11iiiiiiiyxfxyhyyhxyxy解出yi+1即可得到（8-2）式。導數(shù)并以yi和yi+1分別表示y(xi)

4、和y(xi+1)的近似值，就得到 如果采用左矩形公式近似計算上式中的積分，并以yi和yi+1分別表示y(xi)和y(xi+1)的近似值，也可得到),(1iiiiyxhfyy 稱為向前歐拉公式，也稱為顯式歐拉公式。此公式對給定的y0，可以算出y1，然后一步一步算出yi+1。（8-2） （3）從xi到xi+1對式 )(,(xyxfdxdy積分，得 1)(,()()(1iixxiidxxyxfxyxy（8-3） 如果對（8-2）式中的積分采用右矩形公式，將得到計算格式),(111iiiiyxhfyy 稱為向后歐拉公式，也稱為隱式歐拉公式。由于未知量yi+1含在右端，因而不能直接計算，一般使用迭代法進

5、行求解。(8-4) 歐拉法的實質是將連續(xù)問題化為在給定節(jié)點上的差分方程（離散問題）。00)(),(yxyyxfy梯形公式 上述的兩種公式是用矩形求積公式來計算（8-3）式右端的積分而得到的，由于矩形公式精度較低，所以導出的歐拉公式精度也很低。如果用數(shù)值積分的梯形公式計算（8-3）式右端的積分，則有)(,()(,(2)(,(111iiiixxxyxfxyxfhdxxyxfii)(,()(,(2111iiiiiixyxfxyxfhyy（8-5） 這就是梯形公式。梯形公式由于其右端含有未知的yi+1，所以也是隱式方程。 以 f (x0，y0)為斜率過點(x0，y0)做函數(shù)的切線，它與直線x=x1的交

6、點就是y1。依此類推，yi+1是以f (xi，yi)為斜率過點（xi，yi）的切線與x=xi+1的交點，所以歐拉法也稱為歐拉折線法。歐拉方法的幾何意義改進歐拉法 (預估校正法) 顯式歐拉公式計算工作量小，但精度低；梯形公式雖然提高了精度，但為隱式并且計算量較大。為了取長補短，在實際計算中我們將這二種公式結合起來使用。 具體做法是：先用歐拉公式求出一個初步的近似值，記作 1iy，稱為預報值；然后用梯形公式作校正，即用預報值 1iy代替梯形公式（8-5）右邊的yi+1，再直接計算，得到的值yi+1是校正值，這樣建立起來的預估校正公式稱作為改進歐拉公式，即預報 ),(1iiiiyxhfyy校正 寫成

7、迭代形式),(),(2111iiiiiiyxfyxfhyy), 2 , 1 , 0)(,(),(2),()(1) 1(1)0(1kyhxfyxfhyyyxhfyykiiiiikiiiii或表示為（編程時使用）)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy歐拉公式的截斷誤差與精度局部截斷誤差：對于數(shù)值方法，若假定在前n個節(jié)點處的數(shù)值解等于其精確解，即yi= y (xi) (n=0,1,2，,n-1)則 Ri+1 = y(xi+1) -yi+1表示為數(shù)值方法在節(jié)點xi+1處的局部截斷誤差。局部截斷誤差在一定程度上反映了方法的精度。常用泰勒展式來討論局部截斷誤差。 下面

8、討論歐拉法（以向前歐拉公式為例）的局部截斷誤差，將泰勒展開式 )(! 2)()()(21iiiixyhxyhxyxy與歐拉公式 ),(1iiiiyxhfyy對照可得 )()(2)(22111hOxyhyxyRiiii 定義定義 若一種數(shù)值方法的局部截斷誤差為)(11PihOR，則稱該方法具有P 階精度或P 階方法。 推導得歐拉公式具有一階精度（一階方法）；梯形公式有二階精度（二階方法）。 例1 用向前歐拉公式和改進歐拉法公式求解初值問題并比較它們的精度1) 0(102yxyxyy取h=0.1；方程的準確解是 xy21解：（1）用歐拉法計算顯然 yxyyxf2),(x0=a=0，n=10，b=1

9、，y0=1由向前歐拉公式),(111iiiiyxhfyy1 . 1)1021 ( 1 . 01)2(00001yxyhyy1918. 1)1 . 11 . 021 . 1 ( 1 . 01 . 1)2(11112yxyhyy（2）用改進歐拉法公式計算)9 , 2 , 1 , 0( 1 . 0, 1)22(2)2(011111ihyyxyyxyhyyyxyhyyiiiiiiiiiiiii以此求出各結點上的數(shù)值解，計算結果見下12 iixyxi歐拉公式歐拉公式y(tǒng)i預估預估校正公式校正公式y(tǒng)i精確值精確值0.01110.11.11.0959091.0954450.21.1918181.1840971

10、.1832160.31.2771381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859661.4832400.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.732051 同精確解比較，第同精確解比較，第2列歐拉公式的結果大約只有兩位有效數(shù)字，列歐拉公式的結果大約只有兩位有效數(shù)字，而用預估而用預估校正公式的結果則有三位有效

11、數(shù)字。校正公式的結果則有三位有效數(shù)字。 00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前 進 Eul er公 式改 進 Eul er公 式精 確值例2 用歐拉法和改進歐拉法公式解初值問題： 4 . 00 . 0 ,1)0(2xyydxdy解：取 f (x) = y2，y0 = 1，h = 0.1（1）用歐拉法求解 ),(1iiiiyxhfyy1 . 111 . 0121y221. 11 . 11 . 01 . 122y3701. 1221. 11 . 0221. 123y5578. 13701. 11 . 03701. 124y該方

12、程的精確解是 xy11，計算結果如下表。 ixiyiy(xi)| |y(xi)-yi| |10.11.11.11110.011120.21.2211.25000.029030.31.37011.42360.053540.41.55781.66670.1089（2）用改進歐拉法公式求解 )(22)(12)1(12)0(1kiiikiiiiyyhyyhyyy用以上的迭代公式計算，結果如下表： ixiyiy(xi)| |y(xi)-yi| |10.11.11181.11110.000720.21.25201.25000.002030.31.43311.42360.009540.41.67631.66

13、670.0004 兩種方法進行比較，發(fā)現(xiàn)用預估校正公式誤差較小。 在上一節(jié)中，我們得到了一些求微分方程近似解的數(shù)值方法，這些方法的局部截斷誤差較大，精度較低，我們希望得到有更高階精度的方法。 考察差商 ，由微分中值定理，存在點 ，使得hxyxyii)()(1),(),()()(11iiiixxyhxyxy*)()(,()()(1hkxyyhfxyxyiii)(,(*yfk 其中 稱為y(x)在xi，xi+1上的平均斜率。由此考慮，只要對k*提供一種近似算法，便可導出一種計算格式，這就是龍格庫塔方法的基本思想。（8-6）即:ix1ixxy)(xyy *k 如果以y(x)在xi處的斜率作為y(x)

14、在 xi，xi+1上的平均斜率k*，即),()(,)(*iiiiiyxfxyxfxykix1ixxy)(xyy k*這就是歐拉法則得 ),(1iiiiyxhfyy 在xi，xi+1上取兩點xi，xi+p（0 p1）的斜率值k1，k2的線性組合1k1+2k2作為 k*的近似值(1、2為待定常數(shù))，此公式 一般形式可寫成)(22111kkhyyii其中k1 = f (xi，yi)，k2為xi，xi+1內任意一點xi+p = xi+ ph (0 p1) 的斜率f (xi+p，y(xi+p)。 由于y (xi+p)并沒有給出，所以首先應該求y (xi+p)，仿照改進歐拉公式的構造思想，得到下面的公式)

15、,(),()(12122111phkyphxfkyxfkkkhyyiiiiii（8-7）1),(phkyyxphfyyiiiipi這樣構造出的公式為ix1ixxy)(xyy k1k2k* 公式中含有三個參數(shù)1，2和p，如果我們適當選取參數(shù)的值，可以使公式的局部截斷誤差為O(h3)。對k1和k2作泰勒展開iiiyyxfk),(1),(),(112phkyphxfphkyxfkiiipi 代入（8-7）得)()()()()(),(),(),(2221hOyphyhOxyphxyhOyxfphkyxfphyxfiiiiiiyiixii )()()()()(322211hOxyphxyhxyyiiii

16、 （*） 又 y (x)在xi處的二階泰勒展開式為 32)(! 3)()(! 2)()()()(iiiiiixxyxxxyxxxyxyxy 當x = xi+1時， )()(! 3)(33hOxxyi ，有 )()(! 2)()()(321hOxyhxyhxyxyiiii （*） （*）)()(32221hOyphyhyiii )(! 232hOyhyhyiii 比較（*）與（*）的系數(shù)即可發(fā)現(xiàn)，要使公式（7-7）的局部截斷誤差滿足 )()(311hOxyyii，即要求公式具有二階精度只要下列條件成立即可。 211221p（8-8） 滿足條件（8-8）的一簇公式統(tǒng)稱為二階龍格庫塔公式。特別的，當

17、2121, 1p塔公式就成為改進歐拉公式。時，龍格-庫ix1ixxy)(xyy *k1k2k改進歐拉公式就是以y(x)在xi和xi+1處的斜率k1和k2的算術平均值作為y(x)在xi，xi+1上的平均斜率k*來進行計算的。1, 0,2121p 若取 時，龍格-庫塔公式就稱為變形的歐拉公式，其形式為),(),(1222121kyxfkyxfkhkyyhihiiiii（8-9）12khyi21ix此處的就是歐拉方法預報出的中點處的近似解；而 )2,2(12khyhxfkii等于中點的斜率值 則近似)(,(2121iixyxf； 所以公式可以看作用中點斜率近似代替平均斜率k*，因此，公式（8-9）也

18、稱作中點公式。 2kix1ixxy)(xyy 2/hix*k 為了提高精度，考慮在xi，xi+1上 取三點xi，xi+p，xi+q的斜率值分別為k1，k2，k3，將k1，k2，k3的線性組合作為平均斜率k*的近似值，其中xi+q=xi+qh （0pq1）這時的計算公式為)()(3322111kkkhyxyii),(1iiyxfk ),(12phkyphxfkii其中 為了預報點xi+q的斜率k3，一種很自然的想法是用歐拉法預報，即取),(),(13phkyqhxfyxfkiiqiqi)(21skrkqhyyiqi 但是，這樣做效率比較差。因為在區(qū)間xi，xi+q內已有兩個斜率可以使用，所以把k

19、1，k2的線性組合作為xi，xi+q上平均斜率的近似值，當然比用歐拉法預報y(x i+q)精度要好。由此，得到y(tǒng)(xi+q)的預報值。)(,(),(213skrkqhyxfyxfknqiqiqi)(,(),(),()(2131213322111skrkqhyqhxfkphkyphxfkyxfkkkkhyyiiiiiiii于是可取從而得到三點計算公式的形式（8-10）下面列出其中的一種形式)2,()2,2(),()4(62131213211hkhkyhxfkkhyhxfkyxfkkkkhyyiiiiiiii 類似前面的推導，利用泰勒展開適當選擇參數(shù)1，2和p，q，r，s的值使上述公式具有三階精度

20、，即局部截斷誤差為O(h4)。這類公式統(tǒng)稱為三階龍格-庫塔公式。（8-11）ix1ix2hxixy)(xyy *k1k2k3k高階龍格-庫塔法 為了進一步提高精度，在xi，xi+1上可取多個點，預報相應點的斜率值，對這些斜率值加權平均作為平均斜率值。利用泰勒展開，比較相應系數(shù)，從而確定在盡可能高的精度下有關參數(shù)應滿足的條件。 從理論上講，可以構造任意高階的龍格-庫塔公式，但實踐證明，高于四階的龍格-庫塔公式，不但計算量大，而且精確度并不一定提高。在實際計算中，四階龍格-庫塔公式是精度及計算量較理想的公式。經(jīng)典的四階龍格-庫塔公式為)22(643211kkkkhyyii),(1iiyxfk )2

21、,2(12khyhxfkii),(34hkyhxfkii)2,2(23khyhxfkii（8-12）龍格-庫塔法的步長選擇 上面所討論的龍格-庫塔方法是定步長的，單從每一步來說，步長 h 越小，局部截斷誤差就越?。坏S著步長的縮小，不但引起計算量的增加，而且也有可能引起舍入誤差的嚴重積累；而步長h太大又不能達到預期的精度要求，所以怎樣選擇合適的步長h，這在實際計算中是很重要的。 我們以四階龍格-庫塔公式為例，從節(jié)點xi出發(fā)，先以某個h為步長求出一個近似值，記為yi+1(h)，然后將步長折半，即取h/2為步長從xi計算兩步到xi+1，再求得一個近似值yi+1(h/2) 。 通過計算| yi+1(

22、h/2) -yi+1(h) |是否成立，成立則步長選擇h/2，否則繼續(xù)將區(qū)間分成h/4，再次重復上述步驟至滿足精度要求。例 使用高階R-K方法計算初值問題1)0(5 . 002yxyy. 1 . 0h取解： (1) 使用三階R-K方法時1i201yk 12102)21 . 0(kyk1025. 121203)2( 1 . 0(kkyk2555. 1)4(61 . 032101kkkyy1111. 1)2,()2,2(),()4(62131213211hkhkyhxfkkhyhxfkyxfkkkkhyyiiiiiiii其余結果如下:(2) 如果使用四階R-K方法時1i201yk 12102)21

23、 . 0(kyk1025. 1 i xi k1 k2 k3 yi 1.0000 0.1000 1.0000 1.1025 1.2555 1.1111 2.0000 0.2000 1.2345 1.3755 1.5945 1.2499 3.0000 0.3000 1.5624 1.7637 2.0922 1.4284 4.0000 0.4000 2.0404 2.3423 2.8658 1.6664 5.0000 0.5000 2.7768 3.2587 4.1634 1.9993)22(643211kkkkhyyii),(1iiyxfk )2,2(12khyhxfkii),(34hkyhxfk

24、ii)2,2(23khyhxfkii2203)21 . 0(kyk1133. 12304)1 . 0(kyk2351. 1)22(61 . 0432101kkkkyy1111. 1其余結果如下: i xi k1 k2 k3 k4 yi 1.0000 0.1000 1.0000 1.1025 1.1133 1.2351 1.1111 2.0000 0.2000 1.2346 1.3756 1.3921 1.5633 1.2500 3.0000 0.3000 1.5625 1.7639 1.7908 2.0423 1.4286 4.0000 0.4000 2.0408 2.3428 2.3892 2.7805 1.6667 5.0000 0.5000 2.7777 3.2600 3.3476 4.0057 2.0000)22(643211