常微分方程初值問題的預估-校正解法【信息科學與技術(shù)專業(yè)】【畢業(yè)設(shè)計+文獻綜述+開題報告】

1、<p>　　本科畢業(yè)論文（設(shè)計）</p><p><b>　?。ā?0 　屆）</b></p><p>　　常微分方程初值問題的預估-校正解法</p><p>　　所在學院 </p><p>　　專業(yè)班級 信息與計算科學 <

2、;/p><p>　　學生姓名 學號 </p><p>　　指導教師 職稱 </p><p>　　完成日期 年 月 </p><p>　　摘要：求解常微分方程初值問題的有限差分法一般分為顯式方法和隱式方法。顯式方法

3、計算簡單，而隱式方法穩(wěn)定性好、精度高，故實際計算中，通常將兩者結(jié)合起來，利用同階顯式公式先做“預估”，再用隱式公式做一次“校正”，即所謂的預估-校正算法。本文綜述了預估-校正算法的理論基礎(chǔ)、算法構(gòu)造以及算法實現(xiàn)，并通過數(shù)值例子驗證理論分析。</p><p>　　關(guān)鍵詞：有限差分法；預估-校正算法；算例</p><p>　　Predictor - Corrector Algorithms f

4、or Initial-Value Problems for Ordinary Differential Equations </p><p>　　Abstract: Finite difference methods for solving Initial-value problems for ordinary differential equations can be divided into explicit

5、 methods and implicit methods generally. For the explicit methods, it can be calculated easily than the corresponding implicit methods, but the implicit methods has good stability and higher degree of accuracy, so in pra

6、ctical calculation, we usually combine them together, using explicit method as a "predictor", then implicit method as a "corrector", it is called</p><p>　　Key words: finite difference met

7、hods；predictor - corrector algorithms；numerical examples</p><p><b>　　目 錄</b></p><p>　　1 緒 論錯誤!未定義書簽。</p><p>　　1.1 問題的背景錯誤!未定義書簽。</p><p>　　1.2 有限差

8、分法的簡介錯誤!未定義書簽。</p><p>　　2 常微分方程初值問題的預估-校正算法錯誤!未定義書簽。</p><p>　　2.1 預估-校正算法的若干準備錯誤!未定義書簽。</p><p>　　2.1.1 常微分方程初值問題的基本概念錯誤!未定義書簽。</p><p>　　2.1.2 常微分方程初值問題的數(shù)值解法——單步

9、法3</p><p>　　2.1.3 常微分方程初值問題的數(shù)值解法——多步法8</p><p>　　2.2 若干個預估-校正算法11</p><p>　　2.2.1 預估-校正算法的構(gòu)造11</p><p>　　2.2.2 預估-校算法的實現(xiàn)錯誤!未定義書簽。</p><p>　　3 數(shù)值算例21

10、</p><p>　　3.3.1 算例 1............................................................ 21</p><p>　　3.3.2 算例 2............................................................ 22</p><p>　　4 結(jié)

11、論錯誤!未定義書簽。</p><p>　　致 謝錯誤!未定義書簽。</p><p>　　參考文獻錯誤!未定義書簽。</p><p>　　附 錄....................................................................... 28</p><p><b>　　1 緒

12、論</b></p><p>　　1.1 問題的背景</p><p>　　許多有關(guān)微分方程的教材都會提到發(fā)現(xiàn)海王星的故事。海王星的發(fā)現(xiàn)是人類智慧的結(jié)晶，也是常微分方程巨大作用的體現(xiàn)，體現(xiàn)了數(shù)學演繹法的強大威力。1781年發(fā)現(xiàn)天王星后，人們注意到它所在的位置總是和萬有引力定律計算出來的結(jié)果不符。于是有人懷疑萬有引力定律的正確性；但也有人認為，這可能是受另外一顆尚未發(fā)現(xiàn)的行星吸引所

13、致。當時雖有不少人相信后一種假設(shè)，但缺乏去尋找這顆未知行星的辦法和勇氣。23歲的英國劍橋大學的學生亞當斯挑下了這項任務(wù)，他利用引力定律和對天王星的觀測資料建立起微分方程，來求解和推算這顆未知行星的軌道。1843年10月21日，他把計算結(jié)果寄給格林威治天文臺臺長艾利，但艾利不相信“小人物”的成果，對其置之不理。兩年后，法國青年勒威耶也開始從事這項研究。1846年9月18日，他把計算結(jié)果告訴了柏林天文臺助理員卡勒。6日晚，卡勒果然在勒威耶預

14、言的位置上發(fā)現(xiàn)了海王星[1]。</p><p>　　對于數(shù)學，特別是數(shù)學的應(yīng)用，微分方程所具有的重大意義主要在于：很多物理與技術(shù)問題可以轉(zhuǎn)化為微分方程的求解問題，而這些問題僅有很少的一部分能通過初等積分法給出通解或通積分，大多數(shù)積分必須進行數(shù)值計算。所以，一開始就使用數(shù)值方法求解通常更加有效。常用的解常微分方程初值問題的數(shù)值方法有有限差分法和有限元法，其中有限差分法可以分為兩類：（1）單步法，例如Euler方法和

15、 Runge-Kutta方法；（2）多步法，例如線性多步法[2]。我們將同階的顯式公式與隱式公式相比較，發(fā)現(xiàn)前者使用方便，計算量較??；而后者一般需用迭代法求解，計算量大，但其局部截斷誤差較小，穩(wěn)定性較好。兩種方法各有其長處和不足。因此，常常將它們配合起來使用，以發(fā)揮它們的優(yōu)點，彌補各自的不足[3]。而將顯式公式和隱式公式聯(lián)合起來使用，前者提供預測值，后者將預測值加以校正，可以使數(shù)值解更精確。這種算法通常被稱作預估-校正算法（簡稱為PC算

