【正文】 ......................................... 26 主要參考文獻 ................................................. 26 致 謝 ..................................................... 27 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 1 第一章 前 言 案例引入微分方程概念 在科技、工程、經濟管理、生態(tài)、生態(tài)、刑偵等各個領域微分方程有著廣泛的應用。同樣，一塊冷的物體，其溫度上升的速度是與他自身溫度同外界溫度的差值成正比。 微分方程的基本概念 微分方程及微分方程的階 含未知函數的導數 (或微分 )的方程稱為微分方程； 未知函數是一元函數的微分方程，稱為常微分方程 ； 未知函數是多元函數的 微分方程，稱為偏微分方 ； 2d 3 d ( )y x x? ， 2 2d ( )d s gt ? ()和 ()式均是微分方程 . 微分方程中未知函數的導數的最高階數，稱為微分方程的階 . 微分方程 ()是一階的，微分方程 ()是二階的 . 微分方程的解、通解與特解 能使微分方程成為恒等式的函數，稱為微分方程的解 . 例如 3y x c??和 3y x 1??都是 3dy 3xdx? 的解 . 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 2 又如 21 212 CtCgts ???和 212s gt?都是 22dds gt ?的解 . 如果微分方程的解中含任意常數 ,且獨立的 (即不可合并而使個數減少的 )任意常數的個數與微分方程的階數相同，這樣的解為微分方程的通解 . 不包含任意常數的解為微分方程特解 . 微分方程的初值條件及其提法 用以確定微分方程解中任意常數的特定條件，稱為微分方程的初值條件 . 初值條件的提法 ： 當 x=x0 時， y=y0， 微分方程的解的幾何意義 . 微分方程的解的圖形稱為微分方程的積分曲線 .通解的圖形是一族積分曲線，稱為微分方程的積分曲線族 .微分方程的某個特解的圖形就是積分曲線族中滿足給定初值條件的某一特定的積分曲線 . 221 2 1 200e e ( ) 4 0 ( ) 0 , 1 ( xxxxy C C C Cy yy y 39。y39。yxy nn ， 個初值條件：階微分方程需給出一般地，對于 ?常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 3 ，得 ，分別求一階及二階導數將函數 xxxxCCy39。 數值求解微分方程的方法基于有限維近似，這個過程稱為離散化，我們將用代數方程代替微分方程，用代數方程的解近似微分方程的解，對初值問題來說，近似解的值是在求解區(qū)間上一步步地產生的，因此求解常微分方程的數值方法也稱為離散變量法，在由一個離散點的值計算下一個點的值時，一般會產生一定的誤差，這樣新的近似解將落在常微分方程的另一個解上，而這個解與開始所求的解是不同的，解的穩(wěn)定 性決定了這類誤差將隨時間的增大而放大或縮小。 建立數值解法，首先要將微分方程離散化，一般采用以下幾種方法： （ i）用差商近似導數 若用向前差商 1( ) ( )nny x y xh? ? 代替 ()nyx? 代入（ ）中的微分方程，則得 1( ) ( ) ( , ( ) ) ( 0 , 1 , )nnnny x y x f x y x nh? ? ?? 化簡得 1( ) ( ) ( , ( ) )n n n ny x y x h f x y x? ?? 如果用 ()nyx 的近似值 ny 代入上式右端，所得結果作為 1()nyx? 的近似值，記為1ny? ， 則有 1 ( , ) ( 0 , 1 , ) ( )n n n ny y hf x y n? ? ? ? 這樣，問題（ ）的近似值可通過求解下述問題 10( , ) ( 0 , 1 , ) ( 1 . 7 )()n n n ny y h f x y ny y a? ? ? ??? ?? 得到，按式（ ）由初值 0y 可逐次算出 12,yy 。例如，對微分方程兩端積分，得 11( ) ( ) ( , ( ) ) ( 0 , 1 , ) ( 1 . 8 )nnxnn xy x y x f x y x d x n?? ? ? ?? 右邊的積分用矩形公式或梯形公式計算。 編程風格 : 按照常微分方程數值解三個基本步驟 ： 將問題離散 化；建立遞推格式；按步進法計算， 所以求微分方程的數值解的算法框架都是 相同的，不同 的是所使用的遞推形式不同，則 可以用公共子程序來代替，對不同的方法的計算結果用統一的格式來顯示，同時也可以比較不同方法的精確度 [4]。 } double exact_value (double x) //計算各離散點解析解以測數值解精度 { return sqrt( *x + ) 。 i++) { CY[i] = exact_value( X[i] )。 printf (\n)。 printf(\n)。 i++) printf()。 printf(%%%%%,X[i],F[i], Y[i],CY[i],E[i])。 i++) printf (=)。 變量 a,b,x0,y0 用于表示求解區(qū)間的左右端點和給定的初始點； 函數 exact_value() 用于求解問題的精確解； 函數 cal_error() 用于計算各離散點的截斷誤差； 函數 showtable() 用于在屏幕上顯示計算結果； 第三章 歐拉（ Euler）方法 Euler方法 思想 Euler 方法就是用差分方程初值問題（ ）的解來近似微分方程初值問題（ ）的解，即由公式 （ ） 依次算出 ()nyx 的近似值 ( 1, 2, )nyn?? 。向前 Euler 公式是顯式的，可直接求解。顯然 p 越大，方法的精度越高。梯形公式為二階方法。 如果實際計算時精度要求不太高，用公式（ ）求解時，每步可以只迭代一次，由此導出一種新的方法 — 改進 Euler 法。該方法引入一種新的思路，來構造求解 [, ]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? 的不定積分為待定求函數 ()yt 。 第五 章 泰勒級數法 泰勒級數法有著廣泛的應用 ， 并且是比較求解初值問題的各種不同數值方法的標準 ， 它可設計為具有任意指定的精度 。 N 次泰勒方法的最終全局誤差是 1()NOh? 階的，因此可選擇所需大小的 N ，使得誤差足夠小。 第六 章 龍格 庫塔（ Runge— Kutta 法） 龍格 庫塔（ Runge— Kutta）方法基本思想 回到 Euler 方法的基本思想 — 用差商代替導數 — 上來。改進的 Euler 公式可理解為 K 取 ( , )nnf x y ， 11( , )nnf x y??的平均值，其中 1 ( , )n n n ny y hf x y? ?? ，這種處理提高了精度。為此我們分析局部截斷誤差 11()nny x y??? ，因為 ()nny y x? ，所以（ ）可以化為 1 1 1 2 212121( ) ( )( , ( ) ) ( )( , ( ) ) ( 6 .3 )( , ( ) ) ( , ( ) )( , ( ) ) ( )nnn n nnnn n x n ny n ny y x h k kk f x y x y xk f x h y x h kf x y x h f x y xh k f x y x o h???????? ? ? ??????? ? ? ??? ???? ??? 其中 2k 在點 ? ?, ( )nnx y x 作了 Taylor展開?？梢宰C明，在 ? ?1,nnxx? 內只取 2 點的龍格 — 庫塔公式精度最高為 2階。一般地，只有 ky 用來計算 1ky? 。使用預估子和校正子的組合在每一步只與要進行兩次函數 (, )f t y 求值。此時值 1kp? 已知，基于點 11( , )kktf??， ( , )kktf 和新點1 1 1 1 1( , ) ( , ( , ) )k k k k kt f t f t p? ? ? ? ?? 構造 ( , ( ))f t yt 的新點的拉格朗日多項式，然后在區(qū)間11[ , ]kktt??上對該多項式積分，結果得 Simpon 公式： 1 1 1 1( 4 )3k k k k khy y f f f? ? ? ?? ? ? ? （ ） 誤差估計于校正 計算預估子和校正子的數值積分公式的誤差項都是 5()Oh 的，公式（ ）和（ ）的局部截斷誤差為 ( 5 ) 51 1 128( ) ( )90k k ky t p y c h? ? ??? （預估子的截斷誤差） （ ） ( 5 ) 51 1 11( ) ( )90k k ky t y y d h? ? ???? （校正子的截斷誤差） （ 7,5） 設 h 足夠小，使得 (5)()yt在區(qū)間 31[ , ]kktt??上近乎為常數，則可消去式（ ）和式（ ）中的 5 階導數項，結果為 1 1 1 128( ) ( )90k k k ky t p y p? ? ? ?? ? ? （ ） 公式 （ ）給出的預估子誤差聚集基于兩個計算值 1kp? 和 1ky? ，而沒有使用高階導數 (5)()yt，可用它來改進預估值。根據經驗，當小誤差遞減傳播時，結果穩(wěn)定；當小誤差遞增傳播時候結果不穩(wěn)定?？梢宰C明，其的步長 h 應該滿足以下條件 ( , )yh f t y? （ ） 第八章 一階微分方程組與高階微分方程的數值解法 一階微分方程組的數值解法 設有一階微分方程組的初值問題 12( , , , , ) ( 1 , 2 , )()i i mi ioy f x y y y imy a y? ?? ?? ?? （ ） 若記 12( , , , )Tmy y y y? ， 0 10 20 0( , , , )Tmy y y y? ， 12( , , , )Tmf f f f? ，則初值問題（ ）可寫成如下向量形式 常微分方程的高精度求解方法 安徽大學江淮學院 07 計算機 (1)班 18 0( , ) ( 8 . 2 )()y f x yy a y? ??? ?? 如果向量函數 ( , )f xy 在區(qū)域 D ： , ma x b y R? ? ? 連續(xù)，且關于 y 滿足 Lipschitz條件，即存在 L0，使得對 ? ?,x ab?? ，12, my y R? ，都有 1 2 1 2( , ) ( , )f x y f x y L y y? ? ? 那么問題（ ）在 ? ?,ab 上存在唯一解 ()y yx? 。 最后需要指出的是，在化學工程及自動控制等領域中，所涉及的常微分方程組初值問題常常是所謂的“剛性”問題。理論上的分析表明，求解剛性問題所選用的數值方法最好是對步長 h 不作任何限制。 如何解釋這一似乎與傳統的產品生命曲線理論相矛盾的現象昵 ?澳大利亞的斯蒂芬斯和莫賽觀察到購買耐用消費品的人大致可以分為兩類：一類是十分善于接受新事物的，稱為“創(chuàng)新型”顧客，他們往往從產品的廣告，制造商提供的