南大數(shù)值分析課件第九章常微分方程數(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 無關(guān)的常數(shù)無關(guān)的常數(shù) l 使使對任意定義在對任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則上述都成立

2、，則上述ivp存存在唯一解在唯一解。| ),(),(|2121yylyxfyxf 要計算出解函數(shù)要計算出解函數(shù) y(x) 在一系列節(jié)點在一系列節(jié)點 a = x0 x1 xn= b 處的近似值處的近似值),., 1()(nixyyii 節(jié)點間距節(jié)點間距 為步長，通常采用為步長，通常采用等距節(jié)點等距節(jié)點，即取即取 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 步計算是精確的前提下，考步計算是精確的前提下，考慮的截斷誤差慮的截斷誤差 ri = y(xi+1) yi+1 稱為稱為局部截斷誤差局部截斷誤差 /* local truncation error */。定義定義若某算法的局部截斷誤差為若某算法的局部截斷誤差為o(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。 歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：),()()()()()(32112iiiihiiiiiyxhfyhoxyxyhxyyxyr )()(322hoxyih

4、 歐拉法具有歐拉法具有 1 階精度。階精度。ri 的的主項主項/* leading term */亦稱為亦稱為歐拉折線法歐拉折線法 /* eulers polygonal arc method*/ 1 eulers 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)在等式的兩邊，不能直接得到，故同時出現(xiàn)在等式的兩邊，

5、不能直接得到，故稱為稱為隱式隱式 /* implicit */ 歐拉公式，而前者稱為歐拉公式，而前者稱為顯式顯式 /* explicit */ 歐拉公式。歐拉公式。一般先用顯式計算一個初值，再一般先用顯式計算一個初值，再迭代迭代求解。求解。 隱式隱式歐拉法的局部截斷誤差：歐拉法的局部截斷誤差：11)(iiiyxyr)()(322hoxyih 即隱式歐拉公式具有即隱式歐拉公式具有 1 階精度。階精度。 hey! isnt the leading term of the local truncation error of eulers method ? seems that we can make

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

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

8、是單步法單步法 /* single-step method */。方方 法法 1 eulers method顯式歐拉顯式歐拉隱式歐拉隱式歐拉梯形公式梯形公式中點公式中點公式簡單簡單精度低精度低穩(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, ca

9、ll me greedy ok, lets 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 */?？梢宰C明該算法具有可以證明該算法具有 2 階精度，同時可以看到它是個

10、階精度，同時可以看到它是個單單步步遞推格式，比隱式公式的迭代求解過程遞推格式，比隱式公式的迭代求解過程簡單簡單。后面將。后面將看到，它的看到，它的穩(wěn)定性高穩(wěn)定性高于顯式歐拉法。于顯式歐拉法。 )1,., 0(),(,),(211 niyxfhyxfyxfhyyiiiiiiii1 eulers method2 龍格龍格 - 庫塔法庫塔法 /* runge-kutta method */建立高精度的單步遞推格式。建立高精度的單步遞推格式。單步遞推法的單步遞推法的基本思想基本思想是從是從 ( xi , yi ) 點出發(fā)，以點出發(fā)，以某一斜某一斜率率沿直線達(dá)到沿直線達(dá)到 ( xi+1 , yi+1 )

11、 點。歐拉法及其各種變形所點。歐拉法及其各種變形所能達(dá)到的最高精度為能達(dá)到的最高精度為2階階。 考察改進(jìn)的歐拉法，可以將其改寫為：考察改進(jìn)的歐拉法，可以將其改寫為：),(),(2121121211hkyhxfkyxfkkkhyyiiiiii 斜率斜率一定取一定取k1 k2 的的平均值平均值嗎？嗎？步長一定是一個步長一定是一個h 嗎？嗎？2 runge-kutta method首先希望能確定系數(shù)首先希望能確定系數(shù) 1、 2、p，使得到的算法格式有，使得到的算法格式有2階階精度，即在精度，即在 的前提假設(shè)下，使得的前提假設(shè)下，使得 )(iixyy )()(311hoyxyriii step 1:

12、將將 k2 在在 ( xi , yi ) 點作點作 taylor 展開展開)(),(),(),(),(2112hoyxfphkyxphfyxfphkyphxfkiiyiixiiii )()()(2hoxyphxyii 將改進(jìn)歐拉法推廣為：將改進(jìn)歐拉法推廣為：),(),(12122111phkyphxfkyxfkkkhyyiiiiii ),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyx step 2: 將將 k2 代入第代入第1式，得到式，得到 )()()()()()()()(322212211hoxyphxyhyhoxyphxyxyhyyii

13、iiiiii 2 runge-kutta methodstep 3: 將將 yi+1 與與 y( xi+1 ) 在在 xi 點的點的泰勒泰勒展開作比較展開作比較)()()()(322211hoxyphxyhyyiiii )()(2)()()(321hoxyhxyhxyxyiiii 要求要求 ，則必須有：，則必須有：)()(311hoyxyriii21,1221 p 這里有這里有 個未知個未知數(shù)，數(shù)， 個方程。個方程。32存在存在無窮多個解無窮多個解。所有滿足上式的格式統(tǒng)稱為。所有滿足上式的格式統(tǒng)稱為2階龍格階龍格 - 庫庫塔格式塔格式。21, 121 p注意到，注意到， 就是改進(jìn)的歐拉法。就是

14、改進(jìn)的歐拉法。 q: 為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？為獲得更高的精度，應(yīng)該如何進(jìn)一步推廣？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均為待定均為待定系數(shù)，確定這些系數(shù)的系數(shù)，確定這些系數(shù)的步驟與前面相似。步驟與前面相似。 2 runge-kutta method).,(.),(),(),(.1122112321313312122122111 mm mmmmimiiiiiimmiihkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkhyy 最常用為四級最常用為四級4

15、階階經(jīng)典龍格經(jīng)典龍格-庫塔法庫塔法 /* classical runge-kutta method */ ：),(),(),(),()22(34222312221432161hkyhxfkkyxfkkyxfkyxfkkkkkyyiihihihihiiihii 2 runge-kutta method注：注： 龍格龍格-庫塔法庫塔法的主要運算在于計算的主要運算在于計算 ki 的值，即計算的值，即計算 f 的的值。值。butcher 于于1965年給出了計算量與可達(dá)到的最高精年給出了計算量與可達(dá)到的最高精度階數(shù)的關(guān)系：度階數(shù)的關(guān)系：753可達(dá)到的最高精度可達(dá)到的最高精度642每步須算每步須算ki

16、的個數(shù)的個數(shù))(2ho)(3ho)(4ho)(5ho)(6ho)(4ho)(2nho8n 由于龍格由于龍格-庫塔法的導(dǎo)出基于泰勒展開，故精度主要受庫塔法的導(dǎo)出基于泰勒展開，故精度主要受解函數(shù)的光滑性影響。對于光滑性不太好的解，最好解函數(shù)的光滑性影響。對于光滑性不太好的解，最好采用采用低階算法低階算法而將步長而將步長h 取小取小。hw: p.202 #1, 22 runge-kutta methodlab 15. runge-kutta (order four) use runge-kutta method of order four to approximate the solution of

17、 the initial-value problem, , and (1)you are supposed to write a functionvoid rk4 ( double (*f)( ), double a, double b, double y0, int n, file *outfile)to approximate the solution of problem (1) with y= f, x in a, b, and the initial value of y being y0. output the approximating values of y on the n+

18、1 equally spaced grid points from a to b to outfile.inputthere is no input file. instead, you must hand in your function in a *.h file. the rule of naming the *.h file is the same as that of naming the *.c or *.cpp files.output ( represents a space) for each test case, you are supposed to print n+1

19、lines, and each line is in the following format: fprintf( outfile, “%8.4f%12.8en”, x, y );),(yxfy ,bax 0)(yay 2 runge-kutta methodsample judge programsample judge program#include #include #include 98115001_15.h double f0(double x, double y) return (y x*x+1.0); void main( ) file *outfile = fopen(out.

20、txt, w);int n;double a, b, y0; a = 0.0; b = 2.0; y0 = 0.50; n = 10;rk4(f0, a, b, y0, n, outfile);fprintf(outfile, n);fclose(outfile);sample output sample output ( ( represents a space) represents a space) 0.00005.00000000e 0010.20008.29293333e 0010.40001.21407621e+0000.60001.64892202e+0000.80002.127