16、法）。</p><p>　　1.2 有限差分法的簡介[4][5]</p><p>　　有限差分法是微分方程和積分方程數(shù)值解的方法。因其比較簡單，所以在實際應(yīng)用的數(shù)值方法中占很大比例。由于數(shù)字電子計算機只能存儲有限數(shù)據(jù)和作有限次運算，所以任何一種適用于計算機的解題方法，都必須把連續(xù)問題離散化，最終化成有限形式的線性代數(shù)方程組。</p><p>　　有限差分法將求解域

17、劃分為差分網(wǎng)格，用有限網(wǎng)格節(jié)點代替連續(xù)的求解域。對其用Taylor級數(shù)展開，把控制方程中的導數(shù)用網(wǎng)格節(jié)點上的函數(shù)值的差商來代替，從而建立以網(wǎng)格節(jié)點上的值為未知數(shù)的線性代數(shù)方程組。該方法是一種直接將微分問題變?yōu)榇鷶?shù)問題的近似數(shù)值解法，對于矩形區(qū)域，有限差分數(shù)值解是一種行之有效的方法，由于差分法數(shù)學概念直觀，表達簡單，因而成為發(fā)展較早且較成熟的數(shù)值方法。對于微分方程定解問題中的每一個差分方法，基本上按照（1）差分格式的建立；（2）差分格式的

18、求解；（3）算例；（4）差分格式的可解性、收斂性和穩(wěn)定性等4個方面展開。</p><p>　　本文將從以下兩個方面進行展開：首先，在第二部分，討論了常微分方程初值問題的理論，并引入常微分方程穩(wěn)定性的概念。然后給出求解的最簡單數(shù)值方法——Euler方法，實際計算中具有較好數(shù)值穩(wěn)定性的方法——向后Euler方法、梯形方法，以及更成熟、收斂更快速的方法——Runge-Kutta方法、Adams-Bashforth方法、

19、Adams-Moulton方法、Simpson方法、Milne方法和Hamming方法。實際計算中，通常將兩者結(jié)合起來，利用同階顯式公式先做“預估”，再用隱式公式做一次“校正”，甚至為了提高精度，還可利用事后誤差估計“修正”。通過這些操作，可以得到相當理想的預估-校正計算方案，如Adams預估-校正算法和Milne-Hamming預估-校正算法；最后，在第三部分，數(shù)值算例中，我們通過幾個算例驗證了我們的理論分析。</p>

20、<p>　　2 常微分方程初值問題的預估-校正算法</p><p>　　2.1 預估-校正算法的若干準備</p><p>　　2.1.1 常微分方程初值問題的基本概念</p><p>　　討論“初值問題”，通常集中于討論“1階初值問題”：</p><p><b>　?。?.1.1-1）</b></p

21、><p>　　這是因為，任何高階方程或方程組的初值問題，經(jīng)過適當?shù)淖儞Q，都可以化為1階方程組的初值問題，而1階方程組的初值問題寫成向量形式，并把向量形式用標量形式來敘述，就是初值問題的基本形式（2.1.1-1）。</p><p>　　設(shè)（2.1.1-1）中的是在上的連續(xù)函數(shù)；且關(guān)于滿足Lipschitz條件，即存在與無關(guān)的常數(shù)，使得</p><p><b>　

22、　, </b></p><p>　　成立，則，初值問題（2.1.1-1）的解存在、唯一、連續(xù)可微且連續(xù)依賴于初始條件[6]。</p><p>　　2.1.2 常微分方程初值問題的數(shù)值解法——單步法[7]</p><p>　　(1)Euler方法</p><p>　　Euler方法是求解初值問題（2.1.1-1）的一種最簡單、最

23、基本的數(shù)值方法，它不一定作為一種獨立的求解方法在實際中使用，但卻提供了數(shù)值解法的本質(zhì)思想。</p><p>　　利用Taylor展開法，設(shè)（2.1.1-1）的解充分光滑，在已取定的處Taylor展開有</p><p>　　略去以上的項，并分別用，近似表示，，且根據(jù)方程可知，于是得近似計算公式</p><p>　　， （2.1.2-1）<

24、/p><p>　　這就使得從初值和出發(fā)，用式（2.1.2-1）可得的近似值，再由，，用式（2.1.2-1）又可得的近似值，如此繼續(xù)，即可得初值問題（2.1.1-1）的Euler方法，式（2.1.2-1）稱為Euler公式，它是一種迭代公式，也是一種差分公式；它利用前一步的數(shù)值解算出下一步的數(shù)值解，所以Euler公式也稱為顯式單步法。</p><p>　　下面繼續(xù)推導求解初值問題（2.1.1-1

25、）的其他方法。現(xiàn)在考慮利用向前差商</p><p><b>　　和向后差商</b></p><p><b>　　，</b></p><p>　　整理并分別以，近似表示，，且令近似式為等式，即可得近似計算公式</p><p><b>　　和</b></p><

26、p><b>　?。?.1.2-2）</b></p><p>　　顯然，前者仍是Euler公式（2.1.2-1），而后者已有了本質(zhì)的不同，待求解的還隱含在計算公式等號右邊的中，故稱它為隱式Euler公式。對應(yīng)地，稱公式（2.1.2-1）為顯式Euler公式。隱式公式的計算，除非比較特殊，一般只好采用迭代法求解，像求解一元非線性方程的迭代法一樣，比較麻煩。所以我們換為另一種思路：</

27、p><p>　　把顯式（2.1.2-1）和隱式（2.1.2-2）作算術(shù)平均，則得</p><p><b>　?。?.1.2-3）</b></p><p>　　這稱為梯形公式。當然，它仍是隱式公式，但它的精度提高了。下面我們以此為例，來看看對隱式公式如何做迭代求解。對每一點以顯式公式（2.1.2-1）提供迭代初值，對梯形公式構(gòu)造迭代格式</p&

