1、常微分方程數(shù)值解法-歐拉法、改進歐拉法和四階龍格庫塔法常微分方程數(shù)值解法常微分方程數(shù)值解法 常微分方程主要有常微分方程主要有: : (1) (1)變量可分離的方程變量可分離的方程 (2)(2)一階線性微分方程一階線性微分方程( (貝努利方程貝努利方程) ) (3) (3)可降階的一類高階方程可降階的一類高階方程 (4)(4)二階常系數(shù)齊次微分方程二階常系數(shù)齊次微分方程 (5)(5)二階常系數(shù)非齊次微分方程二階常系數(shù)非齊次微分方程 (6)(6)全微分方程全微分方程常微分方程數(shù)值解法常微分方程數(shù)值解法l主要內(nèi)容：l一、 引引 言言l二、二、 建立數(shù)值解法的常用方法建立數(shù)值解法的常用方法l三、三、

2、EulerEuler方法方法l四、四、 幾何意義幾何意義l五、五、 EulerEuler方法的誤差估計方法的誤差估計l六、六、 改進歐拉法改進歐拉法l七、四階龍格庫塔法七、四階龍格庫塔法l七、程序設計要求七、程序設計要求主要內(nèi)容主要內(nèi)容 許多實際問題的數(shù)學模型是微分方程或微分許多實際問題的數(shù)學模型是微分方程或微分方程的定解問題方程的定解問題, ,如物體運動如物體運動, ,電路震蕩電路震蕩, ,化學反化學反映及生物群體的變化等映及生物群體的變化等. . 能用解析方法求出精確解的微分方程為數(shù)不能用解析方法求出精確解的微分方程為數(shù)不多多, ,而且有的方程即使有解析解而且有的方程即使有解析解, ,也可

