【正文】 1( ) ( ) ( ) ( , ( , ( ) )n n n n n ny x y x h y x h f x y x y x? ?? ? ? 再將 ()nyx 的近似值 ny 代入上式右端，所得結(jié)果作為 1()nyx? 的近似值 1ny? ，得到離散化的計算公式 1 ( , )n n n ny y hf x y? ?? 以上三種方法都是將微分方程離散化的常用方法，每一類方法又可導(dǎo)出不同形式的計算公式。 第二章 數(shù)值解法公共程序模塊分析 編程選擇 : 由于并不需要采用 STL 等泛型程序設(shè)計的方法，采用 C++并不會比采用 C 減少太多代碼，況且這里的實際代碼比較簡單，所以為了減少系統(tǒng)的開銷，采用 Tubro C 來實驗 。 公共程序模塊如下 : 這里為了良好地比較，選用可求解析解的一階常微分方程作為討論 : 2= y , x [ 0,1] y ( 0) = 1 dy xdx y () 其解析式為 y ( x) = 1+2x ,h = /* Filename: */ include include 常微分方程的高精度求解方法 安徽大學(xué)江淮學(xué)院 07 計算機(jī) (1)班 6 include define MAX 100 int real = 1。 /******以下代碼根據(jù)待求解的對象的特殊性進(jìn)行賦值和在 main()中選取 *****/ double func (double x, double y) //計算各離散點處導(dǎo)數(shù)值 { return y *x/y 。 } /********************* changing part end**************************/ void cal_error() //計算誤差值以對各種方法進(jìn)行比較 { int i。 i=real。 E[i] = fabs( CY[i] Y[i])。 } void showtable_s() //微分方程組輸出時用 { //內(nèi)容與 showtable() 類似 } // 輸出各離散點處的 X 值， Y 值， 導(dǎo)數(shù)值， 精 確值， 誤差值 // 分別對應(yīng)于 X[k]， Y[k]， F[k]， CY[k]， E[k] void showtable() //優(yōu)化輸出顯示 { int i,j。 for (i=0。 i++) printf(=)。 printf(%5s%8s%15s%15s%18s%14s,k,X[k],F[k],Y[k],CY[k],E[k])。 i=78。 for(i=0。 i++) { printf(\n%5d,i)。 } 常微分方程的高精度求解方法 安徽大學(xué)江淮學(xué)院 07 計算機(jī) (1)班 7 printf(\n)。 i=78。 getchar()。in? 數(shù)組 CY ,E 用于存放離散 點處的精確解的值 ()iyx 和 iy 的整體誤差ii|y y ( x ) | , i = 0 ,1 ,2 , . .. n。這組公式求問題（ ）的數(shù)值解稱為向前 Euler 公式。 向后 Euler 法與 Euler 法形式上相似，但實際計算時卻復(fù)雜得多。向后 Euler公式的右端含有 1ny? ，因此是隱式公式，一般要用迭代法求解，迭代公式通常為 常微分方程的高精度求解方法 安徽大學(xué)江淮學(xué)院 07 計算機(jī) (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 公式（ ）我們看到，當(dāng) 1,2,n? 時公式右端的 ny 都是近似的，所以用它計算的 1ny? 會有累積誤差，分析累積誤差比較復(fù)雜，這里先討論比較簡單的，所謂局部截斷誤差。為了估計它，由 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 階精度。式（ ）說明，向前 Euler 方法是一階方法，因此它的精度不高。 直觀上容易看出，用梯形公式計算數(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 關(guān)于 y滿足 Lipschitz條件，容易看出 ( 1 ) ( ) ( ) ( 1 )1 1 1 12k k k kn n n nhLy y y y??? ? ? ?? ? ? 其中 L為 Lipschitz常數(shù)。但這樣做計算量最大。 改進(jìn) Euler 法 按式（ ）計算問題（ ）的數(shù)值解時，如果每步只迭代一次，相當(dāng)于將 Euler 公式 與梯形公式結(jié)合使用：先用 Euler公式求 1ny? 的一個初步近似值 1ny? ，稱為預(yù)測值，然后用梯形公式校正求得近似值 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ù) 算（ ）y 校 正 式（ ）稱為由 Euler 公式和梯形公式得到的 預(yù)測 — 校正系統(tǒng)，也叫改進(jìn) Euler 法。 第四 章 休恩方法 休恩方法思想 休恩（ Heun）其實也是歐拉方法的改進(jìn)型，因此也叫做歐拉改進(jìn)法。對 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 。休恩方法的一般步驟為： 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?? ? ?? ? ? ?? ? ? () 休恩方法的步長和誤差 積分公式（ ）中梯形公式的誤差項為 常微分方程的高精度求解方法 安徽大學(xué)江淮學(xué)院 07 計算機(jī) (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? ?? ? ?? （ ） 所以休恩方法的精確度可以這樣計算： 設(shè) ()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è) 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 ?? ? ? ? （ ） 其中， 常微分方程的高精度求解方法 安徽大學(xué)江淮學(xué)院 07 計算機(jī) (1)班 12 () 11()( , ( ) ) !jN jkN k kjytT t y t hj ??? ? （ ） ( ) ( )( ) ( , ( ))jjy t f t y t? 表示函數(shù) f 關(guān)于 t 的 ( 1)j? 次全導(dǎo)數(shù)。 N 次泰勒方法 N 次泰勒方法的一般步驟為： 32 3211 2 ! 3 ! !NNkk d h d hdhy y d h N? ? ? ? ? ? ??? ? （ ） 其中在個步 0,1, , 1kM? ??? ?有 ()()jjkd y t? ， 1,2, ,jN? ??? 。如果 N 固定，則理論上可以推導(dǎo)出步長 h ，使之滿足任意想要的最終全局誤差。 N 次泰勒方法的精度：設(shè) ()yt 為初值問題的解 ,如果 3 0( ) [ , ]y t C t b? ，并且0{( , )}Mk k kty? 是休恩方法產(chǎn)生的一個近似值序列，則 常微分方程的高精度求解方法 安徽大學(xué)江淮學(xué)院 07 計算機(jī) (1)班 13 ( ) ( )Nk k ke y t y O h? ? ? 111( ) ( , ) ( )Nk k k k ky t y hf t y O h? ???? ? ? ? （ ） 特別得，區(qū)間重點處的最終全局誤差滿足， ( ( ) , ) ( ) ( )NME y b h y b y O h? ? ? （ ） 這就是說，步長如果減小為原來的 1k，則可期望最終全局誤差將降至大約 1Nk，k 為整數(shù)。實際上，按照微分中值定理 應(yīng)有 1( ) ( ) ( ) , 0 1nnny x y x y x hh ??? ? ?? ? ? ? 注意到方程 ( , )y f x y?? 就有 1( ) ( ) ( , ( ) ) ( )n n n ny x y x hf x h y x h??? ? ? ? ? 不妨記 ( , ( ) )nnK f x h y x h??? ? ?，稱為區(qū)間 ? ?1,nnxx? 上的平均斜率。 向前 Euler 公式簡單地取 ( , )nnf x y 為 K ，精度自然很低。 如上分析啟示我們，在區(qū)間 ? ?1,nnxx? 內(nèi)多取幾個點，將它們的斜率加權(quán)平均作為 K ，就有可能構(gòu)造出精度更高的計算公式。