28、gt;<p>　　這時，如果迭代序列當時收斂于，便可取適當?shù)淖鳛榈慕浦?。這里容易證明，迭代對每一個，當時是收斂的。事實上，因為關(guān)于滿足Lipschitz條件，因此</p><p>　　可見，只要選取，使得，則有</p><p>　　它表示每多迭代一次，誤差總是小一點，故當時，便有，即迭代收斂。</p><p>　　既然迭代格式是收斂的，則取一次迭代（

29、即由顯式Euler公式算出并記為，代入梯形公式右端的，使梯形公式由隱式變?yōu)轱@式）可得如下的格式：</p><p><b>　　(2.1.2-4)</b></p><p>　　這稱為改進的Euler公式（其實質(zhì)是為了使用較高精度的梯形公式，借助顯式Euler公式使梯形公式顯式化）。</p><p>　　(2)Runge-Kutta方法</p

30、><p>　　設(shè)想要構(gòu)造的目標公式的形式為</p><p><b>　　（2.1.2-5）</b></p><p>　　其中，等均為待定參數(shù)，，式中用到個斜率值作加權(quán)平均以提高方法的階，故稱為級的Runge-Kutta方法。</p><p>　　對，式（2.1.2-5）實際就是Euler公式。</p><

31、p><b>　　對，有</b></p><p>　　， （2.1.2-6）</p><p>　　其中，為3個待定系數(shù)。通過確定這3個參數(shù)，使得式（2.1.2-6）具有2階精度。</p><p>　　為此，按局部截斷誤差的要求，分別計算出和，并考慮?，F(xiàn)設(shè)初值問題（2.1.1-1）的和解充分光滑，于是，一方面，由Ta

32、ylor展開有</p><p><b>　?。?.1.2-7）</b></p><p>　　其中以下標“”表示及其導數(shù)在處取值。另一方面，當解為式（2.1.2-6）時，由式（2.1.2-6）中的和作二元Taylor展開，并以下標“”表示及其導數(shù)在處取值，有</p><p><b>　　（2.1.2-8）</b></p

33、><p>　　由于考慮的是局部截斷誤差，應(yīng)有，故有，=。于是，由式（2.1.2-7）與式（2.1.2-8）兩邊相減可知，要使得式（2.1.2-6）具有二階精度，需且只需</p><p><b>　　（2.1.2-9）</b></p><p>　　至此，式（2.1.2-6）和式（2.1.2-9）聯(lián)立構(gòu)成2階精度的單步顯式公式族。由于它使用兩個斜率值，

34、通常稱為2級2階Runge-Kutta公式族。其中，一種常用的形式是在式（2.1.2-9）中取定，則，這時得所謂2級2階中點公式</p><p>　　對，類似上述做法，推導可得含8個未知數(shù)，6個方程的參數(shù)約束條件，從而得3級3Runge-Kutta方法族。其中一種確定的3級3階方法的具體形式為</p><p>　　對，推導可得含13個未知數(shù)，11個方程的參數(shù)約束條件，從而得4級4階Rung

35、e-Kutta方法族。其中，一種最著名、最常用的稱為經(jīng)典4階Runge-Kutta公式</p><p>　　解初值問題（2.1.1-1）的顯式單步法統(tǒng)一表達為</p><p>　　， （2.1.2-10）</p><p>　　其中稱為增量函數(shù)。則收斂問題定義如下：</p><p>　　由初值問題（2.1.1-

36、1）的顯式單步法（2.1.2-10）產(chǎn)生的近似解，如果對任一固定的均有，則稱單步法（2.1.2-10是收斂的。</p><p>　　定義使用于隱式單步法。注意定義中的是固定點，因為當?shù)猛瑫r有。顯然，若方法時收斂的，則在固定點的整體截斷誤差趨于零</p><p>　　設(shè)初值問題（2.1.1-1）的一種單步法（2.1.2-10）是階的，且函數(shù)關(guān)于滿足Lipschitz條件，即存在常數(shù)，使<

37、;/p><p>　　成立，則方法（2.1.2-10）是收斂的，且整體截斷誤差為。</p><p>　　容易驗證，Euler公式，改進Euler公式和Runge-Kutta公式均收斂。</p><p>　　對初值問題（2.1.1-1），如果用某一數(shù)值方法按某一固定步長進行計算時，在節(jié)點值上產(chǎn)生的擾動為，而由且僅由這個擾動引起的以后各節(jié)點值的變化不會逐步擴大，即，就稱這個數(shù)

38、值方法關(guān)于步長是絕對穩(wěn)定的。</p><p>　　對于具體的初值問題，由于各不相同，要逐一討論用于求解的數(shù)值公式的穩(wěn)定性比較困難，因此，在處置問題的理論研究中，采用一種比較特殊的方法，即通過一個簡單的所謂“試驗方程”用于檢驗各個數(shù)值方法的絕對穩(wěn)定性，從而建立一個絕對穩(wěn)定條件“框架”，以供求解具體問題的穩(wěn)定性分析作參照。試驗方程的形式如下：</p><p>　　（2.1.2-11）</

39、p><p>　　其中，是實數(shù)或復數(shù)，且規(guī)定為實數(shù)時，；為復數(shù)時，實部。</p><p>　　一般地，一個單步法（2.1.2-10）用于試驗方程（2.1.2-11），從計算一步得到，它們的關(guān)系表示成形式</p><p>　　其中，記號與所選的方法有關(guān)，則當時，方法（2.1.2-10）就是絕對穩(wěn)定的，從而使得的復變量所在的區(qū)域成為方法（2.1.2-10）的絕對穩(wěn)定域。絕對穩(wěn)

40、定域與實軸之交稱為絕對穩(wěn)定區(qū)間。</p><p>　　2.1.3 常微分方程初值問題的數(shù)值解法——多步法</p><p>　　求解初值問題（2.1.1-1）的另一類數(shù)值方法是多步法，即某一步解得計算不只依賴于前一步解得值，而是依賴于前兩步或兩步以上的解得值。</p><p>　　設(shè)線性公式關(guān)于和為線性，考慮步法，則一般形式為</p><p>

