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

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

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

3、計(jì)算簡(jiǎn)單，而隱式方法穩(wěn)定性好、精度高，故實(shí)際計(jì)算中，通常將兩者結(jié)合起來，利用同階顯式公式先做“預(yù)估”，再用隱式公式做一次“校正”，即所謂的預(yù)估-校正算法。本文綜述了預(yù)估-校正算法的理論基礎(chǔ)、算法構(gòu)造以及算法實(shí)現(xiàn)，并通過數(shù)值例子驗(yàn)證理論分析。</p><p>　　關(guān)鍵詞：有限差分法；預(yù)估-校正算法；算例</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 緒 論錯(cuò)誤!未定義書簽。</p><p>　　1.1 問題的背景錯(cuò)誤!未定義書簽。</p><p>　　1.2 有限差

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

9、法3</p><p>　　2.1.3 常微分方程初值問題的數(shù)值解法——多步法8</p><p>　　2.2 若干個(gè)預(yù)估-校正算法11</p><p>　　2.2.1 預(yù)估-校正算法的構(gòu)造11</p><p>　　2.2.2 預(yù)估-校算法的實(shí)現(xiàn)錯(cuò)誤!未定義書簽。</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、論錯(cuò)誤!未定義書簽。</p><p>　　致 謝錯(cuò)誤!未定義書簽。</p><p>　　參考文獻(xiàn)錯(cuò)誤!未定義書簽。</p><p>　　附 錄....................................................................... 28</p><p><b>　　1 緒

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

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

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

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

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

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

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

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

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

21、><p>　　這是因?yàn)?，任何高階方程或方程組的初值問題，經(jīng)過適當(dāng)?shù)淖儞Q，都可以化為1階方程組的初值問題，而1階方程組的初值問題寫成向量形式，并把向量形式用標(biāo)量形式來敘述，就是初值問題的基本形式（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）的一種最簡(jiǎn)單、最

23、基本的數(shù)值方法，它不一定作為一種獨(dú)立的求解方法在實(shí)際中使用，但卻提供了數(shù)值解法的本質(zhì)思想。</p><p>　　利用Taylor展開法，設(shè)（2.1.1-1）的解充分光滑，在已取定的處Taylor展開有</p><p>　　略去以上的項(xiàng)，并分別用，近似表示，，且根據(jù)方程可知，于是得近似計(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ù)推導(dǎo)求解初值問題（2.1.1-1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 若干個(gè)預(yù)估-校正算法</p><p>　　2.2.1 預(yù)估-校正算法的構(gòu)造[8]</p><p>　　由前面分析可知，具有相同步數(shù)

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

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

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

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

53、后計(jì)算，再使用校正算法（2.2.1-2）計(jì)算出的一次校正值。如果對(duì)重復(fù)上述校正過程，則可得到的次校正值。算法具體為：</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>　　上述計(jì)算程序通常記為，表示計(jì)算值。算法最后是以校正步結(jié)束的，此時(shí)我們已經(jīng)得到了，但是沒有進(jìn)一步計(jì)算的值。因而，在下一步的預(yù)估算式中， 仍用近似地代替。顯然，應(yīng)該比更為精確，因此，如果在上述計(jì)算程序結(jié)束時(shí)，我們計(jì)算出的值，則我們就可以用近似代替，我們記這種算法為，計(jì)算公式具體為：</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>　　注意，每一步的校正均需要計(jì)

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

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

58、;p>　　， （2.2.1-8）</p><p>　　， （2.2.1-9）</p><p>　　因而，預(yù)估方法（2.2.1-3）的局部截?cái)嗾`差是</p><p>　?。?.2.1-10）</p><p>　　考慮到校正算法一般是一個(gè)迭代的過程，校正算法的局部截?cái)嗾`差不能直接從