3、能由于解的也可能由于解的表達式非常復雜而不易計算表達式非常復雜而不易計算, ,因此有必要研究微因此有必要研究微分方程的數(shù)值解法分方程的數(shù)值解法 一、引一、引 言言重點 研究一階常微分方程的初值問題的數(shù)值解)1()(),(:00bxayxyyxfdxdy 其一般形式為其一般形式為假定:.),(),(:,),(論知論知這樣由常微分方程的理這樣由常微分方程的理）條件）條件滿足利普希茨（滿足利普希茨（且關于且關于連續(xù)連續(xù)函數(shù)函數(shù)yyLyxfyxfLipschitzyyxf .)(存在并且唯一存在并且唯一解解xyy )1(初值問題初值問題常微分方程數(shù)值解法常微分方程數(shù)值解法bxxxxxaxyNn 210

4、)(在一系列離散節(jié)點在一系列離散節(jié)點求解求解所謂數(shù)值解法，就是尋所謂數(shù)值解法，就是尋初值問題數(shù)值解的提法為定數(shù)，為定數(shù)，定定如不特別說明，總是假如不特別說明，總是假稱為步長。稱為步長。相鄰兩個節(jié)點的間距相鄰兩個節(jié)點的間距上的近似值上的近似值)2 , 1 , 0(,1210ihhxxhyyyyyinnnNn.2 , 1 , 0,0nnhxxn常微分方程數(shù)值解法常微分方程數(shù)值解法建立微分方程數(shù)值解法,首先要將微分方程離散化.一般采用以下幾種方法:(1) 用差商近似導數(shù)11:() ,()nnnnyy xyy x進一步令dxdyhyynn1)(,11,nnnnnnxyxfxxxyxyyxdxdynn

5、二、建立數(shù)值解法的常用方法二、建立數(shù)值解法的常用方法(2) (2) 用數(shù)值積分近似積分用數(shù)值積分近似積分 1)(,()()(1nnxxnndxxyxfxyxy即即)(,)(:11nnnnxyyxyy 令令進一步進一步),()(,(11nnxxnnyxfhdxxyxfyynn 實際上是矩形法實際上是矩形法寬寬高高), 1 , 0(),(11 ndxyxfdxdxdynnnnxxxx常用方法常用方法(3) (3) 用用TaylorTaylor多項式近似并可估計誤差多項式近似并可估計誤差)( ! 2)( )()( ! 2)( )()()(221nnnnnnnxyhxhyxyyhxhyxyhxyxy

6、)(,)(:11nnnnxyyxyy 令令進一步進一步),(1nnnnyxhfyy )( max2)(211xyhyxybxann 常用方法常用方法)1()(),(:00bxayxyyxfdxdy 式為式為已知初值問題的一般形已知初值問題的一般形用差商近似導數(shù)用差商近似導數(shù)dxdyhyynn 1問題轉化為問題轉化為,.)3 , 2 , 1 , 0()(),(001 nxyyyxhfyynnnnEulerEuler方法的迭代公式方法的迭代公式 三、三、EulerEuler方法方法111:(,)(0)1nnnnyyhKKf xyy作等價變換,有,.)3 , 2 , 1 , 0()(),(001 n

7、xyyyxhfyynnnn1(,)nnKf xy令令EulerEuler方法方法 已知，必有切線方程。已知，必有切線方程。及及由于由于斜率斜率的切線（存在?。?，則的切線（存在?。瑒t出發(fā)取解曲線出發(fā)取解曲線由由 , , , , 000000,0000yxyxfyxfyxxy yyxdxdy ),()( 0000,0000yxfxxyyxxxyydxdy ：由點斜式寫出切線方程由點斜式寫出切線方程 四、幾何意義四、幾何意義）（：，可由切線算出，可由切線算出，則，則為為等步長等步長0001101, yxhfyyyhxxh 2 1 0 , )(n1n1n，）（：點的值點的值在在逐步計算出逐步計算出

8、nyxhfyyxxyynn注意：這是“折線法”而非“切線法”除第一個點是曲線切線外，其他點不是！Y=y(x)ab1x2x幾何意義幾何意義五、Euler方法的誤差估計為簡化分析，先考慮計算一步所產(chǎn)生的誤差，即假設為簡化分析，先考慮計算一步所產(chǎn)生的誤差，即假設)(nnxyy 是精確的，是精確的， 估計誤差估計誤差111)(nnnyxyR這種誤差稱為這種誤差稱為局部截斷誤差局部截斷誤差。 估計截斷誤差的主要方法是估計截斷誤差的主要方法是Taylor展開法，即將函數(shù)展開法，即將函數(shù))(xy在在nx處展開：處展開： )(2)()()(2nnnnxyhxyhxyhxy取一次取一次Taylor多項式近似函數(shù)

9、，得多項式近似函數(shù)，得 )()(1hxyxynn)(2)()(2yhxyhxynn )(21)(,()(2yhxyxhfxynnn 12 )(21),( nnnnnxxyhyxhfy得得Euler方法的方法的局部截斷誤差公式局部截斷誤差公式為為 )(21)(2111yhyxyRnnn 結論：結論：上式說明上式說明Euler公式公式的的局部截斷誤差局部截斷誤差為為)(2hO它的精度很差。它的精度很差。 一般很少用它來求近似值，但是一般很少用它來求近似值，但是Euler法法卻體現(xiàn)了數(shù)值方法的基本思想。卻體現(xiàn)了數(shù)值方法的基本思想。 定義定義8.1 如果某種數(shù)值方法的局部截斷誤差為如果某種數(shù)值方法的局

10、部截斷誤差為1()pO h，則稱該方法是，則稱該方法是p階方法或階方法或具有具有p階精度階精度。顯然。顯然p越大，方法的越大，方法的精度越高。精度越高。 注注：Euler方法具有一階精度，因此它的精度不高。方法具有一階精度，因此它的精度不高。 六六改進的改進的EulerEuler方法方法改進的改進的Euler方法方法)1()(),(:00bxayxyyxfdxdy 式為式為已知初值問題的一般形已知初值問題的一般形利用數(shù)值積分將微分方程離散化利用數(shù)值積分將微分方程離散化得梯形公式得梯形公式: :),(),(211 nnnnyxfyxfh 1)(,()()(1nnxxnndxxyxfxyxy),(

11、),(2111 nnnnnnyxfyxfhyy解決方法：有的可化為顯格式，但有的不行解決方法：有的可化為顯格式，但有的不行梯形方法為隱式算法梯形方法為隱式算法改進的改進的Euler方法方法 ,kyxfyxfhyyyxhfyynEulerknnnnnknnnnn210 ,2, 210)(11)1(1)0(1，對，對法結合，形成迭代算法法結合，形成迭代算法與與梯形公式比歐拉法精度高一些梯形公式比歐拉法精度高一些, ,但計算量較大但計算量較大 實際計算中只迭代一次，這樣建立的預測實際計算中只迭代一次，這樣建立的預測校正系統(tǒng)稱作校正系統(tǒng)稱作改進的歐拉公式改進的歐拉公式。改進的改進的Euler方法方法)

12、,(),(2111 nnnnnnyxfyxfhyy 00121211)(),(),()(2:yxyhKyhxfKyxfKKKhyynnnnnn作等價變換作等價變換),(),(2),(1111 kkkkkkkkkkyxfyxfhyyyxhfyy校正校正預測預測改進的改進的Euler方法方法二、改進的Euler法 梯形方法雖然提高了精度，但算法復雜，計算梯形方法雖然提高了精度，但算法復雜，計算量大。如果實際計算時精度要求不太高，用梯形公量大。如果實際計算時精度要求不太高，用梯形公式求解時，每步可以迭代一次，由此導出一種新的式求解時，每步可以迭代一次，由此導出一種新的方法方法改進改進Euler法法。

13、這種方法實際上就是。這種方法實際上就是將將Euler公式與梯形公式結合使用公式與梯形公式結合使用：先用：先用Euler公式公式求求 1ny的一個初步近似值的一個初步近似值 1ny，稱為，稱為預測值預測值，預測值，預測值 的精度可能很差，再用的精度可能很差，再用梯形公式校正梯形公式校正求得近似值求得近似值 1ny ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyy即即 校正預測改進改進Euler法法 亦稱為由亦稱為由Euler公式公式和和梯形公式梯形公式得到的得到的 預測校正預測校正(Predictor-Corrector)系統(tǒng)。系統(tǒng)。 )(21),( ),(1qp

14、npnnqnnnpyyyyhxhfyyyxhfyy為便于上機編程，常改寫成為便于上機編程，常改寫成l改進改進EulerEuler方法計算框圖方法計算框圖開始開始00, ,xy h b輸輸入入1n 1000001112(,)(,)()pcpcpxxhyyhf xyyyhf xyyyy11,xy輸輸出出1xb結結束束01011,nnxxyyY YN N 1)0()10(2)1 . 0(yxyxyyh步步長長求求解解初初值值問問題題例例解解yxyyxf/2),( （1）用）用Euler方法方法得算式為得算式為 1)0()2(1yyxyhyynnnnn（2）用）用改進的改進的Euler方法方法得算式為

15、得算式為 )2(0.1nnnnpyxyyy )1 . 0(2( 1 . 0pnpnqyxyyy )(211qpnyyy 七、七、 龍格龍格 - 庫塔法庫塔法 /* Runge-Kutta Method */建立高精度的單步遞推格式。建立高精度的單步遞推格式。單步遞推法的基本思想是從單步遞推法的基本思想是從 ( xi , yi ) 點出發(fā)，以某一點出發(fā)，以某一斜率沿直線達到斜率沿直線達到 ( xi+1 , yi+1 ) 點。歐拉法及其各種變點。歐拉法及其各種變形所能達到的最高精度為形所能達到的最高精度為2階階。 考察改進的歐拉法，可以將其改寫為：考察改進的歐拉法，可以將其改寫為：),(),(21

16、21121211hKyhxfKyxfKKKhyyiiiiii 斜率斜率一定取一定取K1 K2 的平均值嗎？的平均值嗎？步長一定是一個步長一定是一個h 嗎？嗎？首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設下，使得的前提假設下，使得 )(iixyy )()(311hOyxyRiii Step 1: 將將 K2 在在 ( xi , yi ) 點作點作 Taylor 展開展開)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii )()()(2hOxyphxyii 將

17、改進歐拉法推廣為：將改進歐拉法推廣為：),(),(12122111phKyphxfKyxfKKKhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx Step 2: 將將 K2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiii Step 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點的點的泰勒泰勒展開作比較展開作比較)()()()(322211hOxyphxyhyyiiii )()(2)()()(321hOxyhx

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

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

20、d注：注： 龍格龍格-庫塔法庫塔法的主要運算在于計算的主要運算在于計算 Ki 的值，即計算的值，即計算 f 的的值。值。Butcher 于于1965年給出了計算量與可達到的最高精年給出了計算量與可達到的最高精度階數(shù)的關系：度階數(shù)的關系：753可達到的最高精度可達到的最高精度642每步須算每步須算Ki Ki 的個數(shù)的個數(shù))(2hO)(3hO)(4hO)(5hO)(6hO)(4hO)(2nhO8n 由于龍格由于龍格-庫塔法的導出基于泰勒展開，故精度主要受庫塔法的導出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用低階算法而將步長采用低階算法而將步長h h 取小。取小。lRunge-KuttaRunge-