41、;　　(2.1.3-1) </p><p>　　并參照單步法局部截斷誤差的定義，推廣到線性多步法，給出局部截斷誤差的定義為</p><p><b>　　或由，寫成</b></p><p><b>　　又由，還可寫成</b></p><p><b>　　（2.1.3-2）<

42、;/b></p><p>　　特別地，當時，公式為顯式公式，當時，公式為隱式公式。</p><p>　　構(gòu)建具體的線性多步法，通常利用Taylor展開原理或插值型數(shù)值求積公式。下面，主要用前者，即去確定一般公式中的，使公式具有要求的精度，如要求具有階精度，即應(yīng)有，從而確定出應(yīng)有的。</p><p>　　實際應(yīng)有中幾個重要的線性多步法是：Adams（阿當姆斯）方

43、法、Milne（米納）方法、Hamming（哈明）方法、Simpson（辛普生）方法，它們均可通過指定的和等操作，解出。</p><p><b>　　取，可得</b></p><p><b>　　代入方程組，得</b></p><p><b>　　特別地，取，則</b></p><

44、p><b>　　，解之得</b></p><p><b>　　代入上式，可得</b></p><p>　　Adams 4步4階顯式公式（也稱Adams-Bashforth方法）</p><p><b>　　取，，得</b></p><p><b>　　(2.1.

45、3-3)</b></p><p><b>　　(2.1.3-4)</b></p><p>　　若取，則，再令，則類似上述做法，解出，于是可得</p><p>　　Adams 3步4階隱式公式（也稱Adams-Moulton方法）</p><p><b>　　取，，得</b></p&

46、gt;<p><b>　　(2.1.3-5)</b></p><p><b>　　(2.1.3-6)</b></p><p>　　類似上述做法，可分別得</p><p>　　Simpson 2步4階隱式公式</p><p><b>　　取，，得</b></p

47、><p><b>　　(2.1.3-7)</b></p><p><b>　　(2.1.3-8)</b></p><p>　　Milne 4步4階顯式公式</p><p><b>　　(2.1.3-9)</b></p><p>　　(2.1.3-10)<

48、;/p><p>　　Hamming 3步4階隱式公式</p><p>　　(2.1.3-11)</p><p>　　(2.1.3-12)</p><p>　　2.2 若干個預估-校正算法</p><p>　　2.2.1 預估-校正算法的構(gòu)造[8]</p><p>　　由前面分析可知，具有相同步數(shù)

49、的隱式方法的精度要比顯式方法高，例如步Adams內(nèi)插法的精度是階，而步Adams外插法的精度是階。此外，隱式方法的穩(wěn)定性也比顯式方法好，特別是當?shù)奶卣髦当容^大時，它的優(yōu)越性就更為突出。下面將討論隱式多步法的計算問題。</p><p>　　步隱式方法一般可以寫成下列形式：</p><p><b>　?。?.2.1-1）</b></p><p>　

50、　其中和已知。對一般的常微分方程（組）（2.1.1-1）而言，公式（2.2.1-1）是一個非線性隱式方程（組）。通?？梢杂玫椒ńo出它的近似解，迭代算式是:</p><p><b>　　(2.2.1-2)</b></p><p>　　其中是給定的初始近似值，。不難驗證，當初始近似值選得恰當，且時，迭代公式(2.2.1-2)是收斂的，即(2.2.1-2)的解收斂到(2

51、.2.1-1)的解，其中是的關(guān)于的Lipschitz常數(shù)。但是值得注意的是：初始近似選取的好壞將直接影響著迭代方法(2.2.1-2)的收斂快慢。解決這個問題的合理方法是用一個顯式線性多步方法計算，例如</p><p><b>　　（2.2.1-3）</b></p><p>　　這樣，將顯式公式（2.2.1-3）和隱式公式（2.2.1-2）聯(lián)合使用，前者提供預測值，而后

52、者將預測值加以校正，使數(shù)值解更精確。由此形成的算法通常被稱作預估-校正算法（簡稱為PC算法，稱公式（2.2.1-3）為P算式，公式（2.2.1-2）為C算式）。由于P算式已經(jīng)給出了的一個比較好的近似值，因此，一般只要對C算式迭代二至三次就可以使得足夠地小。顯然，校正次數(shù)過多的算法是不宜使用的，否則將給數(shù)值計算帶來很大的資源浪費。</p><p>　　PC算法具體實現(xiàn)如下：首先利用預估算法（2.2.1-3）給出，然

53、后計算，再使用校正算法（2.2.1-2）計算出的一次校正值。如果對重復上述校正過程，則可得到的次校正值。算法具體為：</p><p><b>　　P: ，</b></p><p><b>　　E：，</b></p><p><b>　　C：，</b></p><p><b

54、>　?。?.2.1-4）</b></p><p>　　上述計算程序通常記為，表示計算值。算法最后是以校正步結(jié)束的，此時我們已經(jīng)得到了，但是沒有進一步計算的值。因而，在下一步的預估算式中， 仍用近似地代替。顯然，應(yīng)該比更為精確，因此，如果在上述計算程序結(jié)束時，我們計算出的值，則我們就可以用近似代替，我們記這種算法為，計算公式具體為：</p><p><b>　　P

55、:，</b></p><p><b>　　E: ，</b></p><p><b>　　C: ， ，</b></p><p>　　E: （2.2.1-5）</p><p>　　注意，每一步的校正均需要計

56、算一次，這相當于一步預估算式的計算量。因此，合理的預估-校正算法應(yīng)該選取使在校正算法絕對穩(wěn)定區(qū)域內(nèi)，而在預估算法的絕對穩(wěn)定區(qū)域外，并且其值超過預估算法絕對穩(wěn)定的最大的倍數(shù)應(yīng)比校正步數(shù)大，否則就喪失了優(yōu)越性。</p><p>　　預估-校正算法的局部截斷誤差和線性穩(wěn)定性均將依賴于預估算法和校正算法，現(xiàn)在我們分析預估-校正算法（2.2.1-3）和（2.2.1-2）的截斷誤差。已知P算法（2.2.1-3）和C算法（2.

