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