【文章內(nèi)容簡介】 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 7 printf(\n)。 for(i=0。 i=78。 i++) printf (=)。 getchar()。 } 對程序中的諸名稱解釋如下： 符號常量 MAX 表示允許取的離散點的個數(shù)的最大值，用于初始化數(shù)組； 變量 real 表示實際取的離散點的個數(shù)； 數(shù)組 X , Y , F 用于存放各 ix , iy , ( , )iif x y , 0,1, 2,..., 。in? 數(shù)組 CY ,E 用于存放離散 點處的精確解的值 ()iyx 和 iy 的整體誤差ii|y y ( x ) | , i = 0 ,1 ,2 , . .. n。 變量 a,b,x0,y0 用于表示求解區(qū)間的左右端點和給定的初始點； 函數(shù) exact_value() 用于求解問題的精確解； 函數(shù) cal_error() 用于計算各離散點的截斷誤差； 函數(shù) showtable() 用于在屏幕上顯示計算結(jié)果； 第三章 歐拉（ Euler）方法 Euler方法 思想 Euler 方法就是用差分方程初值問題（ ）的解來近似微分方程初值問題（ ）的解，即由公式 （ ） 依次算出 ()nyx 的近似值 ( 1, 2, )nyn?? 。這組公式求問題（ ）的數(shù)值解稱為向前 Euler 公式。 如果在微分方程離散化時，用向后差商代替導數(shù)，即 11 ( ) ( )() nnn y x y xyx h?? ?? ?，則得計算公式 1 1 10( , ) ( 0 , 1 , ) ( 3 . 1 )()n n n ny y h f x y ny y a? ? ?? ? ??? ?? 用這組公式求問題（ ）的數(shù)值解稱為向后 Euler 公式。 向后 Euler 法與 Euler 法形式上相似，但實際計算時卻復雜得多。向前 Euler 公式是顯式的，可直接求解。向后 Euler公式的右端含有 1ny? ，因此是隱式公式，一般要用迭代法求解，迭代公式通常為 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 8 ( 0 )1( 1 ) ( )1 1 1( , ) ( 3 . 2 )( , ) ( 0 , 1 , 2 , )n n n nkkn n n ny y h f x yy y h f x y k? ?? ? ?? ????? ? ??? Euler方法的誤差估計 對于向前 Euler 公式（ ）我們看到，當 1,2,n? 時公式右端的 ny 都是近似的，所以用它計算的 1ny? 會有累積誤差，分析累積誤差比較復雜，這里先討論比較簡單的，所謂局部截斷誤差。 假定用（ ）式時右端的 ny 沒有誤差，即 ()nny y x? ，那么由此算出 1 ( ) ( , ( ) ) ( )n n n ny y x hf x y x? ?? 局部截斷誤差指的是，按（ 3,3）式計算由 nx 到 1nx? 這一步的計算值 1ny? 與精確值1()nyx? 之差 11()nny x y??? 。為了估計它，由 Taylor展開得到的精確值 1()nyx? 是 2 31( ) ( ) ( ) ( ) ( ) ( 3 . 4 )2n n n nhy x y x h y x y x o h? ? ??? ? ? ? （ ）、（ ）兩式相減（注意到 ( , )y f x y?? ）得 2 3211( ) ( ) ( ) ( ) ( 3 . 5 )2n n nhy x y y x o h o h?? ??? ? ? ? 即局部截斷誤差是 2h 階的，而數(shù)值算法的精度定義為： 若一種算法的局部截斷誤差為 1()poh? ，則稱該算法具有 p 階精度。顯然 p 越大，方法的精度越高。式（ ）說明，向前 Euler 方法是一階方法，因此它的精度不高。 改進的 Euler 方法 梯形公式 利用數(shù)值積分方法將微分方程離散化時，若用梯形公式計算式 ()中之右端積分，即 ? ?111( , ( ) ) ( , ( ) ) ( , ( ) )2nnx n n n nx hf x y x d x f x y x f x y x? ????? 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 9 這就是求解初值問題（ ）的梯形公式。 直觀上容易看出，用梯形公式計算數(shù)值積分要比矩形公式好。梯形公式為二階方法。 梯形公式也是隱式格式，一般需用迭代法求解，迭代公式為 ( 0 )1( 1 ) ( )1 1 1( , )( , ) ( , ) ( )2( 0 , 1 , 2 , )n n n nkkn n n n n ny y hf x yhy y f x y f x yk??? ? ?? ??????? ? ?? ?????? 由于函數(shù) ( , )f xy 關于 y滿足 Lipschitz條件，容易看出 ( 1 ) ( ) ( ) ( 1 )1 1 1 12k k k kn n n nhLy y y y??? ? ? ?? ? ? 其中 L為 Lipschitz常數(shù)。因此，當 012hL??時，迭代收斂。但這樣做計算量最大。 如果實際計算時精度要求不太高，用公式（ ）求解時，每步可以只迭代一次，由此導出一種新的方法 — 改進 Euler 法。 改進 Euler 法 按式（ ）計算問題（ ）的數(shù)值解時，如果每步只迭代一次，相當于將 Euler 公式 與梯形公式結(jié)合使用：先用 Euler公式求 1ny? 的一個初步近似值 1ny? ，稱為預測值，然后用梯形公式校正求得近似值 1ny? ，即 ? ?11 1 1()( , ) ( , )2n n n nn n n n n ny y h f x yhy f x y f x y?? ? ?? ? ???? ? ? ???預 算（ ）y 校 正 式（ ）稱為由 Euler 公式和梯形公式得到的 預測 — 校正系統(tǒng)，也叫改進 Euler 法。 為便于編制程序上機，式（ ）常改寫成 1( , )( , ) ( )1 ()2p n n nq n n pn p qy y hf x yy y hf x h yy y y??? ????? ? ???? ???? 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 10 改進 Euler 法是二階方法。 第四 章 休恩方法 休恩方法思想 休恩（ Heun）其實也是歐拉方法的改進型，因此也叫做歐拉改進法。該方法引入一種新的思路，來構(gòu)造求解 [, ]ab 上的初值問題 ( , ( ))y f t y t?? ， 00()yt y? （ ） 要得到解 11[ , ]ty ，可以用微積分基本定理，在 01[ , ]tt 上對 ()yt? 積分得 1100 10( , ( ) ) ( ) ( ) ( )ttf t y t d t y t d t y t y t?? ? ??? （ ） 其中 ()yt? 的不定積分為待定求函數(shù) ()yt 。對 1()yt 求解方程（ ） 1010( ) ( ) ( , ( ) )tty t y t f t y t d t?? ? （ ） 然后可用數(shù)值積分方法逼近（ ）中的定積分，如果采用步長 10h t t?? 的梯形公式，則結(jié)果為 1 0 0 0 1 1( ) ( ) ( ( , ( ) ) ( , ( ) ) )2hy t y t f t y t f t y t? ? ? （ ） 注意公式（ ）的右端包含了待定值 1()yt 。要繼續(xù)求解，需要 1()yt 的一個估計值，歐拉方法的解能夠滿足這一目的，將它代入（ ）后，得到求解 11( , )ty 的公式，稱為休恩（ Heun）方法： 1 0 0 0 1 0 0 0( ) ( ) ( ( , ( ) ) ( , ( , ) ) )2hy t y t f t y t f t y h f t y? ? ? ? （ ） 重復這個過程，得到逼近解曲線 ()y yt? 的一系列點，在每一步中都用歐拉方法作為預報，然 后用梯形公式進行校正，得到最終的值。休恩方法的一般步驟為： 11 ! 1( , ) ,( ( , ) ( , ) )2k k k k k kk k k k k kp y hf t y t t hhy y f t y f t p?? ? ?? ? ? ?? ? ? () 休恩方法的步長和誤差 積分公式（ ）中梯形公式的誤差項為 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 11 3 (2) 3()12kh y c h? （ ） 如果每步中的誤差僅由式（ ）給出，則 M 步后休恩方法的累積誤差將是 3 ( 2 ) 2 ( 2 ) 21 ( ) ( ) ( )1 2 1 2Mkk h b ay c h y c O h? ?? ? ?? （ ） 所以休恩方法的精確度可以這樣計算： 設 ()yt 是初值問題（ ）的解，如果 2 0( ) [ , ]y t C t b? 且 0{( , )}Mk k kty? 是有歐拉方法計算的近似值序列，則 2( ) ( )k k ke y t y O h? ? ? （ ） 311( ) ( , ) ( )k k k k ky t y hf t y O h? ??? ? ? ? （ ） 其中， 1( , ) ( ( , ) ( , ( , ) ) )2k k k k k k k k kht y y f t y f t y h f t y?? ? ? ? ? 特別地，區(qū)間終點處的最終全局誤差滿足， 2( ( ) , ) ( ) ( )ME y b h y b y O h? ? ? （ ） 這就是說，步長如果減小為原來的 1k （ k 為整數(shù)），則可期望最終全局誤差將降至大約 21k。 第五 章 泰勒級數(shù)法 泰勒級數(shù)法有著廣泛的應用 ， 并且是比較求解初值問題的各種不同數(shù)值方法的標準 ， 它可設計為具有任意指定的精度 。 泰勒定理 設 1 0( ) [ , ]Ny t C t b?? ，且 ()yt 在不動點 0[ , ]kt t t b?? 處有 N 次泰勒級數(shù)展開： 1( ) ( ) ( , ( ) ) ( )Nk k N k ky t h y t hT t y t O h ?? ? ? ? （ ） 其中， 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 12 () 11()( , ( ) ) !jN jkN k kjytT t y t hj ??? ? （ ） ( ) ( )( ) ( , ( ))jjy t f t y t? 表示函數(shù) f 關于 t 的 ( 1)j? 次全導數(shù)。求導公式可以遞歸得計算： ()y t f? ? t y t yy f f y f f f?? ?? ? ? ? 22 ( )tt ty y y yy f f y f y f y??? ? ?? ?? ? ? ? 22 ( )tt ty y y y tf f f f f f f f f? ? ? ? ? ( 4 ) 23 3 ( ) 3tt t tt y ty y tyy f f y f y f y? ? ??? ? ? ? 33 ( )y yy yyyf y f y y f y??? ? ?? ?? ? ?