57、2.1-2）的特征多項式分別是</p><p>　　， （2.2.1-6）</p><p><b>　?。?.2.1-7）</b></p><p>　　如果分別用與表示步顯式算法（2.2.1-3）和步隱式算法（2.2.1-1），并假設(shè)它們的階分別為與，則當微分方程的解充分光滑時，我們有</p><

58、;p>　　， （2.2.1-8）</p><p>　　， （2.2.1-9）</p><p>　　因而，預估方法（2.2.1-3）的局部截斷誤差是</p><p>　　（2.2.1-10）</p><p>　　考慮到校正算法一般是一個迭代的過程，校正算法的局部截斷誤差不能直接從

59、中得出。根據(jù)定義，局部截斷誤差表示計算時，所有已知值均取精確值時產(chǎn)生的誤差。將（2.2.1-4）和（2.2.1-5）中的校正公式寫成如下統(tǒng)一形式：</p><p>　　， （2.2.1-11）</p><p><b>　　這里或，。由于</b></p><p><b>　　，</b></p><

60、;p><b>　　及已知的假設(shè)，則</b></p><p>　　， （2.2.1-12）</p><p>　　其中介于和之間，。先設(shè)，此時將（2.2.1-10）代入上式，則有</p><p>　　再設(shè)，并將上式代入（2.2.1-12）中，則可給出的估計。如此類推，我們有</p><p>　　應(yīng)用（2.2.1-9

61、），得局部截斷誤差</p><p>　　， （2.2.1-13）</p><p><b>　　由此，當時，</b></p><p>　　（2.2.1-14）</p><p><b>　　一般地</b></p><p>　　，

62、</p><p><b>　　其中</b></p><p>　?。?.2.1-15）</p><p>　　從（2.2.1-14）看出，當預估算法誤差階比校正算法誤差階低一階或者相等時，最后的誤差將與校正算法誤差同階。</p><p>　　有了局部截斷誤差的估計后，我們還可以進一步給出誤差首項的估計，并對原算法作適當?shù)男拚?/p>

63、使其最優(yōu)。令，，則從方程（2.2.1-10）和（2.2.1-14）知</p><p><b>　　，</b></p><p>　　由此可得誤差首項的估計</p><p><b>　　類似地，我們還有</b></p><p>　　， （2.2.1-16）</p><p>

64、;　　但是它只能在校正步結(jié)束后才能用。為此應(yīng)使用近似式</p><p><b>　　，</b></p><p>　　我們可將（2.2.1-16）改變成在預估步之后就可使用的預估步誤差首項的估計</p><p>　?。?.2.1-17）</p><p>　　將此代回（2.2.1-10），并令，則有</p>&

65、lt;p>　　， （2.2.1-18）</p><p><b>　　其中</b></p><p>　?。?.2.1-19）</p><p>　　這說明給出了比更好的近似。（2.2.1-19）通常被稱為修正算式或M算式。類似地，對校正值也可以作如下修正：</p><p>　?。?

66、.2.1-20）</p><p>　　如果將M算式（2.2.1-19）和（2.2.1-20）結(jié)合到預估-校正算法（2.2.1-4）或（2.2.1-5）中去，具體計算公式是（記為）：</p><p><b>　　P:，</b></p><p><b>　　M: ，</b></p><p><b&

67、gt;　　E: ，</b></p><p><b>　　C: ， ，</b></p><p><b>　　M: ，</b></p><p>　　E: ，如果。 （2.2.1-21）</p><p>　　2.2.2 預估-校正算法的實現(xiàn)

68、</p><p>　　常用的一種預估-校正方案是用4階4步Adams顯式公式做預測（Predictor），4階3步隱式公式做一次校正（Corrector），這就是下列的</p><p>　　Adams預估-校正格式（稱PC格式）</p><p><b>　?。?.2.2-1）</b></p><p>　　實際計算中，在預

69、估和校正之后分別緊跟一個為下一步做準備的函數(shù)值計算（Evaluation）</p><p>　　于是，進一步得Adams預估-校正格式（也稱PECE格式）為</p><p>　　如果再考慮修正技巧，由（2.2.2-1）第一式和第二式，以及注意到Adams的局部截斷誤差式(2.1.3-4)和(2.1.3-6)，于是有</p><p><b>　　，</

70、b></p><p><b>　　兩式相減可得</b></p><p><b>　　從而有</b></p><p><b>　　，</b></p><p>　　若用代替上式中的，并令</p><p><b>　　，</b>&l

71、t;/p><p>　　則這時比好，比好。不過前一式中的是未知的，計算實踐中可改成</p><p>　　于是，可得含修正技巧的Adams預估-校正格式（PMECME格式）：</p><p>　　該算法的穩(wěn)定性比較好，是最常用的預校算法。</p><p>　　另一個預估-校正方案是將4階4步Milne顯式公式做一次預測，4階3步Hamming隱式公式

72、做一次校正，并類似上述的設(shè)計，利用事后誤差估計做修正，可以得下列的所謂</p><p>　　Milne-Hamming預估-校正格式（PMECME）</p><p>　　考慮求解初值問題（2.1.1-1）的線性步方法，即把式（2.1.3-1）的線性多步法改寫成形式略有差異的步線性公式</p><p><b>　　（2.2.2-2）</b><

73、;/p><p>　　其中為常數(shù)，不全為零。如果引入多項式</p><p><b>　　和</b></p><p>　　則步線性公式（2.2.2-2）完全確定了和；反之，如果確定了和，也就確定了一個線性步法， 和分別稱為線性步法（2.2.2-2）的第一特征多項式和第二特征多項式。</p><p>　　如果線性步公式（2.2.2

