1、第9章 常微分方程初值問題數(shù)值解法 9.1 引 言 本章要著重考察的一階方程的初值問題 （1.1）（1.2） 只要函數(shù) 適當(dāng)光滑譬如關(guān)于 滿足利普希茨(Lipschitz）條件 （1.3）理論上就可以保證初值問題（1.1），（1.2）的解 存在并且唯一. 1 所謂數(shù)值解法，就是尋求解 在一系列離散節(jié)點(diǎn) 上的近似值 . 相鄰兩個(gè)節(jié)點(diǎn)的間距 稱為步長. 如不特別說明，總是假定 為定數(shù)，這時(shí)節(jié)點(diǎn)為 . 初值問題（1.1），（1.2）的數(shù)值解法的基本特點(diǎn)是采取“步進(jìn)式”，即求解過程順著節(jié)點(diǎn)排列的次序一步一步地向前推進(jìn). 描述這類算法，只要給出用已知信息 計(jì)算 的遞推公式. 首先，要對(duì)方程（1.1）離散

2、化，建立求數(shù)值解的遞2推公式. 一類是計(jì)算 時(shí)只用到前一點(diǎn)的值 ，稱為單步法. 另一類是用到 前面 點(diǎn)的值 ，稱為 步法. 其次，要研究公式的局部截?cái)嗾`差和階，數(shù)值解 與精確解 的誤差估計(jì)及收斂性，還有遞推公式的計(jì)算穩(wěn)定性等問題. 3 9.2 簡(jiǎn)單的數(shù)值方法與基本概念 9.2.1 歐拉法與后退歐拉法 在 平面上，微分方程（1.1）的解 稱作它的積分曲線. 積分曲線上一點(diǎn) 的切線斜率等于函數(shù) 的值. 如果按函數(shù) 在 平面上建立一個(gè)方向場(chǎng)，那么，積分曲線上每一點(diǎn)的切線方向均與方向場(chǎng)在該點(diǎn)的方向相一致. 基于上述幾何解釋，從初始點(diǎn) 出發(fā)，先依方向場(chǎng)在該點(diǎn)的方向推進(jìn)到 上一點(diǎn) ，然后再從 依方向場(chǎng)的方

3、向推進(jìn)到 上一點(diǎn) ，循此前進(jìn)做出4一條折線 （圖9-1). 一般地，設(shè)已做出該折線的頂點(diǎn) ，過 依方向場(chǎng)的方向再推進(jìn)到 ，顯然兩個(gè)頂點(diǎn) 的坐標(biāo)有關(guān)系 圖9-15即 （2.1）這就是著名的歐拉（Euler）公式. 若初值 已知，則依公式（2.1）可逐步算出 例1 求解初值問題 （2.2）6 解 歐拉公式的具體形式為 取步長 ，計(jì)算結(jié)果見表9-1. 7 初值問題（2.2）有解 ，按這個(gè)解析式子算出的準(zhǔn)確值 同近似值 一起列在表9-1中，兩者相比較可以看出歐拉方法的精度很差. 還可以通過幾何直觀來考察歐拉方法的精度. 假設(shè) ，即頂點(diǎn) 落在積分曲線 上，那么，按歐拉方法做出的折線 便是 過點(diǎn) 的切線（

4、圖9-2）. 圖9-28 從圖形上看，這樣定出的頂點(diǎn) 顯著地偏離了原來的積分曲線，可見歐拉方法是相當(dāng)粗糙的. 為了分析計(jì)算公式的精度，通?？捎锰├照归_將 在 處展開，則有 在 的前提下， .于是可得歐拉法（2.1）的公式誤差 （2.3）稱為此方法的局部截?cái)嗾`差. 9如果對(duì)方程（1.1）從 到 積分，得 （2.4）右端積分用左矩形公式 近似，再以 代替 代替 也得到（2.1），局部截?cái)嗾`差也是（2.3）. 如果在（2.4）中右端積分用右矩形公式 近似，則得另一個(gè)公式 （2.5）稱為后退的歐拉法. 后退的歐拉公式與歐拉公式有著本質(zhì)的區(qū)別，后者是關(guān)于 的一個(gè)直接的計(jì)算公式，這類公式稱作是顯式的；10

5、然而公式（2.5）的右端含有未知的 ，它實(shí)際上是關(guān)于 的一個(gè)函數(shù)方程，這類公式稱作是隱式的. 隱式方程（2.5）通常用迭代法求解，而迭代過程的實(shí)質(zhì)是逐步顯示化. 設(shè)用歐拉公式 給出迭代初值 ，用它代入（2.5）式的右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得 然后再用 代入（2.5）式，又有 11如此反復(fù)進(jìn)行，得 （2.6）由于 對(duì) 滿足利普希茨條件（1.3）. 由（2.6）減（2.5）得 由此可知，只要 迭代法（2.6）就收斂到解 . 關(guān)于后退歐拉方法的公式誤差，從積分公式看到它與歐拉法是相似的. 12 9.2.2 梯形方法 為得到比歐拉法精度高的計(jì)算公式，在等式（2.4）右端積分中若用梯形求積公式近似

6、，并用 代替 代替 ，則得 （2.7）稱為梯形方法. 梯形方法是隱式單步法，可用迭代法求解. 同后退的歐拉方法一樣，仍用歐拉方法提供迭代初值，則梯形法的迭代公式為 13（2.8） 為了分析迭代過程的收斂性，將（2.7）式與（2.8）式相減，得 于是有 式中 為 關(guān)于 的利普希茨常數(shù). 14如果選取 充分小，使得 則當(dāng) 時(shí)有 ，這說明迭代過程（2.8）是收斂的. 15 9.2.3 單步法的局部截?cái)嗾`差與階 初值問題(1.1)，(1.2)的單步法可用一般形式表示為 （2.9）其中多元函數(shù) 與 有關(guān)，當(dāng) 含有 時(shí)，方法是隱式的，若不含 則為顯式方法，所以顯式單步法可表示為 （2.10） 稱為增量函數(shù)