59、中得出。根據(jù)定義，局部截?cái)嗾`差表示計(jì)算時(shí)，所有已知值均取精確值時(shí)產(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è)，此時(shí)將（2.2.1-10）代入上式，則有</p><p>　　再設(shè)，并將上式代入（2.2.1-12）中，則可給出的估計(jì)。如此類推，我們有</p><p>　　應(yīng)用（2.2.1-9

61、），得局部截?cái)嗾`差</p><p>　　， （2.2.1-13）</p><p><b>　　由此，當(dāng)時(shí)，</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）看出，當(dāng)預(yù)估算法誤差階比校正算法誤差階低一階或者相等時(shí)，最后的誤差將與校正算法誤差同階。</p><p>　　有了局部截?cái)嗾`差的估計(jì)后，我們還可以進(jìn)一步給出誤差首項(xiàng)的估計(jì)，并對(duì)原算法作適當(dāng)?shù)男拚?/p>

63、使其最優(yōu)。令，，則從方程（2.2.1-10）和（2.2.1-14）知</p><p><b>　　，</b></p><p>　　由此可得誤差首項(xiàng)的估計(jì)</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）改變成在預(yù)估步之后就可使用的預(yù)估步誤差首項(xiàng)的估計(jì)</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算式。類似地，對(duì)校正值也可以作如下修正：</p><p>　?。?

66、.2.1-20）</p><p>　　如果將M算式（2.2.1-19）和（2.2.1-20）結(jié)合到預(yù)估-校正算法（2.2.1-4）或（2.2.1-5）中去，具體計(jì)算公式是（記為）：</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 預(yù)估-校正算法的實(shí)現(xiàn)

68、</p><p>　　常用的一種預(yù)估-校正方案是用4階4步Adams顯式公式做預(yù)測(cè)（Predictor），4階3步隱式公式做一次校正（Corrector），這就是下列的</p><p>　　Adams預(yù)估-校正格式（稱PC格式）</p><p><b>　　（2.2.2-1）</b></p><p>　　實(shí)際計(jì)算中，在預(yù)

69、估和校正之后分別緊跟一個(gè)為下一步做準(zhǔn)備的函數(shù)值計(jì)算（Evaluation）</p><p>　　于是，進(jìn)一步得Adams預(yù)估-校正格式（也稱PECE格式）為</p><p>　　如果再考慮修正技巧，由（2.2.2-1）第一式和第二式，以及注意到Adams的局部截?cái)嗾`差式(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>　　則這時(shí)比好，比好。不過前一式中的是未知的，計(jì)算實(shí)踐中可改成</p><p>　　于是，可得含修正技巧的Adams預(yù)估-校正格式（PMECME格式）：</p><p>　　該算法的穩(wěn)定性比較好，是最常用的預(yù)校算法。</p><p>　　另一個(gè)預(yù)估-校正方案是將4階4步Milne顯式公式做一次預(yù)測(cè)，4階3步Hamming隱式公式

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

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

74、-2）的第一特征多項(xiàng)式的零點(diǎn)的模均不超過1，并且模等于1的零點(diǎn)為單零點(diǎn)，則稱步公式（2.2.2-2）滿足根條件[9][10]。</p><p>　　可以驗(yàn)證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）的解對(duì)任意的，有</p><p><b>　　，</b></p><p><b>　　且初始條件滿足</b></p><p><b>　　， </b></p><p>　　則稱線性步法（2.2.2-2）是

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

78、>　　把線性步法（2.2.2-2）用于求解試驗(yàn)方程</p><p><b>　　， ，</b></p><p><b>　　這時(shí)線性步法簡(jiǎn)化為</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>　　對(duì)于給定的，如果特征方程（2.2.2-5）的根都滿足，則稱線性多步法(2.2.2-4)關(guān)于此是絕對(duì)穩(wěn)定的。</p><p>　　設(shè)為復(fù)平面內(nèi)的一個(gè)集合，如果對(duì)全部，線性多步法(2.2.2-4)都是絕對(duì)穩(wěn)定的，則稱此線性多步法有絕對(duì)穩(wěn)定域；與實(shí)軸之交稱為此線性多步法的絕對(duì)

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

81、式的絕對(duì)穩(wěn)定性好。</p><p>　　（2）Milne方法對(duì)任何它都不是決對(duì)穩(wěn)定的。但與Hamming方法組成預(yù)估-校正公式，則有較好的穩(wěn)定性。</p><p>　?。?）Simpson方法的絕對(duì)穩(wěn)定域?yàn)榭占Ｋ苑椒ú唤^對(duì)穩(wěn)定。實(shí)數(shù)對(duì)應(yīng)，不管取什么符號(hào)，誤差將連續(xù)增長(zhǎng)，所以在實(shí)際計(jì)算中不推薦Simpson方法。</p><p>　　上述各種預(yù)估-校正算法通常需要

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

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

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

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

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

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

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

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

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

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