74、-2）的第一特征多項式的零點的模均不超過1，并且模等于1的零點為單零點，則稱步公式（2.2.2-2）滿足根條件[9][10]。</p><p>　　可以驗證Adams顯式公式、Adams隱式公式、Simpson隱式公式和Hamming隱式公式均滿足根條件。</p><p>　　線性步法（2.2.2-2）是（）階相容的，則其收斂的充分必要條件是滿足根條件。</p><p&

75、gt;　　線性步法（2.2.2-2）穩(wěn)定的充分必要條件是它滿足根條件[11][12]。</p><p>　　現(xiàn)在假設(shè)由初值和線性步法求得的解表示如下：</p><p>　　， （2.2.2-3）</p><p><b>　　其中，。</b></p><p>　　設(shè)初值問題（2.1.1-1）中

76、連續(xù)，且關(guān)于滿足Lipschitz條件。如果線性多步法（2.2.2-3）的解對任意的，有</p><p><b>　　，</b></p><p><b>　　且初始條件滿足</b></p><p><b>　　， </b></p><p>　　則稱線性步法（2.2.2-2）是

77、收斂的[13]。</p><p>　　如果對應(yīng)于線性多步法（2.2.2-2）的第一特征多項式的零點都在單位圓內(nèi)或在單位圓周上，并在單位圓周上的根為單根，則稱線性多步法（2.2.2-2）滿足根條件。</p><p>　　設(shè)步線性公式（2.2.2-2）對應(yīng)地特征多項式，那么線性多步法（2.2.2-3）收斂的充分必要條件是線性多步法（2.2.2-2）滿足根條件。</p><p

78、>　　把線性步法（2.2.2-2）用于求解試驗方程</p><p><b>　　， ，</b></p><p><b>　　這時線性步法簡化為</b></p><p><b>　　(2.2.2-4)</b></p><p>　　記，相應(yīng)于(2.2.2-4)的特征方程為&

79、lt;/p><p><b>　?。?.2.2-5）</b></p><p>　　對于給定的，如果特征方程（2.2.2-5）的根都滿足，則稱線性多步法(2.2.2-4)關(guān)于此是絕對穩(wěn)定的。</p><p>　　設(shè)為復平面內(nèi)的一個集合，如果對全部，線性多步法(2.2.2-4)都是絕對穩(wěn)定的，則稱此線性多步法有絕對穩(wěn)定域；與實軸之交稱為此線性多步法的絕對

80、穩(wěn)定區(qū)間。</p><p>　　下面引述若干在應(yīng)用中可供參考的數(shù)據(jù)[14]</p><p>　　步數(shù)的Adams顯式與隱式方法的絕對穩(wěn)定區(qū)間如下表：</p><p>　　其中，Adams顯式方法在時就是單步Euler方法，所得絕對穩(wěn)定區(qū)間與單步法的結(jié)論一致；Adams隱式方法在時就是梯形公式，所得結(jié)論也一致；此外，同階的Adams方法中，隱式公式的絕對穩(wěn)定性比顯式公

81、式的絕對穩(wěn)定性好。</p><p>　?。?）Milne方法對任何它都不是決對穩(wěn)定的。但與Hamming方法組成預估-校正公式，則有較好的穩(wěn)定性。</p><p>　?。?）Simpson方法的絕對穩(wěn)定域為空集。所以方法不絕對穩(wěn)定。實數(shù)對應(yīng)，不管取什么符號，誤差將連續(xù)增長，所以在實際計算中不推薦Simpson方法。</p><p>　　上述各種預估-校正算法通常需要

82、通過編程上機計算。MATLAB提供了多個常微分方程求解的函數(shù)，采用自適應(yīng)變步長的求解方法， 即當解的變化較慢時采用較大的計算步長， 從而使得計算速度很快， 當方程的解變化得較快時， 積分步長會自動地變小， 從而使得計算的精度很高[15]。</p><p><b>　　3 數(shù)值算例</b></p><p><b>　　3.1 算例1</b>&l

83、t;/p><p><b>　　求解初值問題</b></p><p>　　當步長h=0.2時，用Euler法、改進Euler法、4階Adams預估-校正算法、Hamming預估-校正算法求該初值問題，數(shù)值解與精確解的絕對誤差如表3-1所示，數(shù)值解與精確解的圖的比較如圖3-1所示。從表3-1和圖3-1可以看出改進Euler法比Euler法的精度高；4階Adams預估-校正算法

84、和Hamming預估-校正算法比改進Euler法的精度高得多，而兩者之間精度差不多。</p><p>　　表3-1 算例1數(shù)值解與精確解之間的絕對誤差</p><p>　　圖3-1算例1的數(shù)值解與精確解的比較</p><p><b>　　3.2 算例2</b></p><p><b>　　求解初值問題<

85、/b></p><p>　　用4階Adams預估-校正算法和Hamming預估-校正算法求初值問題的結(jié)果，當步長h=1/5時，數(shù)值解與精確解以及它們的絕對誤差如表3-2所示，數(shù)值解與精確解的圖的比較如圖3-2所示；當h=1/16時，數(shù)值解與精確解的圖的比較如圖3-3所示。從表3-2、圖3-2和圖3-3可以看出，在4階Adams預估-校正算法中，當h=1/5或h=1/16時算法都收斂于精確值；而在Hammin

86、g預估-校正算法中，當h=1/16時算法收斂于精確值，當h=1/5時算法不收斂于精確值。這個現(xiàn)象說明Hamming預估-校正算法對步長的依賴性比4階Adams預估-校正算法強。</p><p>　　表3-2 算例2數(shù)值解與精確解以及它們的絕對誤差</p><p>　　圖3-2 當h=1/5算例2的數(shù)值解與精確解的比較</p><p>　　圖3-3 當h=1/16