21、20268e+0001.00002.64082269e+0001.20003.17989417e+0001.40003.73234007e+0001.60004.28340950e+0001.80004.81508569e+0002.00005.30536300e+0003 收斂性與穩(wěn)定性收斂性與穩(wěn)定性 /* convergency and stability */ 收斂性收斂性 /* convergency */定義定義 若某算法對于任意固定的若某算法對于任意固定的 x = xi = x0 + i h，當(dāng)，當(dāng) h0 ( 同時同時 i ) 時有時有 yi y( xi )，則稱該算法是，則稱該算

22、法是收斂收斂的。的。 例：例：就初值問題就初值問題 考察歐拉顯式格式的收斂性。考察歐拉顯式格式的收斂性。 0)0(yyyy 解：解：該問題的精確解為該問題的精確解為 xeyxy 0)( 歐拉公式為歐拉公式為iiiiyhyhyy)1 (1 0)1 (yhyii 對任意固定的對任意固定的 x = xi = i h ，有，有iixhhxihyhyy )1()1(/10/0 ehhh /10)1(lim)(0ixxyeyi 3 convergency and stability 穩(wěn)定性穩(wěn)定性 /* stability */例：例：考察初值問題考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用

23、歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進(jìn)歐拉法改進(jìn)歐拉法 歐拉隱式歐拉隱式歐拉顯式歐拉顯式 節(jié)點節(jié)點 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

24、 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 perturbation on the conditions will cause a crash on the conclusions!3 convergency and stability定義定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都算中都逐步衰減逐步衰減，則稱該算法是，則稱該算法是絕對穩(wěn)定的絕對穩(wěn)定的 /*absolutely stable */。一般分析時為簡單起見，只考