92、，2008．</p><p>　　[8]余德浩，湯華中．微分方程數(shù)值解法[M]．北京：科學(xué)出版社，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]林成森．?dāng)?shù)值計(jì)算方法[M]．北京：科學(xué)出版社，1998．</p><p>　　[12]李立康，於崇華，朱政華．微分方程數(shù)值解法[M]．上海：復(fù)旦大學(xué)出版社，1999．</p><p>　　[13]關(guān)治，陸金甫．?dāng)?shù)值分析基礎(chǔ)[M]．北京：高等教育出版社，1998．</p><p>

95、　　[14]清華大學(xué)應(yīng)用數(shù)學(xué)系現(xiàn)代應(yīng)用數(shù)學(xué)手冊(cè)編寫組．現(xiàn)代應(yīng)用數(shù)學(xué)書冊(cè)（計(jì)算方法分冊(cè)）[M].北京：北京出版社，1990．</p><p>　　[15]王沫然，肖勁松．MATLAB5.X與科學(xué)計(jì)算[M]．北京：清華大學(xué)出版社，2000．</p><p><b>　　附錄</b></p><p>　　算例1在MATLAB運(yùn)行窗口中輸入</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運(yùn)行窗口中輸入</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法,改進(jìn)Euler法,4階Adams預(yù)估-校正算法,Hamming預(yù)估-校正算法,解微分方程</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; % 離散點(diǎn)數(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); % 采用改進(jì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預(yù)校算法</p><p>　　yt1 = f(x1); % 在精確解上取離散點(diǎn)</p>

104、;<p>　　subplot(2,2,1) % 將窗口分為一行四列，在第一個(gè)圖中作出各種方法的解的曲線</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); % 在精確解上取離散點(diǎn)</p><p>　　subplot(1,4,2) % 將窗口分為一行四列，在第二個(gè)圖中作出各種方法的解的曲線</p>

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

107、/p><p>　　xlim([0 2])</p><p>　　yt113=f(x113); % 在精確解上取離散點(diǎn)</p><p>　　subplot(1,4,3) % 將窗口分為一行四列，在第三個(gè)圖中作出各種方法的解的曲線</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); % 在精確解上取離散點(diǎn)</p><p>　　subplot(1,4,4) % 將窗口分為一行四列，在第四個(gè)圖中作出各種方法的解的曲線</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為取值范圍的一個(gè)端點(diǎn)</p><p><b>　　% N為區(qū)間的個(gè)數(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為取值范圍的一個(gè)端點(diǎn)</p><p><b>　　% N為區(qū)間的個(gè)數(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 離散點(diǎn)數(shù)</p><p>　　% 根據(jù)取決于KC = 1/0的誤差估計(jì)決定是否使用改正公式</p><p>　　% Hamming預(yù)估-校正公式為：</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>　　% 程序引導(dǎo)：</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); % 步長(zhǎng)</p><p>　　xt0 = x0+2*h;</p><p>　　[x,y] = Runge_Kutta(f,x0,xt0,y0,3); % 用 Runge-Kutta 得初始三點(diǎn)值</p><p>　　x = [x(1:3)' x(4