87、時算例2的數(shù)值解與精確解的比較</p><p><b>　　4 結(jié)論</b></p><p>　　常微分方程是建立物理、生物科學以及工程學的模型時所使用的最重要的數(shù)學工具之一。本文考慮的是求解常微分方程的數(shù)值方法。首先討論了常微分方程初值問題的理論，并到引入微分方程穩(wěn)定性的概念。給出求解的最簡單數(shù)值方法——Euler方法，雖然它不是有效的數(shù)值方法，但研究它有利于引入

88、和理解許多重要思想。隨后談?wù)撛趯嶋H計算中具有較好數(shù)值穩(wěn)定性的方法——向后Euler方法、梯形方法。最后著重于更成熟、收斂更快速的方法——Runge-Kutta方法、Adams-Bashforth方法、Adams-Moulton方法、Simpson方法、Milne方法和Hamming方法。由于顯式公式與隱式公式和各有其長處和不足，故實際計算中，通常將兩者結(jié)合起來，利用同階顯式公式先做“預估”，再用隱式公式做一次“校正”，甚至為了提高精度，

89、還可利用事后誤差估計“修正”。通過這些操作，可以得到相當理想的預估-校正計算方案，如Adams預估-校正算法和Milne-Hamming預估-校正算法。最后通過MATLAB編程計算，給出數(shù)值例子，驗證了理論分析。</p><p><b>　　參考文獻</b></p><p>　　[1]Kendall Atkinson，韓渭敏著，王國榮，徐兆亮，孫劼譯．數(shù)值分析導論(第

90、三版)[M]．北京：人民郵電出版社，2009．</p><p>　　[2]孫志忠，袁慰平，聞?wù)鸪酰當?shù)值分析(第二版)[M]．南京：東南大學出版社，2002．</p><p>　　[3]李榮華，馮果忱．微分方程數(shù)值解法[M]．北京：高等教育出版社，1995．</p><p>　　[4] Arieh Iserles．微分方程數(shù)值分析基礎(chǔ)教程[M]．北京：清華大學出版社，

91、2005．</p><p>　　[5] 李瑞遐，何志慶．微分方程數(shù)值方法[M]．上海：華東理工大學出版社，1999．</p><p>　　[6]Arieh Lserles著，劉曉艷，劉學深譯．微分方程數(shù)值分析基礎(chǔ)教程[M]．北京：清華大學出版社，2005．</p><p>　　[7]鄭咸義，姚迎新，雷秀仁，陸子強，黃鳳輝．應(yīng)用數(shù)值分析[M]．廣州：華南理工大學出版社

92、，2008．</p><p>　　[8]余德浩，湯華中．微分方程數(shù)值解法[M]．北京：科學出版社，2003．</p><p>　　[9]Richard L.Burden and J.Douglas．NUMERICAL ANALYSIS[M]．Wadsworth Group．Brooks，2001． </p><p>

93、　　[10]Ascher U.M. and Petzold L.R．Computer methods for ordinary differential equations and differential-algebraic equations[M]．Society for Industrial 

94、Mathematics，1998．</p><p>　　[11]林成森．數(shù)值計算方法[M]．北京：科學出版社，1998．</p><p>　　[12]李立康，於崇華，朱政華．微分方程數(shù)值解法[M]．上海：復旦大學出版社，1999．</p><p>　　[13]關(guān)治，陸金甫．數(shù)值分析基礎(chǔ)[M]．北京：高等教育出版社，1998．</p><p>

95、　　[14]清華大學應(yīng)用數(shù)學系現(xiàn)代應(yīng)用數(shù)學手冊編寫組．現(xiàn)代應(yīng)用數(shù)學書冊（計算方法分冊）[M].北京：北京出版社，1990．</p><p>　　[15]王沫然，肖勁松．MATLAB5.X與科學計算[M]．北京：清華大學出版社，2000．</p><p><b>　　附錄</b></p><p>　　算例1在MATLAB運行窗口中輸入</p

96、><p>　　y=dsolve('(Dy)-y+x^2-1=0','y(0)=0.5','x')</p><p><b>　　輸出結(jié)果為</b></p><p><b>　　y =</b></p><p>　　1+x^2+2*x-1/2*exp(x)<

97、;/p><p>　　算例2在MATLAB運行窗口中輸入</p><p>　　y=dsolve('(Dy)+8*y-4*x^2+7*x+1=0','y(0)=1','x')</p><p><b>　　輸出結(jié)果為</b></p><p><b>　　y =</b&

98、gt;</p><p>　　1/2*x^2-x+exp(-8*x)</p><p><b>　　主程序：</b></p><p>　　% 分別用Euler法,改進Euler法,4階Adams預估-校正算法,Hamming預估-校正算法,解微分方程</p><p>　　clear, clf</p><p