7、，例如對(duì)歐拉法（2.1）有 它的局部截?cái)嗾`差已由（2.3）給出.16 對(duì)一般顯式單步法則可如下定義. 定義1 設(shè) 是初值問題(1.1)，(1.2)的準(zhǔn)確解，稱 （2.11）為顯式單步法（2.10）的局部截?cái)嗾`差. 之所以稱為局部的，是假設(shè)在 前各步?jīng)]有誤差. 當(dāng) 時(shí)，計(jì)算一步，則有 17 根據(jù)定義，顯然歐拉法的局部截?cái)嗾`差 即為（2.3）的結(jié)果. 這里 稱為局部截?cái)嗾`差主項(xiàng). 顯然 .的公式誤差. 所以，局部截?cái)嗾`差可理解為用方法（2.10）計(jì)算一步的誤差，也即公式（2.10）中用準(zhǔn)確解 代替數(shù)值解產(chǎn)生18 定義2 設(shè) 是初值問題(1.1)，(1.2)的準(zhǔn)確解，若存在最大整數(shù) 使顯式單步法(2

8、.10)的局部截?cái)嗾`差滿足 （2.12）則稱方法（2.10）具有 階精度. 若將（2.12）展開式寫成 則 稱為局部截?cái)嗾`差主項(xiàng). 以上定義對(duì)隱式單步法（2.9）也是適用的. 19 對(duì)后退歐拉法（2.5）其局部截?cái)嗾`差為 這里 ，是1階方法，局部截?cái)嗾`差主項(xiàng)為 . 20 對(duì)梯形法（2.7）有 所以梯形方法（2.7）是二階的，其局部誤差主項(xiàng)為 . 21 9.2.4 改進(jìn)的歐拉公式 梯形方法雖然提高了精度，但其算法復(fù)雜，在應(yīng)用迭代公式（2.9）進(jìn)行實(shí)際計(jì)算時(shí)，每迭代一次，都要重新計(jì)算函數(shù) 的值，而迭代又要反復(fù)進(jìn)行若干次，計(jì)算量很大，而且往往難以預(yù)測(cè). 為了控制計(jì)算量，通常只迭代一兩次就轉(zhuǎn)入下一步的

9、計(jì)算，這就簡(jiǎn)化了算法. 具體地，先用歐拉公式求得一個(gè)初步的近似值 , 稱之為預(yù)測(cè)值，預(yù)測(cè)值 的精度可能很差，再用梯形公式（2.7）將它校正一次，即按（2.8）式迭代一次得，這個(gè)結(jié)果稱校正值，而這樣建立的預(yù)測(cè)-校正系統(tǒng)通常稱為改進(jìn)的歐拉公式： 22預(yù)測(cè)校正（2.13）或表為下列平均化形式 23 例2 用改進(jìn)的歐拉方法求解初值問題（2.2）. 解 改進(jìn)的歐拉公式為 仍取 ，計(jì)算結(jié)果見表9-2. 同例1中歐拉法的計(jì)算結(jié)果比較，改進(jìn)歐拉法明顯改善了精度. 2425 9.3 龍格-庫塔方法 9.3.1 顯式龍格-庫塔法的一般形式 上節(jié)給出了顯式單步法的表達(dá)式其局部截?cái)嗾`差為對(duì)歐拉法 ，即方法為 階，若用

10、改進(jìn)歐拉法，它可表為 （3.1）26此時(shí)增量函數(shù) （3.2）它比歐拉法的 ，增加了計(jì)算一個(gè)右函數(shù) 的值，可望 . 若要使得到的公式階數(shù) 更大， 就必須包含更多的 值. 從方程（1.1）等價(jià)的積分形式（2.4），即 （3.3）若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求積公式27精度提高，必然要增加求積節(jié)點(diǎn)，為此可將（3.3）的右端用求積公式表示為點(diǎn)數(shù) 越多，精度越高，上式右端相當(dāng)于增量函數(shù) ，為得到便于計(jì)算的顯式方法，可類似于改進(jìn)歐拉法（3.1），（3.2），將公式表示為 （3.4）其中 （3.5）28這里 均為常數(shù). (3.4)和(3.5)稱為 級(jí)顯式龍格-庫塔（Runge-Kutta）法，

