1、數值計算方法數值計算方法第七章第七章常微分方程的數值解法常微分方程的數值解法計算機科學與工程系27.1 引言引言n常微分方程求解常微分方程求解 n微分方程的分類微分方程的分類 n一階常微分方程初值問題的數值解法一階常微分方程初值問題的數值解法 計算機科學與工程系37.1.1 常微分方程求解常微分方程求解n微分方程的定義微分方程的定義n包含自變量、未知函數以及未知函數的導數或微分包含自變量、未知函數以及未知函數的導數或微分的方程叫做微分方程的方程叫做微分方程n在求解微分方程時，必須附加某種定解條件在求解微分方程時，必須附加某種定解條件 n微分方程和定解條件一起組成定解問題微分方程和定解條件一起組

2、成定解問題 n微分方程的解析解法微分方程的解析解法n在數學分析中，求解微分方程采用的方法就是解析在數學分析中，求解微分方程采用的方法就是解析法，求得解析解為未知函數法，求得解析解為未知函數y(x)n如果想求微分方程某點的函數值，只需將某點的點如果想求微分方程某點的函數值，只需將某點的點值值x代入所求的未知函數代入所求的未知函數y(x)，即可求得某點的函數，即可求得某點的函數值值 計算機科學與工程系47.1.1 常微分方程求解常微分方程求解n微分方程的數值解法微分方程的數值解法n在計算方法中求解微分方程采用的方法就是差分法。在計算方法中求解微分方程采用的方法就是差分法。其數值解為某點函數值其數值

3、解為某點函數值y(xi)的近似值的近似值yin即給定一組離散點值即給定一組離散點值x0，x1，x2，xk，就可求得一就可求得一組組 y(x0)，y(x1)，y(x2)，y(xk)，此值為精確值，此值為精確值，但實際求得的是近似值但實際求得的是近似值y0，y1，y2，yk 計算機科學與工程系57.1.2 微分方程的分類微分方程的分類n未知變量的個數未知變量的個數n在微分方程中只含有一個自變量在微分方程中只含有一個自變量x的方程稱作常微的方程稱作常微分方程，相當于一元函數的范疇。其應用廣泛，求分方程，相當于一元函數的范疇。其應用廣泛，求解方法也比較成熟解方法也比較成熟 n在微分方程中，含有多個自變