99、>　　format long</p><p>　　x0 = 0; % 積分區(qū)域開始值</p><p>　　xt = 2; % 積分區(qū)域終止值</p><p>　　y0 = 0.5; % y初值</p><p>　　N = 11; % 離散點數(shù)</p><p>　　fun = inline('y-x^2+1

100、','x','y'); % 需要求解的微分方程</p><p>　　f = inline('1+x.^2+2*x-1/2*exp(x)','x'); % 方程的精確解析解</p><p><b>　　% 算例2</b></p><p>　　% x0 = 0; </p&

101、gt;<p>　　% xt = 3; </p><p>　　% y0 = 1; </p><p>　　% N = 16; % N = 16; </p><p>　　% fun = inline('-8*y+4*x^2-7*x-1','x','y'); </p><p>　　% f

102、= inline('1/2*x.^2-x+exp(-8*x)','x'); </p><p>　　[x1,y1]=Euler_1(fun,x0,xt,y0,N); % 采用Euler法</p><p>　　[x2,y2]=Euler_2(fun,x0,xt,y0,N); % 采用改進Euler法</p><p>　　[x113,y11

103、3]=ode113(fun,[x0:(xt-x0)/(N-1):xt],y0);</p><p>　　% 采用Adams-Bashforth-Moulton法</p><p>　　[xH,yH] = Hamming(fun,x0,xt,y0,N); % 采用Hamming預校算法</p><p>　　yt1 = f(x1); % 在精確解上取離散點</p>

104、;<p>　　subplot(2,2,1) % 將窗口分為一行四列，在第一個圖中作出各種方法的解的曲線</p><p>　　plot(x1,yt1,'*r',x1,y1,'og')</p><p>　　legend('精確解','Euler解')</p><p>　　title('

105、Euler法')</p><p><b>　　hold on</b></p><p>　　xlim([0 2])</p><p>　　yt2 = f(x2); % 在精確解上取離散點</p><p>　　subplot(1,4,2) % 將窗口分為一行四列，在第二個圖中作出各種方法的解的曲線</p>

106、<p>　　plot(x2,yt2,'*r',x2,y2,'og')</p><p>　　legend('精確解','改進Euler解')</p><p>　　title('改進Euler法')</p><p><b>　　hold on</b><

107、/p><p>　　xlim([0 2])</p><p>　　yt113=f(x113); % 在精確解上取離散點</p><p>　　subplot(1,4,3) % 將窗口分為一行四列，在第三個圖中作出各種方法的解的曲線</p><p>　　plot(x113,yt113,'*r', x113, y113,'og

108、9;)</p><p>　　legend('精確解','ABM解')</p><p>　　title('4階Adams法')</p><p><b>　　hold on</b></p><p>　　xlim([0 2])</p><p>　　ytH

109、= f(xH); % 在精確解上取離散點</p><p>　　subplot(1,4,4) % 將窗口分為一行四列，在第四個圖中作出各種方法的解的曲線</p><p>　　plot(xH,ytH,'*r',xH,yH,'og')</p><p>　　legend('精確解','Hamming解')<

110、;/p><p>　　title('Hamming法')</p><p><b>　　hold off</b></p><p>　　xlim([0 2])</p><p>　　x1,y1,yt1,err1=abs(yt1-y1),</p><p>　　x2,y2,yt2,err2=abs

111、(yt2-y2),</p><p>　　x113,y113,yt113,err113=abs(yt113-y113),</p><p>　　xH,yH,ytH,errH=abs(ytH-yH)</p><p><b>　　子程序：</b></p><p>　　% Euler_1()</p><p>

112、;　　function [x,y]=Euler_1(fun,x0,xt,y0,N)</p><p>　　% fun為一階微分方程的函數(shù)</p><p>　　% x0,y0為初始條件</p><p>　　% b為取值范圍的一個端點</p><p><b>　　% N為區(qū)間的個數(shù)</b></p><p&g

113、t;　　x=zeros(1,N);y=zeros(1,N);</p><p>　　x(1)=x0;y(1)=y0;</p><p>　　h=(xt-x0)/(N-1);</p><p>　　for n=1:N-1</p><p>　　x(n+1)=x(n)+h;</p><p>　　y(n+1)=y(n)+h*feva

114、l(fun,x(n),y(n));</p><p><b>　　end</b></p><p>　　% Euler_2()</p><p>　　function [x,y]=Euler_2(fun,x0,xt,y0,N)</p><p>　　% fun為一階微分方程的函數(shù)</p><p>　　%

115、x0,y0為初始條件</p><p>　　% b為取值范圍的一個端點</p><p><b>　　% N為區(qū)間的個數(shù)</b></p><p>　　x=zeros(1,N);y=zeros(1,N);</p><p>　　x(1)=x0;y(1)=y0;</p><p>　　h=(xt-x0)/(N

116、-1);</p><p>　　for n=1:N-1</p><p>　　x(n+1)=x(n)+h;</p><p>　　z0=y(n)+h*feval(fun,x(n),y(n));</p><p>　　y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));</p&g

117、t;<p><b>　　end</b></p><p>　　% Hamming()</p><p>　　function [x,y]=Hamming(f,x0,xt,y0,N,KC)</p><p>　　% Hamming——解常微分方程</p><p>　　% 用hamming方法解向量微分方程 y’(t

118、) = f(t,y(t))</p><p>　　% x向量范圍[x0,xt], 初值 y0,N 離散點數(shù)</p><p>　　% 根據(jù)取決于KC = 1/0的誤差估計決定是否使用改正公式</p><p>　　% Hamming預估-校正公式為：</p><p>　　% p(n+1)=y(n-3)+4*h/3*(2*f(n-2)-f(n-1)+

119、2*f(n))</p><p>　　% m(n+1)=p(n+1)+112/121*(c(n)-p(n))</p><p>　　% c(n+1)=1/8*[9*y(n)-y(n-2)+3*h*(-f(n-1)+2*f(n)+f(x(n+1),m(n+1)))]</p><p>　　% y(n+1)=c(n+1)-9/121*(c(n+1)-p(n+1))

120、</p><p><b>　　% 程序引導：</b></p><p>　　% p1=y(1)+4*h/3*(2*f(2)-f(3)+2*f(4))</p><p>　　% m1=p1+112/121*(c-p)</p><p>　　% c1=1/8*[9*y(4)-y(2)+3*h*(-f(3)+2*f(4)

121、+f(x(5),m1))]</p><p>　　% y(5)=c1-9/121*(c1-p1)</p><p>　　if nargin < 6,</p><p>　　KC = 1; % 缺省使用更正公式</p><p><b>　　end </b></p><p>　　if nargin

122、 < 5 | N <= 0</p><p>　　N = 100; % 最大迭代數(shù)缺省為100</p><p><b>　　end </b></p><p>　　if nargin < 4</p><p>　　y0 = 0; % 初值缺省為0</p><p><b>　　

123、end </b></p><p>　　h = (xt-x0)/(N-1); % 步長</p><p>　　xt0 = x0+2*h;</p><p>　　[x,y] = Runge_Kutta(f,x0,xt0,y0,3); % 用 Runge-Kutta 得初始三點值</p><p>　　x = [x(1:3)' x(4