11、簡(jiǎn)稱R-K方法. 當(dāng) 時(shí)，就是歐拉法，此時(shí)方法的階為 . 當(dāng) 時(shí)，改進(jìn)歐拉法(3.1)，(3.2)也是其中的一種，下面將證明階 . 要使公式(3.4)，(3.5)具有更高的階 , 就要增加點(diǎn)數(shù) . 下面只就 推導(dǎo)R-K方法. 29 9.3.2 二階顯式R-K方法 對(duì) 的R-K方法，由（3.4），（3.5）可得到如下的計(jì)算公式 （3.6）這里 均為待定常數(shù)，希望適當(dāng)選取這些系數(shù)，使公式階數(shù) 盡量高. 根據(jù)局部截?cái)嗾`差定義，（3.6）的局部截?cái)嗾`差為 （3.7）30這里 . 為得到 的階 ，要將上式各項(xiàng)在 處做泰勒展開，由于 是二元函數(shù)，故要用到二元泰勒展開，各項(xiàng)展開式為 其中 （3.8）31將以

12、上結(jié)果代入（3.7）則有 要使公式（3.6）具有 階，必須使 32（3.9）即 （3.9）的解是不唯一的. 令 ，則得 這樣得到的公式稱為二階R-K方法，如取 ，則這就是改進(jìn)歐拉法（3.1）.33若取 ，則 . 得計(jì)算公式 （3.10）稱為中點(diǎn)公式，相當(dāng)于數(shù)值積分的中矩形公式. （3.10）也可表示為 34 的R-K公式(3.6)的局部誤差不可能提高到 . 把 多展開一項(xiàng)，從（3.8）的 看到展開式中 的項(xiàng)是不能通過選擇參數(shù)消掉的. 實(shí)際上要使 的項(xiàng)為零，需增加3個(gè)方程，要確定4個(gè)參數(shù) ，這是不可能的. 故 的顯式R-K方法的階只能是 ，而不能得到三階公式. 35 9.3.3 三階與四階顯式R

13、-K方法 要得到三階顯式R-K方法，必須 . 此時(shí)（3.4），（3.5）的公式表示為 （3.11）其中 及 均為待定參數(shù)，公式（3.11）的局部截?cái)嗾`差為 只要將 按二元函數(shù)泰勒展開，使 ，可36得待定參數(shù)滿足方程 （3.12）這是8個(gè)未知數(shù)6個(gè)方程的方程組，解也不是唯一的. 可以得到很多公式. 37滿足條件（3.12）的公式（3.11）統(tǒng)稱為三階R-K公式. 一個(gè)常見的公式為 此公式稱為庫塔三階方法. 繼續(xù)上述過程，經(jīng)過較復(fù)雜的數(shù)學(xué)演算，可以導(dǎo)出各種四階龍格-庫塔公式，下列經(jīng)典公式是其中常用的一個(gè)： 38 四階龍格-庫塔方法的每一步需要計(jì)算四次函數(shù)值 ，可以證明其截?cái)嗾`差為 . 39 例3

14、設(shè)取步長 ，從 直到 用四階龍格-庫塔方法求解初值問題 解 這里，經(jīng)典的四階龍格-庫塔公式（3.13）具有形式 4041 比較例3和例2的計(jì)算結(jié)果，顯然以龍格-庫塔方法的精度為高. 雖然四階龍格-庫塔方法的計(jì)算量（每一步要4次計(jì)算函數(shù) ）比改進(jìn)的歐拉方法（它是一種二階龍格-庫塔方法，每一步只要2次計(jì)算函數(shù) ）大一倍，但由于這里放大了步長 ，表9-3和表9-2 所耗費(fèi)的計(jì)算量幾乎相同. 龍格-庫塔方法的推導(dǎo)基于泰勒展開方法，因而它要求所求的解具有較好的光滑性質(zhì). 反之，如果解的光滑性差，那么，使用四階龍格-庫塔方法求得的數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法. 42 9.3.4 變步長的龍格-

15、庫塔方法 單從每一步看，步長越小，截?cái)嗾`差就越小，但隨著步長的縮小，在一定求解范圍內(nèi)所要完成的步數(shù)就增加了. 步數(shù)的增加不但引起計(jì)算量的增大，而且可能導(dǎo)致舍入誤差的嚴(yán)重積累. 因此同積分的數(shù)值計(jì)算一樣，微分方程的數(shù)值解法也有個(gè)選擇步長的問題. 在選擇步長時(shí)，需要考慮兩個(gè)問題： 1 怎樣衡量和檢驗(yàn)計(jì)算結(jié)果的精度？ 2 如何依據(jù)所獲得的精度處理步長？ 考察經(jīng)典的四階龍格-庫塔公式（3.13）. 43 從節(jié)點(diǎn) 出發(fā)，先以 為步長求出一個(gè)近似值，記為 ，由于公式的局部截?cái)嗾`差為 ，故有 （3.14）然后將步長折半，即取 為步長從 跨兩步到 ，再求得一個(gè)近似值 ，每跨一步的截?cái)嗾`差是 ，因此有 （3.15）比較（3.14）式和（3.15）式我們看到，步長折半后