【文章內(nèi)容簡介】 這些結(jié)果在表 81中，可見計算結(jié)果的精度， Euler法 與 后退 Euler法 差不多，與準(zhǔn)確值相比較 Euler法 偏小，而 后 退 Euler法 偏大； 中點法 與 梯形法 精度同為 2階，但梯形法 更好一些，這跟它們局部截斷誤差的符號，階數(shù)和系 數(shù)的大小是完全一致的。 表見下屏： W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 825 表格 81 表 81 y ?= ?y, y(0)=1 的數(shù)值解 (h=) x 精 確 解 歐拉法 后退歐拉 中點法 梯形法 .1 .904837 .900000 .909091 .900000 .904762 .2 .808731 .810000 .826446 .820220 .818594 .3 .740818 .729000 .751315 .736000 .740633 .4 .670320 .656100 .683013 .627800 .670096 .5 .060531 .590490 .620921 .601440 .606278 .6 .548812 .531441 .654474 .552512 .548537 .7 .496585 . 478298 .5131458 .490938 .496295 .8 .449329 .430467 .466507 .454324 .449029 .9 .406570 .387421 .424098 .400073 .406264 1 .367879 .348679 .385543 .374310 .367573 W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 826 表格 82 而 表 82是分別取了不同的 h= ,h=， h=, h=，還是利用這些公式，經(jīng)過若干步的計算（ h越 小，計算量越大）算到 y(1)的近似值， 可見 ： 隨著 h的減小， y(1)的近似值的精度在提高， ，即 y(1)準(zhǔn)確。 Y ′=y, y(0)=1的解 y(1)的近似值 (y(1)=) h 歐拉法 后退歐拉法 中點法 梯形法 .348678 .385543 .374310 .367573 .366033 .369711 .367944 .367877 .367700 .368052 .367879 .367876 .367800 .367800 .367881 .368020 （緊接下屏） W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 827 表 82計算結(jié)果說明（續(xù)） 但 h太小，到 h=，這與前面所說 h越小， p階越高，應(yīng)該局部截斷誤差越小，因而計算精度更高矛盾了，為什么會產(chǎn)生這種情況呢？ 這是由于 h太小而引起計算量大因而造成了舍入誤差和截斷誤差的積累，這種情況由于初值問題不同可能會影響更大，偏離更嚴(yán)重， 如下面的例 2 。這種問題實際上是穩(wěn)定性問題，我們將會討論方法的穩(wěn)定性，由此得出對 h有一定的要求的穩(wěn)定性制區(qū)域。 W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 828 例 2 求解初值問題 y′=20y， y(0)=1，計算 y(1)的近似值。 解 ： 類似于例 1，用歐拉法、后退歐拉法、中點法、梯形法求解，得到如下 表 。 表 y ?= ?20,y(0)= ?1的解 y(1)的近似值 (y(1)=?8) h 歐拉法 后退歐拉法 中點法 梯形法 .1 1 .169351E4 .514229E+6 0 .01 .203704E9 .120746E7 .413244E+7 .192743E8 .001 .168300E8 .251090E8 .484136E+5 .205979E8 .0001 .202081E8 .210331E8 .475130+3 .206103E8 由 表 ，盡管 中點法 的階數(shù)與 梯形法 相同，比 歐拉法 和 后退歐拉法 的階數(shù)高，計算結(jié)果的精度卻很糟 糕。此外，盡管 歐拉法 與 后退歐拉法 的階數(shù)相同，但 歐 拉法 計算結(jié)果的精度，當(dāng) h= 后退歐拉法 差。 W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 829 167。 2 龍格 庫塔（ Rungekutta）方法 龍格 庫塔方法的基本思想 ： ))(,(*:*,),()(,(1 1 )(8 )(,()()(),()10( )()()(,)()(1111hxyhxfkkxxhxyhxfhxyhxhfxyxyyxfyhxyhxyxyhxyxynnnnnnnnnnnnnnn??????????????????????????????即記作上的平均斜率）稱作區(qū)間這里）得到解于是由微分方程由微分中值定理開始我們從研究差商 因此只要對平均斜率 k*提供一種算法，由（ 811）式 便相應(yīng)地得到一種微分方程的數(shù)值計算公式。 （緊接下屏） W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 830 龍格 庫塔方法的基本思想 (續(xù)） )(21*21 kkk ?? 改進(jìn)歐拉公式比歐拉公式精度高的原因，也就在于確 定平均斜率時，多取了一個點的斜率值。因此它啟發(fā)我們，如果設(shè)法在 [xi,xi+1]上多預(yù)報幾個點的斜率值，然后將它們 加權(quán)平均作為 k*的近似值，則有可能構(gòu)造出更高精度的計 算公式，這是 龍格 庫塔 方法的基本思想。 用這個觀點來研究歐拉公式與改進(jìn)歐拉公式，可以發(fā)現(xiàn) 歐拉公式由于僅取 xn一個點的斜率值 f (xn,yn)作為平均斜率 k* 的近似值，因此精度很低。而改進(jìn)歐拉公式（ 810）卻 是利用了 xn與 xn1兩個點的斜率值 k1 = f (xn,yn)與 k2=f (xn+1,yn+hk1)取算術(shù)平均作為平均斜率 k*的近似值。 其中 k2是通過已知信息 yn利 用歐拉公式求得的。 W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 831 二階龍格 庫塔公式 首先推廣改進(jìn)歐拉公式， 考察區(qū)間 [xn,xn+1]內(nèi)任一點： )10( ????? llhxx nln 我們希望用 xn和 x n+1兩個點的斜率值 k1和 k2加權(quán)平均作為 平均斜率 k*的近似值 ： )(* 221112211 kckchyykckck nn ????? ?即取其中 c1， c2為待定常數(shù)，同 改進(jìn)歐拉公式一樣，這里仍?。? ),(1 nn yxfk ?問題在于怎樣預(yù)測 xn+l處的斜率值 k2。 仿照改進(jìn)歐拉公式，先用 歐拉公式提供 y (xn+l)的預(yù)測值 1l h kyy nln ??? 然后再用預(yù)測值 yn+l通過計算 f 產(chǎn)生斜率值 k2=f (xn+l ,yn+l)， 這樣設(shè)計出的計算公式具有形式 ： W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 832 二階龍格 庫塔公式（續(xù) 1） 1 2 )(8 ),(),()(12122111?????????????l h kyxfkyxfkkckchyynlnnnnn 公式（ 812）中含有三個待定參數(shù) c1,c2和 l，我們希望適 當(dāng)選取這些參數(shù)值，使得公式（ 812）具有二階精度，亦即使： )()( 311 hOyxy nn ?? ??現(xiàn)在仍假定 yn=y(xn)，即 yn是準(zhǔn)確的，將 y(xn+1)與 yn+1都在 xi處作泰勒展開： 1 3 )(8 )()(2)()()( 321 hOxyhxyhxyxy nnnn ????????W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 833 二階龍格 庫塔公式（續(xù) 2） )()()()()],(),(),([),(),()(),(22121hOxylhxyhOyxfyxfyxflhyxfl h kyxfkxyyxfknnnnnnynnxnnnlnnnn????????????????由： 代入（ 812）式，得： 1 4 )(8 )()()()()( 322211 hOxylhcxycchxyy nnnn ????????? 比較（ 813）與（ 814）兩式，要使公式具有二階精度， 只有滿足下列條件： 15)(8 211221????????lccc 這里一共有三個待定參數(shù)，但只需滿足兩個條件，因此有一個自由度，于是滿足條件（ 815）的參數(shù)不止一組，而是一族，相應(yīng)的公式（ 812）也有一族，這些公式統(tǒng)稱為 二 階龍格 庫塔公式 ，簡稱 二階 RK公式 。 W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 834 特別，當(dāng) l=1即 xn+l=xn+1時， c1=c2=1/2，二階 RK公式就 是 改進(jìn)歐拉公式 。 如果取 l=1/2，則 c1=0,c2=1， 這時二階 RK公式稱為 變形的歐拉公式 ，其形式見左邊： 16)(8 2, ),( 1212121????????????????????????khyxfkyxfkhkyynnnnnn 從表面上看，變形的歐拉公式僅含一個斜率值k2，但 k2是 通過 k1計算出來的，因此每完成一步，仍然需要兩次計算函數(shù) f 的值，工作量和改進(jìn)歐拉公式相同。 二階龍格 庫塔公式（續(xù) 3） W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 835 構(gòu)造二階 RK公式的主要步驟 綜上所述，構(gòu)造二階 RK公式主要由以下幾步產(chǎn)生： 1)在區(qū)間 [xn,xn+1]上取二點，預(yù)報相應(yīng)點 的斜率值； 2) 對此兩斜率值加權(quán)平均作為平均斜率 值的近似值； 3) 將 yn+1與 y(xn+1)都在 xi處作泰勒展開， 為使公式達(dá)到二階精度，比較相應(yīng)系 數(shù)，建立有關(guān)參數(shù)所應(yīng)滿足的方程組； 4) 解此方程組得一族二階 RK公式。 W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 836 高階 RK公式 為了進(jìn)一步提高精度，在 [xn,xn+1]上除 xn和 xn+l外再增加 一點 xn+m=xn+mh(l ? m ? 1)，并用 xn,xn+l,xn+m三處的斜率值k1,k2,k3加權(quán)平均作為 k*的近似值，這時計算公式為： )( 3322111 kckckchyy nn ?????其中 k1,k2仍用（ 812）式所取的形式。 為了預(yù)測 xn+m處的斜率值 k3，要定出 xn+m處所對應(yīng)的 yn+m，可以看作在區(qū)間 [xn,xn+m]上使用二階 RK公式，從而 得到 y(xn+m)的預(yù)測值： )( 2211 kbkbmhyy nmn ????于是，再通過計算函數(shù)值 f 得到 : ))(,(),( 22112 kbkbmhymhxfyxfk nnmnmn ????? ??W Y 阜師院數(shù)科院 第八章 常微分方程數(shù)值解法 837 高階 RK公式（續(xù) 1） 這樣設(shè)計出的計算公式具有形式 : 1 7 )(8 ))(,(),(),()(221131213322111???????????