4、量在微分方程中，含有多個自變量x，y，z，t的方的方程稱為偏微分方程，它相當于多元函數的范疇，其程稱為偏微分方程，它相當于多元函數的范疇，其求解比較困難求解比較困難 00)(),(yxyyxfdxdy計算機科學與工程系67.1.2 微分方程的分類微分方程的分類n例如：對一個各向同性的物體，溫度例如：對一個各向同性的物體，溫度u(x, y, z, t)在在物體內物體內的分布滿足偏微分方程的分布滿足偏微分方程 n求解的條件求解的條件n初值問題：已知函數的導函數和給定的初始條件，初值問題：已知函數的導函數和給定的初始條件，按照假設的步長值按照假設的步長值h，確定一系列的值，確定一系列的值x0，x1，

5、x2，然后求得一系列的函數值，然后求得一系列的函數值y(x0)，y(x1)，y(x2)，的近似值的近似值y0，y1，y2，直到求得所需，直到求得所需點值點值xk處的函數值處的函數值y(xk)的近似值的近似值yk為止為止 ckatuazuyuxu222222221計算機科學與工程系77.1.2 微分方程的分類微分方程的分類n邊值問題：已知函數的導函數和給定的邊界條件，邊值問題：已知函數的導函數和給定的邊界條件，求解滿足微分方程式的某點函數值求解滿足微分方程式的某點函數值n未知函數的個數未知函數的個數n一階常微分方程：只含有一個未知變量一階常微分方程：只含有一個未知變量x，未知函，未知函數的導數階

6、數為一階的，而且只有一個未知函數，數的導數階數為一階的，而且只有一個未知函數，簡言之，方程中只含有簡言之，方程中只含有x、y、y的方程的方程 00)(),(yxyyxfdxdy計算機科學與工程系87.1.2 微分方程的分類微分方程的分類n一階常微分方程組：只含有一個未知變量一階常微分方程組：只含有一個未知變量x，未知，未知函數的導數階數為一階的，而未知函數的個數是多函數的導數階數為一階的，而未知函數的個數是多個，簡言之，方程中含有個，簡言之，方程中含有x，y1，y2，yn，y1，y2，yn的方程的方程 2002100121222111)()(),(),(yxyyxyyyxfyyyxfy計算機科

7、學與工程系97.1.2 微分方程的分類微分方程的分類n未知函數導數的階未知函數導數的階n一階常微分方程：只含有一階常微分方程：只含有x，y，y時，就是最常見時，就是最常見的微分方程的微分方程 n高階常微分方程：除了含有高階常微分方程：除了含有x，y，y以以外，還含有外，還含有未知函數未知函數y的高階導數的高階導數y，y， 0000)()( ) ,( yxyyxyyyxfy計算機科學與工程系107.1.2 微分方程的分類微分方程的分類n高階常微分方程組：凡是由高階常微分方程組：凡是由x，y1，y1， y1 ， ，y2，y2，y2，等組成的微分方，等組成的微分方程程 計算機科學與工程系117.1.

8、3一階常微分方程初值問題的數值解法一階常微分方程初值問題的數值解法 n數值方法的基本思想數值方法的基本思想n在解的存在區間上取在解的存在區間上取n + 1個節點個節點 這里這里 ，i = 0,1, , n，稱為由稱為由xi 到到xi+1的步長的步長 n在這些節點上采用離散化方法，將上述初值問題化在這些節點上采用離散化方法，將上述初值問題化成關于離散變量的相應問題成關于離散變量的相應問題 n把這個相應問題的解把這個相應問題的解yn作為作為y(xn)的近似值的近似值 bxxxxan210iiixxh1計算機科學與工程系127.1.3一階常微分方程初值問題的數值解法一階常微分方程初值問題的數值解法n

9、定理：常微分方程初值問題定理：常微分方程初值問題 設設x0a, b，f(x,y)對對x連續并且關于連續并且關于y滿足李普滿足李普希茲條件：存在常數希茲條件：存在常數L，使，使 對所有對所有xa, b以及任何實數以及任何實數y1，y2都成立，則都成立，則初值問題的解初值問題的解y(x)存在并且唯一存在并且唯一 00)(),(yxyyxfdxdy| ),(),(|2121yyLyxfyxf計算機科學與工程系137.1.3一階常微分方程初值問題的數值解法一階常微分方程初值問題的數值解法n常微分方程初值問題的數值解法常微分方程初值問題的數值解法 n一步法：一步法是在計算一步法：一步法是在計算yn+1時

10、，只用到時，只用到xn+1，xn和和yn，即前一步的值。因此，有了初值以后就可以逐，即前一步的值。因此，有了初值以后就可以逐步往下計算，其代表是龍格庫塔法步往下計算，其代表是龍格庫塔法 n多步法：多步法就是在計算多步法：多步法就是在計算yn+1時，除了用到時，除了用到xn+1，xn和和yn以外，還要用到以外，還要用到xn-p，yn-p (p=1,2,k)，即，即前面前面k步的值，其代表是亞當斯法步的值，其代表是亞當斯法 計算機科學與工程系147.2 歐拉法歐拉法n歐拉法歐拉法n兩步歐拉公式兩步歐拉公式n改進歐拉法改進歐拉法計算機科學與工程系157.2.1.歐拉法歐拉法n公式的導出公式的導出n由

11、于由于y (x0) = y0已給定，因而可以算出已給定，因而可以算出 n設設x1 x0= h充分小，則近似地有充分小，則近似地有 記記n從而我們可以取從而我們可以取),()( 000yxfxy),()( )()(00001yxfxyhxyxy,n, ,ixyyii10 )(),(0001yxhfyy計算機科學與工程系167.2.1.歐拉法歐拉法 作為作為y (x1)的的近似值近似值 n利用利用y1及及f (x1, y1)又可以算出又可以算出y(x2)的近似值的近似值 n一般地，在任意點一般地，在任意點xn+1 = (n + 1)h+x0處處y(x)的近似值的近似值由下式給出由下式給出),(11

12、12yxhfyy),(1nnnnyxhfyy計算機科學與工程系177.2.1.歐拉法歐拉法n歐拉法的幾何意義歐拉法的幾何意義n一階常微分方程一階常微分方程 的解的解 是通過點是通過點 的一條曲線，稱之為微分的一條曲線，稱之為微分方程的積分曲線。積分曲線上每一點方程的積分曲線。積分曲線上每一點 的切線的切線斜率斜率 等于函數等于函數 在這點的值在這點的值 00( , )()yf x yy xy ( )yy x00(,)xy( , )x y( )y x( , )f x y計算機科學與工程系187.2.1.歐拉法歐拉法計算機科學與工程系197.2.1.歐拉法歐拉法n過過 做以做以 為切線斜率的方程為

13、切線斜率的方程 當當 時，得時，得 ，取，取 n過過 做以做以 為切線斜率的方程為切線斜率的方程 當當 時，得時，得 ，取，取 n一般地，過一般地，過 做以做以 為切線斜率的為切線斜率的方程方程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計算機科學與工程系207.2.1.

14、歐拉法歐拉法n當當 時，得時，得 ，取，取n從從 出發逐個算出出發逐個算出 ，對應的數值解，對應的數值解 n一般取一般取 ，得歐拉公式，得歐拉公式n歐拉公式的幾何意義：歐拉公式的幾何意義： 用一條初始點重合的折線，用一條初始點重合的折線，來近似表示微分方程的解（積分曲線）來近似表示微分方程的解（積分曲線） 1nxx11(,)()nnnnnnyyf xyxx1()nny xy0 x12,nx xx12,ny yy1nnxxh1(,)nnnnyyhf xy( )yy x計算機科學與工程系217.2.1.歐拉法歐拉法n歐拉法的其它幾種解釋歐拉法的其它幾種解釋 n歐拉法的數學推導（泰勒展開法）歐拉法的

15、數學推導（泰勒展開法）n將將y(xn+1)在在xn處做二階泰勒展開處做二階泰勒展開n當當h充分小時，忽略高次項得充分小時，忽略高次項得n因此，有歐拉公式因此，有歐拉公式21()()()()()2!nnnnnhy xy xhy xhy xy22()()2!nhyO h1(,)nnnnyyhf xy計算機科學與工程系227.2.1.歐拉法歐拉法n歐拉法的數值積分推導歐拉法的數值積分推導 n將方程將方程 兩端從兩端從xn到到xn+1積分，有積分，有n利用左矩形公式代入算出積分項利用左矩形公式代入算出積分項n離散化，可得離散化，可得 ，有歐拉公式，有歐拉公式( )( , )y xf x y11( )d

16、( , ( )dnnnnxxxxy xxf x y xx11()()( , ( )dnnxnnxy xy xf x y xx1()ny x1(,)nnnnyyhf xy1)(,()(,(nnxxnnxyxhfdxxyxf計算機科學與工程系237.2.1.歐拉法歐拉法n例：用歐拉法求解初值問題例：用歐拉法求解初值問題n解：歐拉公式解：歐拉公式 取步長取步長h=0.1，n=0,1,9時，有時，有1)0(10,2yxyxyy)2(1nnnnnyxyhyy計算機科學與工程系247.2.1.歐拉法歐拉法1 . 1)1021 ( 1 . 01)2(000001yxyhyyn191818. 1)1 . 11

17、 . 021 . 1 ( 1 . 01 . 1)2(111112yxyhyyn計算機科學與工程系257.2.1.歐拉法歐拉法計算機科學與工程系267.2.1.歐拉法歐拉法n例：取步長例：取步長 h=0.2 ，用歐拉法解初值問題，用歐拉法解初值問題 1)0(2yxyyy6 . 0 , 0 x計算機科學與工程系277.2.1.歐拉法歐拉法n解：用歐拉法求解公式，得解：用歐拉法求解公式，得 取步長取步長h=0.2時，時， ，有，有 n=0 )(),(21nnnnnnnnyxyhyyxhfyy6 . 0 , 0 x221000010.2( 1 0 1 )0.8yyhyx y 1n 22211110.8

18、0.20.80.2 0.80.6144yyhyx y2n 計算機科學與工程系287.2.1.歐拉法歐拉法n例例: 用歐拉法求初值問題用歐拉法求初值問題 當當h = 0.02時在區間時在區間0, 0.10上的數值解上的數值解 01)(219 . 000 xxyyxy計算機科學與工程系297.2.1.歐拉法歐拉法n解解: 把把 代入歐拉法計算公式，就得代入歐拉法計算公式，就得 具體計算結果如下表具體計算結果如下表 yxyxf219 . 0),(5, 1, 021018. 01219 . 01nyxyxhyynnnnnn計算機科學與工程系307.2.1.歐拉法歐拉法計算機科學與工程系317.2.1.

19、歐拉法歐拉法n局部截斷誤差和階局部截斷誤差和階n定義：在定義：在yn準確的前提下，即準確的前提下，即yn=y(xn)時，用數值方時，用數值方法計算法計算yn+1的誤差的誤差 稱為該數值方法計算稱為該數值方法計算yn+1時的局部截斷誤差時的局部截斷誤差 n定義：數值方法的局部截斷誤差為定義：數值方法的局部截斷誤差為O(hp+1)，則稱這，則稱這種數值方法的階數為種數值方法的階數為p 11)(nnyxy計算機科學與工程系327.2.1.歐拉法歐拉法n歐拉公式的截斷誤差歐拉公式的截斷誤差n寫出寫出y(xn+1)的泰勒展開式為的泰勒展開式為n由歐拉法得由歐拉法得n兩式相減得兩式相減得 )(! 2)(

20、)()()(21nnnnnxyhxhyxyhxyxy)( )(),(1nnnnnnxhyxyyxhfyy)()(2)(2211hOxyhyxynnn 計算機科學與工程系337.2.1.歐拉法歐拉法n提高計算精度的方法提高計算精度的方法 n增加展開項數增加展開項數n增加增加y(xn+1)函數的泰勒展開式的階數，即逐次取其高次項，函數的泰勒展開式的階數，即逐次取其高次項，這樣就可以減少截斷誤差。其誤差系列為這樣就可以減少截斷誤差。其誤差系列為 O(h2), O(h3), O(h4), n 縮小步長值縮小步長值h n由于差分間隔減小，即由于差分間隔減小，即h值變得很小，可使計算精度提高。值變得很小，

21、可使計算精度提高。在理論上步長值在理論上步長值h越小越好，而且越小越好，而且h趨近于趨近于0最好。但是，最好。但是，h越小會使計算次數增加，浪費計算時間，而且還會使舍越小會使計算次數增加，浪費計算時間，而且還會使舍入誤差增大入誤差增大 計算機科學與工程系347.2.2 兩步歐拉公式兩步歐拉公式n隱式歐拉公式隱式歐拉公式n假設使用向后差商假設使用向后差商 替代方程中的導數替代方程中的導數項項y(xn+1)，再離散化，即可導出下列格式，再離散化，即可導出下列格式n用顯式歐拉公式給出初始值，再由隱式公式進行迭用顯式歐拉公式給出初始值，再由隱式公式進行迭代，因此得到代，因此得到n如果迭代過程收斂，則如

22、果迭代過程收斂，則 為隱式方程的解為隱式方程的解 )()(11nnxyxyh),(111nnnnyxhfyy,.1 , 0),(),()(11)1(1)0(1kyxhfyyyxhfyyknnnknnnnn)(11limknknyy計算機科學與工程系357.2.2 兩步歐拉公式兩步歐拉公式n隱式歐拉公式的局部截斷誤差隱式歐拉公式的局部截斷誤差 211()()2nnnhy xyyx 計算機科學與工程系367.2.2 兩步歐拉公式兩步歐拉公式n兩步歐拉公式兩步歐拉公式 n為了提高精度，改用中心差商為了提高精度，改用中心差商 替代替代 中的導數項，并離散化，有兩步歐拉中的導數項，并離散化，有兩步歐拉公

23、式公式n兩步歐拉公式的局部截斷誤差兩步歐拉公式的局部截斷誤差111 ()()2nny xy xh()(,)nnny xf xy112(,)nnnnyyhf xy231()()()()()( )2!3!nnnnnhhy xy xhy xhy xy xy231()()()()()( )2!3!nnnnnhhy xy xhy xhy xy xy計算機科學與工程系377.2.2 兩步歐拉公式兩步歐拉公式n兩步歐拉公式的局部截斷誤差是兩步歐拉公式的局部截斷誤差是 ，是二階方，是二階方法法 311()()2()( )3nnnhy xy xhy xy311()()2()( )3nnnhy xy xhy xy

24、311()2(,)( )3nnnnhy xyhf xyy112(,)nnnnyyhf xy311()( )3nnhy xyy3()O h計算機科學與工程系387.2.3 改進歐拉法改進歐拉法 n梯形法梯形法 n將方程將方程 的兩端從的兩端從xn到到xn+1求積分得求積分得n用左矩形法進行數值積分則用左矩形法進行數值積分則n因此有因此有),(yxfy 1)(,()()(1nnxxnndxxyxfxyxy1),()(,(nnxxnnyxhfdxxyxf),()(1nnnnyxhfyxy計算機科學與工程系397.2.3 改進歐拉法改進歐拉法n用左矩形法計算右端的積分與用歐拉法計算出的結用左矩形法計算

25、右端的積分與用歐拉法計算出的結果完全相同果完全相同 n用梯形公式計算上式右端的積分，可期望得到較高用梯形公式計算上式右端的積分，可期望得到較高的精度，這時的精度，這時n代入，從而有代入，從而有)(,()(,(2)(,(111nnnnxxxyxfxyxfhdxxyxfnn)(,()(,(2)()(111nnnnnnxyxfxyxfhxyxy計算機科學與工程系407.2.3 改進歐拉法改進歐拉法n將上式中的將上式中的y(xn)，y(xn+1)分別用分別用yn，yn+1代替，作為代替，作為離散化的結果導出下列計算公式離散化的結果導出下列計算公式n將這個結果代入，并將其中的將這個結果代入，并將其中的y

26、(x1)用用y1近似代替，近似代替，則得則得n一般地當求出一般地當求出yn以后，要求以后，要求yn+1，則可歸結為解方程，則可歸結為解方程 ),(),(2111nnnnnnyxfyxfhyy),(),(21110001yxfyxfhyy),(),(2111nnnnnnyxfyxfhyy計算機科學與工程系417.2.3 改進歐拉法改進歐拉法n梯形公式的局部截斷誤差為梯形公式的局部截斷誤差為 )( )( 2)()()(,()(,(2)()()(1111111nnnnnnnnnnnnxyxyhxyxyxyxfxyxfhxyxyyxy)()( 61)( 21)( )()(4321hOxyhxyhxhy

27、xyxynnnnn)()( 21)( )( )( 321hOxyhxhyxyxynnnn)()( 121)()( 21)( )( 22)( 61)( 21)( )(43423211hOxyhhOxyhxhyxyhxyhxyhxhyyxynnnnnnnnn計算機科學與工程系427.2.3 改進歐拉法改進歐拉法n歐拉法與梯形法則相結合，計算公式為歐拉法與梯形法則相結合，計算公式為n先用歐拉法由先用歐拉法由(xn, yn)得出得出y(xn+1)的初始近似值的初始近似值 ，然后進行迭代，反復改進這個近似值，直到然后進行迭代，反復改進這個近似值，直到 （ 為所允許的誤差）為止，并把為所允許的誤差）為止，

28、并把 取作取作y(xn+1)的近似值的近似值yn+1 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn)0(1ny)(1)1(1knknyy)1(1kny計算機科學與工程系437.2.3 改進歐拉法改進歐拉法n改進歐拉法改進歐拉法 n在實際上，當在實際上，當h取值較小時，讓梯形法的迭代公式只取值較小時，讓梯形法的迭代公式只迭代一次就結束迭代一次就結束 n先用歐拉公式求得一個初步近似值先用歐拉公式求得一個初步近似值 ，稱之為預，稱之為預報值報值 n預報值的精度不高，用它替代梯形法右端的預報值的精度不高，用它替代梯形法右端的yn

29、+1，再，再直接計算得出直接計算得出yn+1，并稱之為校正值，這時得到預報，并稱之為校正值，這時得到預報校正公式校正公式 1ny1111(,) (,)(,)2nnnnnnnnnnyyh f xyhyyf xyf xy計算機科學與工程系447.2.3 改進歐拉法改進歐拉法n將預報校正公式稱為改進歐拉公式將預報校正公式稱為改進歐拉公式 。這個公式還。這個公式還可寫為可寫為 n改進歐拉法的程序流程圖改進歐拉法的程序流程圖),(),(2121121211kyhxhfkyxhfkkkyynnnnnn計算機科學與工程系457.2.3 改進歐拉法改進歐拉法計算機科學與工程系467.2.3 改進歐拉法改進歐拉

30、法n改進公式的截斷誤差改進公式的截斷誤差 n由于由于)( )(,(),(1nnnnnxhyxyxhfyxhfk )()( )(,)( )(,)(,)(,)(,)(,)(,),(221112nnnnnnnnnnnnnnnnnnnxyhxhyxyxfyxyxyxfxhxyxhfxyxfykxyxfxhxyxfhkxyhxhfkyhxhfk計算機科學與工程系477.2.3 改進歐拉法改進歐拉法n代入，可得代入，可得ny(xn+1)的二階泰勒展開式為的二階泰勒展開式為n因此有因此有 )(2)( )()(2)( 2)( 2221nnnnnnnnxyhxhyxyxyhxyhxyhyy)()(311hOyx

31、ynn)()( 21)( )()(321hOxyhxhyxyxynnnn計算機科學與工程系487.2.3 改進歐拉法改進歐拉法n例：例： 在區間在區間0, 1.5上，取上，取h = 0.1，求解，求解 1)0(2yyxydxdy計算機科學與工程系497.2.3 改進歐拉法改進歐拉法n解：用歐拉法計算公式如下解：用歐拉法計算公式如下 用迭代一次的改進歐拉法計算公式如下用迭代一次的改進歐拉法計算公式如下 1 . 0, 1201hyyxyhyynnnnnnnnnnyxyhyy2)0(11 . 0, 1222),(),(20)0(11)0(1)0(111hyyxyyxyhyyxfyxfhyynnnnn

32、nnnnnnnn計算機科學與工程系507.2.3 改進歐拉法改進歐拉法計算機科學與工程系517.2.3 改進歐拉法改進歐拉法n例：用歐拉預報例：用歐拉預報-校正公式求解初值問題校正公式求解初值問題 要求步長要求步長 ，計算，計算 的近似值的近似值n 1) 1 (0sin2yxyyy2 . 0h)4 . 1 (),2 . 1 (yy計算機科學與工程系527.2.3 改進歐拉法改進歐拉法n解：由題意知解：由題意知 ，改進歐拉法的具，改進歐拉法的具體形式為體形式為 由由 知知 ，取步長，取步長 ，計算如下，計算如下xyyyxfsin,2cpnnppncnnnnpyyyxyyhyyxyyhyy21si

33、nsin11221) 1 (0 yy10 x2 . 0h計算機科學與工程系537.2.3 改進歐拉法改進歐拉法n=0時時22000022011sin1 0.21 1 sin10.631706sin1 0.20.6317060.631706 sin1.20.79927211 0.6317060.7992720.71548922pcpppcyyhyyxyyhyyxyyy 計算機科學與工程系547.2.3 改進歐拉法改進歐拉法n=1時時21111221222sin0.7154890.20.7154890.715489 sin1.20.476964sin0.7154890.20.4769640.476

34、964 sin1.40.575259110.4769640.5752590.52611222pcpppcyyhyyxyyhyyxyyy715489. 0)2 . 1 (1 yy526112. 0)4 . 1 (2 yy計算機科學與工程系557. 3 龍格庫塔法龍格庫塔法n泰勒級數法泰勒級數法n龍格龍格庫塔法的基本思想庫塔法的基本思想n二階龍格庫塔法和三階龍格庫塔法二階龍格庫塔法和三階龍格庫塔法 n三階龍格庫塔法三階龍格庫塔法 n四階龍格庫塔法四階龍格庫塔法n龍格庫塔法公式推導小結龍格庫塔法公式推導小結 n步長的自動選擇步長的自動選擇 計算機科學與工程系567.3.1 泰勒級數法泰勒級數法n方法

35、實現方法實現n如果初值問題如果初值問題 的精確解的精確解y(x)在在x0, x上存在上存在k + 1階導數且連續，階導數且連續，那么有泰勒公式那么有泰勒公式00)(),(yxyyxfdxdyknkknnnRxykhxhyxyxy)(!)( )()()(1計算機科學與工程系577.3.1 泰勒級數法泰勒級數法n其中截斷誤差為其中截斷誤差為n略去截斷誤差，并用近似值略去截斷誤差，并用近似值 代替真值代替真值 則得則得n用上述公式解初值問題的數值方法稱為泰勒級數法用上述公式解初值問題的數值方法稱為泰勒級數法 11)1(1)()()!1(nnkkkkxxhOykhR)(kny)()(nkxy)(21!

36、 2knknnnnykhyhhyyy 計算機科學與工程系587.3.1 泰勒級數法泰勒級數法n當當 k = 3時，公式變為時，公式變為n這時的截斷誤差是這時的截斷誤差是nnnnnyhyhhyyy 6232114)4(43)()(! 4nnxxhOyhR計算機科學與工程系597.3.1 泰勒級數法泰勒級數法n例：導出用三階泰勒級數法解方程例：導出用三階泰勒級數法解方程 的計算公式的計算公式解：解：22yxy22),(yxyxfy)(222222yxyxyyxfy )3()(242) (22222222yxyxxyyyyfy )32)(8)53(4462422222222)4(yxyxyyxxyy

37、yyyyyyyyyfy 計算機科學與工程系607.3.1 泰勒級數法泰勒級數法n 因此因此n 而而n 其中其中 表示表示f(x, y)對對x的的k階偏導數在階偏導數在x = xn點上點上的值的值 nnnnnfhfhhfyy 32161211)3(43)(! 4nnnxxfhR)(knf計算機科學與工程系617.3.1 泰勒級數法泰勒級數法n泰勒級數法特點泰勒級數法特點n只要初值問題的真解充分光滑，就可以獲得精確度只要初值問題的真解充分光滑，就可以獲得精確度較高的數值解較高的數值解 n須計算須計算y(x)的各階導數，這當的各階導數，這當f (x, y)的表達式復雜時的表達式復雜時是很繁瑣的是很繁

38、瑣的n因此泰勒級數法一般只用于求因此泰勒級數法一般只用于求“表頭表頭”（即開頭幾（即開頭幾點的數值解，如點的數值解，如y1, y2, y3, y4等）等） 計算機科學與工程系627.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想n歐拉法的基本出發點歐拉法的基本出發點 n如果在區間如果在區間xn, xn+1中，只考查一個點中，只考查一個點xn時，則可時，則可得到求解得到求解yn+1的歐拉公式（一點公式）的歐拉公式（一點公式） n如果在區間如果在區間xn, xn+1中，考查兩個點中，考查兩個點xn，xn+1時，則時，則可得到求解可得到求解yn+1的改進歐拉公式（兩點公式）的改進歐拉公式（兩點公式

39、） ),(111nnnnyxhfkkyy),(),(2121121211kyhxhfkyxhfkkkyynnnnnn計算機科學與工程系637.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想n一點公式與兩點公式的比較一點公式與兩點公式的比較 n在區間在區間xn, xn+1內，考查一個點內，考查一個點xn或兩個點或兩個點xn，xn+1 n求解求解yn+1的實質是在的實質是在yn的基礎上增加一個增量的基礎上增加一個增量yn。因此，求解因此，求解yn+1的通式可寫為：的通式可寫為：*1kyynn),(*111nnnnyxhfkkkkyy),(),()(21*121211kyhxhfkyxhfkkkk

40、kyynnnnnn計算機科學與工程系647.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想n如果如果k1=k2，則，則k*=k1；如果；如果k1k2，則，則k*k1；如果；如果k1k1 n考查兩個點的計算公式比只考查一個點的計算公式考查兩個點的計算公式比只考查一個點的計算公式求解精度高求解精度高 計算機科學與工程系657.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想計算機科學與工程系667.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想n預報值的啟示預報值的啟示n在計算新點函數值在計算新點函數值yn+1時，首先計算預報值，然后時，首先計算預報值，然后再修正計算校正值再修正計算校正值 n

41、多預報幾個增量值（或斜率值）多預報幾個增量值（或斜率值）k1, k2, ，然后將，然后將它們加權平均，則有可能得到具有更高精度的計算它們加權平均，則有可能得到具有更高精度的計算公式公式 n初值問題初值問題00)(),(yxyyxfdxdy計算機科學與工程系677.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想 等價于等價于n龍格庫塔法的基本思路就是用龍格庫塔法的基本思路就是用f(x, y)在幾個不同在幾個不同點的數值加權平均代替點的數值加權平均代替 的值，而使的值，而使截斷誤差的階數盡可能高。也就是說，取不同點的截斷誤差的階數盡可能高。也就是說，取不同點的斜率加權平均作為平均斜率，從而提高方

42、法的階數斜率加權平均作為平均斜率，從而提高方法的階數 10),(,()()(,()()(11hxyhxhfxydxxyxfxyxynnnxxnnnn)(,(hxyhxfnn計算機科學與工程系687.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想n常微分方程的離散化表示形式常微分方程的離散化表示形式 是是f(x, y)在一些點的線性組合，解釋為在一些點的線性組合，解釋為f(x, y)的平均斜率，即表示成的平均斜率，即表示成 其中其中),(1hyxhyynnnn),(hyxnniriinnkchyx1),(),(1nnyxfk rikhyhxfkjijijnini,.,3 , 2,11計算機科學

43、與工程系697.3.2 龍格龍格庫塔法的基本思想庫塔法的基本思想n上式中計算了上式中計算了r個個f(x, y)的函數值，稱該式為的函數值，稱該式為r級的龍級的龍格庫塔法格庫塔法 n龍格庫塔法的主要思路龍格庫塔法的主要思路n用用f(x, y)在某些點上的函數值的線性組合，得出在某些點上的函數值的線性組合，得出y(xn+1)的近似值的近似值yn+1。增加調用。增加調用f(x, y)的次數，可以的次數，可以提高精度的階數。或者說，在提高精度的階數。或者說，在xn, xn+1這一步內多這一步內多預報幾個點的斜率值，然后將其加權平均作為平均預報幾個點的斜率值，然后將其加權平均作為平均斜率，則可構造出有更

44、高精度的計算公式斜率，則可構造出有更高精度的計算公式 計算機科學與工程系707.3.3 二階龍格庫塔法二階龍格庫塔法n根據預報校正的思想，首先在區間根據預報校正的思想，首先在區間xn, xn+1內尋內尋找兩個點，即：找兩個點，即：n如果兩點的斜率分別為如果兩點的斜率分別為k1與與k2，則通過線性組合，則通過線性組合，可得到如下的預報校正系統可得到如下的預報校正系統 ) 10(pphxxxxnpnpnn);(),(),()(12122111phkyyphxxyxfkyxfkkkhyynpnnpnpnpnnnnn計算機科學與工程系717.3.3 二階龍格庫塔法二階龍格庫塔法n從公式中可以看出，先計

45、算從公式中可以看出，先計算k1，再計算，再計算yn+p，也就是，也就是先算預報斜率，再算校正斜率先算預報斜率，再算校正斜率 n選擇合適的選擇合適的1，2，p等系數，使得下述公式等系數，使得下述公式 的局部截斷誤差為的局部截斷誤差為O(h3) n將將k1, k2在同一點（在同一點（xn, yn）泰勒展開，則有）泰勒展開，則有11)(nnyxy)( ),(1nnnxyyxfk計算機科學與工程系727.3.3 二階龍格庫塔法二階龍格庫塔法n將將k1，k2的表達式代入的表達式代入)()( )( )(),(),(),(),()(),(),(222),(),(12hOxphyxyhOyxfyxfyxfph

46、yxfhOxyyfxfphyxfphkyphxfknnnnynnnnxnnyxyxnnnnnnnn)()( )( )()()()( )( )( )()(322212222122111hOxyphxyhxyhOxphyxyxyhxyhkkyynnnnnnnnn計算機科學與工程系737.3.3 二階龍格庫塔法二階龍格庫塔法n再對再對y(xn+1)在在xn點進行泰勒展開點進行泰勒展開n令令y n1=y(xn1)，逐項比較，令逐項比較，令h、h2項的系數相等，項的系數相等，便得到系數之間的關系便得到系數之間的關系n當取當取 ， ， ，則有，則有)()( 21)( )()(321hOxyhxhyxyxy

47、nnnn211221p41143232p計算機科學與工程系747.3.3 二階龍格庫塔法二階龍格庫塔法 這種方法稱為休恩公式這種方法稱為休恩公式 n當取當取 ， ， ，此即為改進歐拉公式，此即為改進歐拉公式 )32,(),()3(4113221211hkyxfkyxfkkkhyynnnnnn1p211212),(),(2121121211kyhxhfkyxhfkkkyynnnnnn計算機科學與工程系757.3.3 二階龍格庫塔法二階龍格庫塔法n當取當取 ， 時，時， ，這時二階龍格庫塔，這時二階龍格庫塔公式稱為變形歐拉公式公式稱為變形歐拉公式n截斷誤差為截斷誤差為O(h3) 21p0112)2

48、,(),(1212121khyxfkyxfkhkyynnnnnn計算機科學與工程系767.3.4三階龍格庫塔法三階龍格庫塔法n如果在區間如果在區間xn, xn+1內尋找三個點作為預報校正內尋找三個點作為預報校正系統，則所得到的校正值更為精確系統，則所得到的校正值更為精確 n任意選取三個點為任意選取三個點為 n將上述三個點，分別代入計算斜率的公式，即可先將上述三個點，分別代入計算斜率的公式，即可先后得到斜率值后得到斜率值k1，k2，k3 n再進行線性組合，構成綜合的斜率再進行線性組合，構成綜合的斜率k*，求得較精確，求得較精確的的yn+1的值的值 321110kqpxkpxkxqnpnn對應斜率

49、為對應斜率為對應斜率為計算機科學與工程系777.3.4三階龍格庫塔法三階龍格庫塔法n因此，構造的三階龍格庫塔原始公式為因此，構造的三階龍格庫塔原始公式為 n對對k1，k2，k3分別進行泰勒展開，再將展開式代入分別進行泰勒展開，再將展開式代入yn+1式式 n對對y(xn+1)進行三階泰勒展開進行三階泰勒展開 )(,(),(),()(2131213322111skrkqhyqhxfkphkyphxfkyxfkkkkhyynnnnnnnn計算機科學與工程系787.3.4三階龍格庫塔法三階龍格庫塔法n將將yn+1與與y(xn+1)相比較，可得七個待定系數的函數關相比較，可得七個待定系數的函數關系，得如

50、下三組解的結果系，得如下三組解的結果n三階龍格三階龍格庫塔公式庫塔公式 )2,(2,2),()4(62131213211hkhkyhxfkhkyhxfkyxfkkkkhyynnnnnnnn計算機科學與工程系797.3.4三階龍格庫塔法三階龍格庫塔法n截斷誤差是截斷誤差是O(h4) 計算機科學與工程系807.3.5 四階龍格庫塔法四階龍格庫塔法n四階龍格庫塔法公式四階龍格庫塔法公式n在給定區間在給定區間xn, xn+1內，尋找四個考查點，以便求內，尋找四個考查點，以便求得四個斜率值，即得四個斜率值，即n根據預報校正系統以及加權平均的思想，可得如根據預報校正系統以及加權平均的思想，可得如下的原始四

51、階龍格庫塔法計算公式下的原始四階龍格庫塔法計算公式 43211110krqxkqpxkpxkxrnqnpnn對應斜率為對應斜率為對應斜率為對應斜率為計算機科學與工程系817.3.5 四階龍格庫塔法四階龍格庫塔法)(,()(,(),(),()(3214213121443322111ckbkakrhyrhxfkvkukqhyqhxfkphkyphxfkyxfkkkkkhyynnnnnnnnnn計算機科學與工程系827.3.5 四階龍格庫塔法四階龍格庫塔法n常用的四階公式常用的四階公式n截斷誤差為截斷誤差為O(h5) 342312143211,2,212,21),()22(6hkyhxfkkhyhx

52、fkkhyhxfkyxfkkkkkhyynnnnnnnnnn計算機科學與工程系837.3.5 四階龍格庫塔法四階龍格庫塔法n龍格庫塔四階法的特點龍格庫塔四階法的特點n龍格庫塔四階公式精度較高，可滿足一般工程計龍格庫塔四階公式精度較高，可滿足一般工程計算的要求算的要求 n每次計算每次計算yn+1時，只用到前一步的計算結果時，只用到前一步的計算結果yn，因此，因此在已知在已知y0的條件下，可以自動的進行計算的條件下，可以自動的進行計算 n可以在計算過程中隨時改變步長可以在計算過程中隨時改變步長h n缺點是每前進一步需要多次調用函數缺點是每前進一步需要多次調用函數f(x, y)，工作，工作量較大，并

53、且誤差不容易估計量較大，并且誤差不容易估計 計算機科學與工程系847.3.5 四階龍格庫塔法四階龍格庫塔法計算機科學與工程系857.3.5 四階龍格庫塔法四階龍格庫塔法n例：寫出四階龍格例：寫出四階龍格-庫塔法求解初值問題庫塔法求解初值問題 的計算公式，并取步長的計算公式，并取步長 ，計算，計算 的近的近似值似值 n解：由題意知解：由題意知 ，所以有，所以有 1)0(1yxyy2 . 0h)4 . 0(y1,xyyxf計算機科學與工程系867.3.5 四階龍格庫塔法四階龍格庫塔法11234121132222,16(,)()() 122220.20.2110.90.9122(,)()() 122

54、220.20.20.90.911022nnnnnnnnnnnnnnnnnnnnnnnnhyykkkkkf xyyxhhhhkf xykykxyyxxyxhhhhkf xykykxyyxx 433.910.911,10.20.910.9110.2 10.8180.8181nnnnnnnnnnnnyxkf xh yhkyhkxhyyxxyx 計算機科學與工程系877.3.5 四階龍格庫塔法四階龍格庫塔法 代入，有代入，有 由于由于 ，取，取 ，取步，取步長長 ，有，有 2 . 0181267. 0818733. 022643211nnnnxykkkkhyy1)0(0 yy2 . 0, 0010hx

55、xx2 . 0h018733. 1)2 . 0(1 yy070324. 1)4 . 0(2 yy計算機科學與工程系887.3.5 四階龍格庫塔法四階龍格庫塔法n例：已知一階常微分方程例：已知一階常微分方程 求解：當求解：當x=0.5時，時，y(0.5)=?1 . 00)0( 1 , 01hyxyxy計算機科學與工程系897.3.6 龍格庫塔法公式推導小結龍格庫塔法公式推導小結n求解公式的格式彼此相似求解公式的格式彼此相似 n一階一階n二階二階n三階三階n四階四階11hkyynn2121211kkhyynn6164613211kkkhyynn6162626143211kkkkhyynn計算機科學

56、與工程系907.3.6 龍格庫塔法公式推導小結龍格庫塔法公式推導小結n每個每個k的系數之和為的系數之和為1nk值的求法與求解順序值的求法與求解順序 n在每階的計算式中，斜率在每階的計算式中，斜率k1都相同，即都相同，即 。其它其它ki的則用的則用ki-1，ki-2等求得。求解等求得。求解k的順序為的順序為k1，k2，k3，k4 n積分公式的形式積分公式的形式 n一階龍格庫塔法為矩形積分形式一階龍格庫塔法為矩形積分形式n二階龍格庫塔法為梯形積分形式二階龍格庫塔法為梯形積分形式n三階龍格庫塔法為拋物線形積分形式三階龍格庫塔法為拋物線形積分形式n四階龍格庫塔法為中點修正的拋物線積分形式四階龍格庫塔法

57、為中點修正的拋物線積分形式 ),(1nnyxfk 計算機科學與工程系917.3.6 龍格庫塔法公式推導小結龍格庫塔法公式推導小結n龍格庫塔法的斜率形式為龍格庫塔法的斜率形式為 n龍格庫塔法的積分形式為龍格庫塔法的積分形式為 n計算誤差分析計算誤差分析n一階龍格庫塔法的計算誤差為一階龍格庫塔法的計算誤差為O(h2)n二階龍格庫塔法的計算誤差為二階龍格庫塔法的計算誤差為O(h3)n三階龍格庫塔法的計算誤差為三階龍格庫塔法的計算誤差為O(h4)n四階龍格庫塔法的計算誤差為四階龍格庫塔法的計算誤差為O(h5)*1khyynn1),(1nnxxnnnndxyxfyy計算機科學與工程系927.3.7 步長

58、的自動選擇步長的自動選擇n步長對微分方程數值解的影響步長對微分方程數值解的影響n步長越小，截斷誤差就越小步長越小，截斷誤差就越小 n但是，隨著步長的減小，在一定的求解區間內所需但是，隨著步長的減小，在一定的求解區間內所需的步數就要增多，這樣會引起計算量的增大，并且的步數就要增多，這樣會引起計算量的增大，并且會引起舍入誤差的大量積累與傳播會引起舍入誤差的大量積累與傳播 n變步長求解方法變步長求解方法n從結點從結點xn出發，以出發，以h為步長求一個近似值記為為步長求一個近似值記為yn+1 (h) n由于局部截斷誤差為由于局部截斷誤差為O(h5)，因此有，因此有5)(11)(chyxyhnn計算機科

59、學與工程系937.3.7 步長的自動選擇步長的自動選擇 假定系數假定系數c變化很慢，近似常數，并且在變化很慢，近似常數，并且在h很小時，很小時，c與與h無關無關 n將步長折半，即取將步長折半，即取h/2為步長，從為步長，從xn跨兩步到跨兩步到xn+1，再求得一個近似值再求得一個近似值yn+1(h/2)，每跨一步的截斷誤差是，每跨一步的截斷誤差是O(h/2)5，因此有，因此有n當步長折半后，誤差大約減少當步長折半后，誤差大約減少1/16，即有，即有 5)2(11)2(2)(hcyxyhnn161)()()(11)2(11hnnhnnyxyyxy計算機科學與工程系947.3.7 步長的自動選擇步長

60、的自動選擇n由此得到下列的驗后誤差估計式由此得到下列的驗后誤差估計式n可以通過檢查步長折半前后兩次計算結果的偏差可以通過檢查步長折半前后兩次計算結果的偏差 來判斷所選取的步長是否合適來判斷所選取的步長是否合適 n對于給定的精度對于給定的精度，如果，如果，反復將步長折半進，反復將步長折半進行計算，直至行計算，直至為止，這時取步長折半后的為止，這時取步長折半后的“新新值值”作為結果作為結果 )(151)()(1)2(1)2(11hnhnhnnyyyxy|)(1)2(1hnhnyy計算機科學與工程系957.3.7 步長的自動選擇步長的自動選擇n如果如果為止，這為止，這時取步長加倍前的時取步長加倍前的