计算方法第6章 - 常微分方程数值方法

第六章 常微分方程数值方法

连续问题的离散处理——寻求微分方程的解y(t)在某些离散点上的值——

y(t)在ti处的近似值yi；

记号：y(ti)——所求函数在ii处的（准确）函数值， yi ——计算得到的ii处的函数（近似）值，

ti——节点， hi?ti?ti?1——步长；

b?a等步长：h?hi?；以下若不另作说明，一般总记h为等步长。

n § 6.1 初值问题的数值方法 考察微分方程初值问题：

a?t?b?y?(t)?f(t,y)?y(a)?y0?关的常数L，使

(6.1)

由微分方程理论可知：若函数f(t,y(t))关于y满足Lipschtz条件，即存在与y无

f(t,y1)?f(t,y2)?Ly1?y2初值问题的解存在且唯一；

6.1.1 Euler法及其变形 1、Euler法 由Taylor 展开式：

?t?[a,b],y1,y2?[c,d]

h2h2y??(?i)?y(ti)?hf(ti,y(ti))?y??(?i) y(ti?1)?y(ti)?hy?(ti)?22作局部化假设： yi?y(ti)，并略去O(h2)项，便有 y(ti?1)?yi?hf(ti,yi)

将此式右端作为y(ti?1)的近似yi?1，便得到公式

yi?1?yi?hf(ti,yi)?????? Euler公式

两者的差（即略去的O(h2)项）称为局部截断误差，记作E(ti,h)：

h2E(ti,h)?y??(?i)

2Euler公式也可由其他方法导出，例如由第四章数值导数公式，可有：

y(ti?1)?y(ti)y??(?i)?y(ti)??h,?i?[ti,ti?1]；

h2解出y(ti?1)，并由 y?(ti)?f(ti,y(ti)) 替换，便可得Euler公式。

又如，根据Newton-leibniz公式：

y(ti?1)?y(ti)??

ti?1tiy?(t)dt

1

将y?(t)以“0次插值多项式”替换，即以

y?(t)?y?(ti)?y??(?i)(t?ti)

代入积分，得到数值积分（左矩形）公式，及其误差：

ti?1h2y??(?i) ?y?(t)dt?hy?(ti)?ti2又得到与Taylor 展开式相同的表达式，从而又导出Euler公式。

r 几何意义——Eule折线法

若一个公式的局部截断误差为O(hp?1)，则称该公式的精度为p阶，或该公式为p阶公式。

Euler公式是1阶公式。

注意，以上的截断误差是在局部化假设的前提下得到的，即认定yi?y(ti)。倘若在每一步都按局部化假设，我们有Euler公式的

nh2b?a总体截断误差： E(y(b),h)?y(b)?yn??y??(?i)?hy??(?)

22i?12、后退Euler法

若取数值导数公式：

y(ti?1)?y(ti)h?y??(?i) y?(ti?1)?h2与前相同的推导过程，可以得到

h2y??(?i) y(ti?1)?y(ti)?hf(ti?1,y(ti?1))?2在局部化假设的前提下截去局部截断误差

h2y??(?i) E(ti,h)??2便得到后退Euler法公式：

yi?1?yi?hf(ti?1,yi?1)

注意到此公式中的右端也有yi?1，需要求解关于yi?1的方程才能得到yi?1。因此将这类公式称为隐式公式，而将可以通过直接计算得到的公式称为显式公式。后退Euler公式是一阶隐式公式，Euler公式是1阶显式公式

6.1.2 多步法 1、 梯形公式： 在式 y(ti?1)?y(ti)??ti?1tiy?(t)dt右端的积分中，取梯形积分公式，有

hh3y(ti?1)?y(ti)?[y?(ti?1)?y?(ti)]?y???(?i)

212由此，并据微分方程，可得：

2

h梯形公式 yi?1?yi?[f(ti,yi)?f(ti?1,yi?1)]

2h3局部截断误差： E(ti,h)??y???(?i)

12这是一个2阶隐式公式。 2、 Simpson公式：

在式 y(ti?1)?y(ti?1)??y?(t)dt右端的积分中，取Simpson积分公式，有

ti?1ti?1hh5(5)y(ti?1)?y(ti?1)?[y?(ti?1)?4y?(ti)?y?(ti?1)]?y(?i)

390由此，并据微分方程，可得：

hSimpson公式 yi?1?yi?1?[f(ti?1,yi?1)?4f(ti,yi)?f(ti?1,yi?1)]

3h5(5)y(?i)； 局部截断误差： E(ti,h)??90 与以前的公式不同，用Simpson公式计算yi?1，必须有前2步的函数值：yi?1和yi。因此这种方法称为2步方法，而为启动此算法所需的最初的2个函数值：

y0,y1称为表头。更一般的，若计算yi?1必须有至少前2步函数值，则这种方法称为多步法。具体地，若计算yi?1必须有前k步的函数值，则这种方法称为k步方法，而为启动该方法所需的最初的k个函数值：y0,y1,?,yk?1称为该方法的表头。与此相对，以前的方法计算yi?1，只须前1步的函数值yi，便称为单步方法。

因此，Simpson公式是2步、4阶、隐式方法。 3、 Adams方法（线性多步法） 在式 y(ti?1)?y(ti)??ti?1tiy?(t)dt右端的积分中，若取具有k+1节点的插值多

项式近似替代y?(t)作为被积函数，导出初值问题的求解方法称为Adams方法。

（1）显式Adams方法——Adams-Bashforth公式

取tj,j?i,i?1,?,i?k处的y?j?f(tj,yj)构造插值多项式取代y?(t)：

y?(t)?p(t)?R(t)

iip(t)?R(t)?其中?(t)?j?i?k?l(t)y?(tjj)??lj(t)f(tj,yj)j?i?k1y(k?2)(?i)?(t)(k?1)!

j?i?k?(t?tij)?(t?ti?k)??(t?ti)。由于?(t)?0,?t?[ti,ti?1]，有

3

y(ti?1)?y(ti)??iti?1tiy?(t)dt??ti?1ti[?lj(t)f(tj,yj)]dt??j?i?kiti?1ti1y(k?2)(?i)?(t)dt(k?1)!?j?i?k?[?ti?1tilj(t)dt]f(tj,yj)?ti?11y(k?2)(?i)??(t)dtti(k?1)!

由此（以后为方便计，记fi?f(ti,yi)），可得显式的多步法Adams-Bash forth公式及其局部截断误差

h(b0fi?b1fi?1???bkfi?k) AE(ti,h)?rhk?2y(k?2)(?i)yi?1?yi?这是k?1步、k?1阶的显式公式。

下表是k?0,1,?,4的 A,bj(j?0,1,?,k),r 的数值（注意：?bj?A）

k A b0 1 3 23 55 1901 b1 -1 -16 -59 -2774 b2 5 37 2616 b3 -9 -1274 b4 251 r 12 0 1 2 3 4

1 2 12 24 720 512 38 251720 95288 以下是 k?2 的公式推导过程：

p(t)?(t?ti?1)(t?ti?2)(t?ti)(t?ti?2)(t?ti)(t?ti?1)fi?fi?1?fi?2

(ti?ti?1)(ti?ti?2)(ti?1?ti)(ti?1?ti?2)(ti?2?ti)(ti?2?ti?1)作变量代换：t?ti?sh，以s为变量，当t的变化区间为[ti,ti?1]时，s的变化区间为[0,1]，且 dt?hds，有

(s?1)(s?2)s(s?1)fi?s(s?2)fi?1?fi?2]ds022

h?(23fi?16fi?1?5fi?2)121?ti?1tip(t)dt?h?[ti?111(4)h4(4)E(ti,h)?y(?i)?(t?ti)(t?ti?1)(t?ti?2)dt?y(?i)?s(s?1)(s?2)dsti03!3!

1?h4y(4)(?i)8 4

（2）隐式Adams方法——Adams-Moulton公式

若取tj,j?i?1,i,i?1,?,i?k?1处的y?j?f(tj,yj)构造插值多项式取代y?(t)，与前一样的方法，可得隐式的多步法Adams-Moulton公式及其局部截断误差

h??(b0f(ti?1,yi?1)?b1?fi???bkfi?k?1) AE(ti,h)?r?hk?2y(k?2)(?i)yi?1?yi?这是k步、k?1阶的隐式公式。

下表是k?0,1,?,4的 A,b?： ,?,k),r? 的数值（注意：?b?j(j?0,1j?A）

k A ? b0b1? 1 8 19 646 ? b2? b3? b4r? 0 1 2 3 4 1 2 12 24 720 1 1 5 9 251 -1 -5 -264 1 106 -19 ?12 ?112 ?124 ?19720 ?3160 ?述评：从表可见，对相同的k，A相同，而b??b,r?r，特别是：jjb?jA ?1，

(j?0,1,?,k)，而有若干j,bjA?1，因而在存在计算误差时，由前步导致的误差

显然隐式公式要比显式公式小（显式公式对前步的误差会被放大，而显隐式公式则不会），而且局部截断误差也是隐式公式要比显式公式小，结论：隐式公式的稳定性一般比显式公式好。

6.1.3 待定系数法

利用 Taylor展开式比较有关项的系数，可以直接导出公式——待定系数法。 例：求以下数值公式的系数?1,?0,?1,?2,使公式具有尽可能高的精度：

yi?1??yi?h(?0fi??1fi?1??2fi?2)

解：由于yi?1?y(ti?1)，因此由Taylor展开式，

111y(ti?1)?y(ti?h)?yi?hyi??h2yi???h3yi????h4yi(4)?O(h5)

2!3!4!同时，对所求公式右端各项也作相应的Taylor展开，并乘以相应的系数：

5

?yi??yih?0fi???yi ?h?0yi?

h?1fi?1?h?1y?(ti?h)?h?1yi??h2?1yi???h?2fi?2131h?1yi????h4?1yi(4)?O(h5) 23!8?h?2y?(ti?2h)?h?2yi??2h2?2yi???2h3?2yi????h4?2yi(4)?O(h5)

3!由于期望yi?1?y(ti?1)尽可能准确，比较各对应项的系数，可得方程组：

???????????1?0??1??2?1???1???23?0121 ，解之可得：????16 ； ??1?2?2?2?1121??2?5?1?2?2?1612?2注意到局部截断误差是E(ti,h)?y(ti?1)?yi?1，因此

13?1164(4)854(4)? E(ti,h)?h4yi(4)??hyi?hyi??O(h5)?h4yi(4)?O(h5)

4!3!128?3!12?显然，这就是k?2的显式Adams-Bashforth公式。

例：求以下数值公式的系数?1,?0,?1,?2,使公式具有尽可能高的精度：

yi?1??yi?3?h(?0fi??1fi?1??2fi?2)

与上例完全一样，比较系数，可得公式：

4 yi?1?yi?3?h(2fi?fi?1?2fi?2)

3145(5)hy?O(h6) 局部截断误差：E(ti,h)?45这个公式称为Milne公式，稍后我们将用到它。

综上所述，从公式的构造过程可见，微分方程初值问题的算法公式基本上源于：Taylor 展开式、数值积分或数值微分，而插值方法一般也是通过数值积分或数值微分完成的。

? Taylor 展开式——待定系数法直接来自于Taylor 展开式，通过比较系

数，形成线性方程组，取得公式与局部截断误差；

? 数值积分——前面的多步法等，基本上都源于此，例如梯形、Simpson方法等，而Adams方法用的是y(ti?1)?y(ti)??ti?1tiy?(t)dt，当y?(t)用不同

的点生成的不同的插值多项式代换时，便得到不同的公式，包括显式和隐式公式；

y(ti?1)?y(ti)h?y??(?i) 导出Euler公式，而由 ? 数值微分——由 y?(ti)?h2

6

y?(ti?1)?y(ti?1)?y(ti)h?y??(?i) 可导出后退Euler公式，若取3点数

h2值微分公式可导出更多公式：

1h2??3y(t0)?4y(t1)?y(t2)??y???(?)y?(t0)?2h3?y(ti?1)?4yi?3yi?1?2hf(ti?1,yi?1)?2hy???(?)33

1??y(t0)y?(t1)?2h?h2?y(t2)??y???(?)6 3hy(ti?1)?yi?1?2hf(ti,yi)?y???(?)31h2?y(t0)?4y(t1)?3y(t2)??y???(?)y?(t2)?2h3 312h?y(ti?1)??4yi?yi?1?2hf(ti?1,yi?1)??y???(?)39其中第二个公式就是“中点公式”；当然它们也可以由待定系数法取得：例如第3个公式：用待定系数法求以下隐式公式的系数及局部截断误差：

yi?1?Ayi?Byi?1?Chf(ti?1,yi?1)

由Taylor 展开式：

Ay(t)?i By(ti?1)?Ayi11By?Bhy??Bh2y???Bh3y????O(h4) ii2i6i1Chf(ti?1,y(ti?1))?Chy?(ti?1)?Chy??Ch2y???Ch3y????O(h4)ii2i11yi?hyi??h2yi???h3yi????O(h4) 比较：y(ti?1)?26为使公式具有尽可能高的精度，比较系数，得方程组：

?A?B?1???B?C?1?B?C?122??A?43???B??13 ?C?2?3??因此，公式是 y(ti?1)?yi?1?与之相应，有局部截断误差：

1?4yi?yi?1?2hf(ti?1,yi?1)?， 312?113???123????E(ti,h)?h3yi?????hy?hy??O(h4)??h3yi(4)?O(h4)

6239?63?

7

6.1.4 问题的性态与算法的稳定性

1、 问题的性态——问题对原始数据的敏感性 原始数据?————问题的应有的结果F(?)

~———问题的结果F(?~) 原始数据扰动?问题的条件数 ——相对误差之比的上确界：

~)][F(?)?F(?F(?)?m?Cond(F) ~][????~)][F(?)?F(?由微积分知识：

~][???F(?)?F?(?)??F?(?)?? F(?)F(?)初值问题：固定步长h，初值y0； 计算结果：y(t)?F(y0) 因此：

F(y0)?y(t)?y(t0?h)?y(t0)?hy?(t0)?O(h2)?y0?hf(t0,y0)?O(h2)?y0?hf(t0,y0)?(t0,y0) F?(y0)?1?hfy

Cond(F)?y(t0)1?hfy?(t0,y0)y0?hf(t0,y0)?y(t0)(1?hfy?(t0,y0))

y(t0?h)??0病态由于y(t0)?y(t0?h)，当h?0?fy?(t0,y0)? ；

??0良态或：y(t0)?y0扰动~y(t0)?y0??

~ y(t)?y?hf(t,y)?O(h2)y(t)?~y(t)?h~y?(t)?O(h2)

00000 ?~y(t0)?hf(t0,y0??)?O(h2) 又 f(t0,y0??)?f(t0,y0)??fy?(t0,y0)?O(?2)

得： ??~y(t)?y(t)??(1?hf?(t,y))

y00因此： fy?(t0,y0)?0??? 良态

??? 病态

fy?(t0,y0)?0例：y?(t)?2y(t) 方程的解 y(t)?ce2t 病态 （图示——） 例：y?(t)??2y(t) 方程的解 y(t)?ce?2t 良态 （图示——）

y(t)11?2,fy??2y(t)?， 例：y?(t)?[y(t)]2?ttt11 当y? 病态， 当y? 良态。（图示——）

2t2t

8

2、 方法的稳定性——方法：由初始误差导致的后期误差的可控性。

y(t?h)?yt(?h)h2?y???(?) 中点法：由数值导数公式 y??2h63h?y)?y(t?h)?2?hy()?t???得 y(t?h ( ) 可导出 3h3E(t,h)?y???(?)?O(h3) 中点法 yi?1?yi?1?2hf(ti,yi),3 此为2步2阶公式，但稳定性差。

例：y???10y,y(0)?1

真解：y(t)?e?10t,y(1)?0.453999?10?4

(2)，y1?y0?hf(t0,y0) （即一步Euler法）

(1)取2种不同的表头：y1?y(t1),下表中列出对不同的步长：h?1N(0) 计算y(1)的近似值yN，yN 按后退Euler法

(1)(2)计算获得，yN按中点法由初值及表头y1(1)获得，yN按中点法由初值及表头y1(2)

(1)(2)yN?yN获得，误差放大则是比值(1)。 (2)y1?y1h (0)后退Euler法yN (1) yN(2) yN误差放大 5741 11909 11121 10798 10?1 0.976562?10?3 ?0.266?103 ?0.238?104 10?2 10?3 10?4 0.725657?10?4 0.477118?10?4 0.456273?10?4 ?0.179?101 ?0.179?10?2 0.619?10?4 ?0.594?102 ?0.566?100 ?0.546?10?2 说明：

1、总体误差，由前可知——后退Euler法：

NN?h2?h???E?E(y(b),h)??E(ti.h)????y(?)??(b?a)y??(?) i??22i?1i?1??因此，对于后退Euler法，若h?10?2，（?(e?10t)???(0,100],?t?[0,1]）可计算得

1111E?10?2?100?；若h?10?4，可计算得 E?10?4?100?10?2 显然实际

2222情况良好；

而对于中点法，其总体误差

?h3?h2 E?E(y(b),h)??E(ti.h)???y???(?i)??(b?a)y???(?)

i?1i?1?3?3NN考虑表头y0,y1都用准确值（注意在浮点数系中运算，仍是有误差的，这可认为是初始误差），若h?10?2，（?(e?10t)????(0,1000],?t?[0,1]）可计算得E?10?510而h?10，可计算得E??13，

?43，而实际计算并非如此。

9

2、对于同是中点法：由此表可见，两个中点法的不同结果是由不同的表头导致的，实际的不同仅是y1(1),y1(2)的不同，而它们之间的误差（或不同）导致的y(1)计算结果的误差被放大5000倍甚至1万余倍。说明中点法并不是好的方法，尽管它是一个2阶方法，但实际计算不如后退Euler法这样的1阶方法。

定义：k步方法，若存在常数K,h0,使?h?h0

yi?yi?Kmaxyj?yj,0?j?k?1i?k,k?1,?

其中yi及yi是按此方法分别由表头?y0,y1,?,yk?1?及?y0,y1,?,yk?1?计算得到，则称此方法是稳定的。

此定义因未限制h0，故实际可要求h0?0，因此本定义又——“渐近稳定性”或“0稳定性”。

实际计算关心——对确定步长h，一个算法，每步计算误差是否被放大？ ——无法讨论一般方程——研究对象：典型方程：

y?(t)??y(t)a?t?b,y(a)?y0 其中：Re(?)?0。原因：

对n维常微分方程组：Y?(t)?F(t,Y(t))，DYF(t,Y(t))为n阶矩阵，它的特征值?反映了矩阵的性态，当Re(?)?0时，方程是良态，若Re(?)?0时则是病态（工程中，通常称之为“稳定性”， Re(?)?0时，系统稳定，否则，不稳定）。

定义：对于典型方程y?(t)??y(t)，及确定的h??h，一个k步方法，若由表头?y0,y1,?,yk?1?及?y0,y1,?,yk?1?计算得到yi及yi，存在常数C,0?C?1使

yi?yi?Cmaxyj?yj,0?j?k?1i?k,k?1,?

则称此方法对h是绝对稳定的，全体h的集合称为该方法的绝对稳定域。

事实上，也可将此定义改写成

yi?yi?Cmaxyi?j?yi?j,1?j?ki?k,k?1,?

对于单步法则为：

yi?yi?Cyi?1?yi?1,对典型方程，单步法 总有

yi?1?p(?h)yi?p(h)yi

例如：Euler法：yi?1?yi??hyi?(1??h)yi， p(h)?1?h

11yi， p(h)? 后退Euler法：yi?1?yi??hyi?1,yi?1?

1?h(1??h)(1??h)?h2y ，p(h)?2?h； 梯形法：yi?1?yi?(yi?yi?1),yi?1?i2?h2(1??h)2

10

i?1,2,?

由于单步法 yi?1?yi?1?p(h)(yi?yi)因

p(h)?1 则方法绝对稳定：

Euler法：1?h?1 ，在复平面上 是以 –1为中心，1为半径的圆内部；

1后退Euler法：?1， 即1?h?1，在复平面上 是以 1为中心，1为

1?h半径的圆的外部。

为稍后介绍方便，对一般的k步方法，写成

kk??j?0jyi?j?1?h??jfi?j?1?0j?0(?)

即： ?0yi?1??1yi????kyi?k?1??0fi?1??1fi????kfi?k?1?0

r yi?1?yi?hfi?0 例如 Eule法：

后退Euler法：yi?1?yi?hfi?1?0

中点法： yi?1?yi?1?2hfi?0

h梯形法： yi?1?yi?(fi?1?fi)?0

2hSimpson法： yi?1?yi?1?(fi?1?4fi?fi?1)?0 ；

3定义——算法(?)的特征多项式：

?(?,h)??(?)?h?(?)

其中 ?(?)???j?k?j??0?k??1?k?1????k?1???k

?(?)???j?k?j??0?k??1?k?1????k?1???k

j?0kj?0k分别称为算法（*）的第一、第二特征多项式。

定理：对h，若算法(?)的特征多项式的全体零点 ?i(i?1,2,?,k)都

满足?i?1，则 h 使该算法绝对稳定，全体这样的h的集合，称为该算法的

绝对稳定域。

例如：中点法的特征多项式 ?(?,h)??2??1h2? ，则?(?,h)?0 的

根?1,2?h?1?h2，由于

h?0:?1?h?1?h2?1;h?0:?1?h?1?h2??1，

可见，中点法总是不稳定的。

6.1.5 预估-校正方法

原因：显式方法——简单；但一般稳定性不好； 隐式方法——求解复杂；但一般稳定性好；

11

将两者结合——问题：如何组合

1、误差估计 —— 两个方法获得同一 y(ti)的不同近似yi，以此估计误差 （1） 同阶、同步长的两个不同方法 方法a：y(ti?1)?yi?1??hp?1y(p?1)(?1) 方法b：y(ti?1)?yi?1??hp?1y(p?1)(?2)

由于步长h一般较小，因此可认为 y(p?1)(?1)?y(p?1)(?2)，（实际上就是估计此因为一旦获得此估计值，便可得到E(t,h)?y(ti?1)?yi?1的估计值） y(p?1)(?1)的值，将两式相减，得

yi?1?yi?1?(???)hp?1y(p?1)(?) 即

hp?1y(p?1)(?)?因此

y(ti?1)?yi?1?yi?1?yi?1

(???)?(???)?yi?1?yi?1?

（2） 不同阶、同步长的两个不同方法，设 p?q

方法a：y(ti?1)?yi?1??hp?1y(p?1)(?1) 方法b：y(ti?1)?yi?1??hq?1y(q?1)(?2)

由于hq?1是hp?1的高阶小量，所以两者相加（减）时，后者可略去；将以上两式相减，可得

yi?1?yi?1??hp?1y(p?1)(?1) 因此：

y(ti?1)?yi?1?yi?1?yi?1

（3） 同阶、不同步长的同一方法，设步长 h?h

步长h：y(ti?1)?yi?1??hp?1y(p?1)(?1)

步长h：y(ti?1)?yi?1??hp?1y(p?1)(?2) 仿（1），将两式相减，得

yi?1?yi?1??[1?(h)p?1]hp?1y(p?1)(?1)

h所以

y(ti?1)?yi?1?(yi?1?yi?1) ；

[1?(h)p?1]h2、预估-校正法 （ Predictor-预估, Evaluation-计算, Correct-校正） （1） Heun 方法——改进Euler法

12

预估：Euler法— pi?1?yi?hf(ti,yi)

h校正：梯形法— yi?1?yi?[f(ti,yi)?f(ti?1,pi?1)]

2保持二阶精度： y(ti?1)?yi?1?O(h3)

证明：记y?为采用原始的梯形法获得的结果，即：y(ti?1)?y??O(h3)，由本节初所作的假设：函数f(t,y)关于变量y满足Lipschitz条件，

hh??y(ti?1)?yi?1?y(ti?1)??yi?[f(ti,yi)?f(ti?1,y?)]?[f(ti?1,pi?1)?f(ti?1,y?)]?

22??h??h ?y(ti?1)??yi?[f(ti,yi)?f(ti?1,y?)]??[f(ti?1,pi?1)?f(ti?1,y?)]

2??2h ?O(h3)?Lpi?1?y?

2h ?O(h3)?Lpi?1?y(ti?1)?y(ti?1)?y??O(h3)

2 (2) Adams-Bashforth-Moulton 方法

??取k?3的Adams-Bashforth显式公式作预估，Adams-Moulton隐式公式校正：

h pi?1?yi?[?9fi?3?37fi?2?59fi?1?55fi]

24h yi?1?yi?[fi?2?5fi?1?19fi?9f(ti?1,pi?1)]

24对应局部截断误差公式：

2515(5)195(5)hy(?1),y(ti?1)?yi?1??hy(?2)， y(ti?1)?pi?1?720720根据前已介绍的误差估计方法，有

19(yi?1?pi?1)。 y(ti?1)?yi?1??270由于在计算过程中，右端的值均已得到，故在每一步，均可对误差进行监控。在计算中，对于给定的相对误差界Re，

若

19yi?1?pi?1?Re270yi?1?Sm19yi?1?pi?1Re?270yi?1?Sm100hthen?h

2then2h?h

此处 Re（Relative error）——相对误差界, Ae（Absolute error）——绝对误差界， Sm (=Small)是一个与Ae相当的小值，目的是避免yi?1?0时出现的估计失真

~y?y（Ae实际上是函数y(t)的数量级，只须注意。 ?y）

Re?y?~y?yRe若相对误差过大，须步长减半（h2?h），新的进程中，将由yi,yi?1等原 有的值计算y(ti?h)的近似y1，而根据本方法，它是fi,f1,fi?1,f3的组合，

2i?i?i?222 13

对此中的原未取得的fi?12与fi?32的值，可采用插值法得到，但应注意，本方法是

4阶方法，因此，用插值法取得这些数值时应保证误差也应至少是O(h5)的，而5点插值多项式的余项： R(t)?f[t0,t1,t2,t3,t4,t](t?t0)(t?t1)(t?t2)(t?t3)(t?t4),且 t?tj?O(h)，因此用此5点插值方法可以实现误差是O(h5)的。为此，取最

近相邻的5点构造插值多项式，可得所求点处f(t,y)的近似值： f1(3fi?4?2f0?9fi?0?f6?1fi 5)3i?32i?0i?12821f1?(?5fi?4?28fi?3?70fi?2?140fi?1?35fi)； i?1282(3) Milne-Simpson方法

?以Milne公式为预估，Simpson公式做校正。

4pi?1?yi?3?h[2fi?2?fi?1?2fi]

3h(?tp yi?1?yi?1?[fi?1?4fi?f )]i1,?i13他们的局部截断误差公式分别为：

145(5)y(ti?1)?pi?1?hy(?1),

451y(ti?1)?yi?1??h5y(5)(?2)；

90与前相同，两者相减有：

295(5)hy(?1) yi?1?pi?1?90可导出预估公式-Milne公式的误差估计

28(yi?1?pi?1) y(ti?1)?pi?1?29与数值积分的Romberg公式的导出思想一样，现在设法将这部分的误差“归还”给预估的pi?1。由于在得到pi?1后，希望对pi?1进行修正，使修正后的值更接近

y(ti?1)再进入校正公式参与计算，这样得到的校正值有望更符合Simpson公式的实际效果。但由于尚未计算yi?1，只能用 yi?pi 近似替换yi?1?pi?1（实际上，是y(5)(?i)?y(5)(?i?1))，即：yi?pi?yi?1?pi?1，因此有

4?p?y?h[2fi?2?fi??i?312fi]?i?13?28? ?mi?1?pi?1?(yi?pi)

29?h?y?y?[fi?1?4fi?f(ti?1,mi?1)]i?1?i?13?

14

注：这是4步方法，必须有表头y0,y1,y2,y3才能开始计算。预估步与校正步当i?3便可执行，但在第一次的修正步计算m4中，由于此前无p3，因此只能令m4?p4，也即第一次因为无法修正所以不进行修正。

（4） 修正Milne-Hamming方法

解初值问题的Hamming公式：

13yi?1?(9yi?yi?2)?h[?fi?1?2fi?f(ti?1,yi?1)]

881E(ti,h)?y(ti?1)?yi?1??h5y(5)(?)

40以Milne公式为预估，Hamming公式进行校正。与前相同，两者的局部截断误差相减有：

yi?1?pi?1?1215(5)hy(?1) 360可导出预估公式-Milne公式的误差估计

112(yi?1?pi?1) y(ti?1)?pi?1?121校正公式-Hamming公式的误差估计：

9(yi?1?pi?1) y(ti?1)?yi?1??121与前一样，现在设法将误差“归还”给预估的pi?1，以及yi?1；由于通常将最终结果记作yi?1，因此将Hamming校正公式计算结果记作ci?1；与前Milne-Simpson方法一样，在修正pi?1时，ci?1尚未产生，因此只能使用ci，从而形成公式：

4?p?y?h[2fi?2?fi?1?2fi]i?3?i?13?112?mi?1?pi?1?(ci?pi)?121 ?

13?ci?1?(9yi?yi?2)?h[?fi?1?2fi?f(ti?1,mi?1)]88?9?y?c?(ci?1?pi?1)i?1i?1?121?注：仍与Milne-Simpson方法一样，这是一个4步方法，必须有表头

y0,y1,y2,y3才能开始计算。预估步与校正步当i?3便可执行，但在第一次的修正步计算m4中，由于此前无p3，因此只能令m4?p4。

注：如果从求解隐式方程的迭代法来看，这校正步可以认为是迭代了一步，而迭代的初值则是预估值，或是预估的一个修正值。因此从迭代收敛的角度看，应对初值提出一定的条件，通过研究，可以得到如下结果，各方法对步长的限制：

A-B-M4方法： h?

0.75

fy?(t,y)15

Milne-Simpson方法： h?修正Milne-Hamming方法：h?6.1.6 Runge-Kutta方法

0.45

fy?(t,y)0.69 。

fy?(t,y)将Heun 方法（改进Euler法）改写成如下形式：

1?y?y?(K1?K2)i?i?12? ?K1?hf(ti,yi)?K?hf(t?h,y?K)2ii1??此为一单步、2阶、显式公式（比较梯形：2阶隐式）。若注意此处的K1,K2，它们是两个不同的增量，Heun 方法则是从yi出发加上若干个增量的加权平均值后得到yi?1，将此推广到一般形式，yi出发加上m个增量的加权平价取得yi?1，便有：m-级Runge-Kutta公式

?yi?1?yi??1K1??2K2????mKm?K?hf(t,y)1ii???K2?hf(ti??2h,yi??21K1)?????????Km?hf(ti??mh,yi??m1K1??m2K2???mm?1Km?1)为导出该公式，先引入二元函数Taylor展开式：

f(t?h,y?k)?f(t,y)?(h??1???k)f(t,y)???(h?k)nf(t,y)?t?yn!?t?y

1???(h?k)n?1f(t??1h,y??2k),0??1,?2?1(n?1)!?t?y以二级Runge-Kutta公式为例：

?yi?1?yi??1K1??2K2? ?K1?hf(ti,yi)

?K?hf(t??h,y??K)ii1?2先考虑准确值y(ti?1)的展开：

1y(ti?1)?y(ti?h)?yi?hyi??h2yi???O(h3)

2其中

y??f

y???因此

df(t,y(t))?ft??fy?y??ft??ffy? dt 16

1 y(ti?1)?yi?hfi?h2(ft??ffy?)?O(h3)

2另一方面，考虑近似值yi?1的展开，由二元函数Taylor展开式，二级Runge-Kutta公式中，

K2?hfi??hft???K1fy??O(h2) 故应有：

??yi?1?yi??1hfi??2hfi??2?h2ft???2?h2ffy??O(h3)

将准确值y(ti?1)与近似值yi?1展开式作比较，应

??1??2?1?11 ???????22?22?若令 ?2?c，则

?1?1?c,??12c,??12c，

取c?1，则 ?1??2?1,????1， ——Heun 方法（改进Euler法）。

22取c?1，则 ?1?0,?2?1,????1 ——变形Euler法：

2??yi?1?yi?K2? K1?hf(ti,yi)??K1hK?hf(t?,y?)2ii?22?（注意此公式与中点法的区别）

以上两公式均为二阶二级Runge-Kutta公式。

常用的四阶四级Runge-Kutta公式（称为经典、古典或标准Runge-Kutta公式——Classical Runge-Kutta）

1?y?y?(K1?2K2?2K3?K4)i?i?16??K1?hf(ti,yi)?Kh? K2?hf(ti?,yi?1)?22??K2hK?hf(t?,y?)3ii?22?K4?hf(ti?h,yi?K3)??几何意义：若将上述RK?4?4公式改写为：

h yi?1?yi?(f1?2f2?2f3?f4)

6其中fj(j?1,2,3,4)为分别为原式中Kj(j?1,2,3,4)中相应的函数fj值。比较

17

在区间 [ti,ti?1]上的数值积分Simpson公式

ti?1hyi?1?yi??f(t,y)dt?yi?(fi?4fm?fi?1)

ti6其中 fi,fm,fi?1 分别表示f(t,y)在区间 [ti,ti?1]的左、中、右端点的值；它们分别有：fi?f1,fm?f2?f3,fi?1?f4；因此可简略地得到该RK?4?4公式的误差公式是 E(ti,h)??误差估计：

在实际计算中，由于Runge-Kutta方法是一个高精度的单步方法。因此，通常作为求解常微分方程的一个独立的方法使用。在计算时，常常要对计算值进行误差估计。以四阶四级的标准Runge-Kutta公式为例：

若用同一公式的不同步长方法进行误差估计，例如用步长h和h计算获得

2y(ti?1)的不同近似进行，为此必须计算11个f(t,y)值（其中：yi?yi?1 计算4 个f(t,y)值，yi?yi?121h5y(5)(?i) 2880。 ?yi?1计算8个f(t,y)值，而K1?hf(ti,yi)是重复的）

这比不估计误差增加7个f(t,y)值的计算量。

一个比较简单的方法是Runge-Kutta-Fehlberg方法

25140821971K1?K3?K4?K5),E(t,h)?O(h5)2162565410451666562856192yi?1?yi?(K1?K3?K4?K5?K6),E(t,h)?O(h6)13512825564305055K1?hf(ti,yi)yi?1?yi?(11h,yi?K1)44339K3?hf(ti?h,yi?K1?K2)8323212193272007296K4?hf(ti?h,yi?K1?K2?K3)132197219721974393680845K5?hf(ti?h,yi?K1?8K2?K3?K4)2165134104183544185911K6?hf(ti?h,yi?K1?2K2?K3?K4?K5)22725654104401128219712E(ti,h)?yi?1?yi?1?K1?K3?K4?K5?K63604275752405055K2?hf(ti?以此作误差估计，在计算过程中仅需计算6个f(t,y)值的计算量。

18

小结：

常微分方程求解的理论问题，主要是问题的性态和方法的稳定性。算法包括基本算法的构造（或：形成）及组合。前者就是通过Taylor展开式——Runge-Kutta法、待定系数法形成如Milne法等，或数值积分（例如：Adams方法，梯形法，Simpson法等），或数值微分（或称数值导数）形成；当然同一方法实际上可以通过多种方法导出。后者就是利用前面已得到的基本方法，进行适当组合，一方面利用显式方法简单，作为预估，利用隐式方法比较稳定，作为校正；进一步，同阶的不同方法，或同一方法不同步长的局部截断误差比较，将误差中的估计项——y(p?1)(?) 变换成可计算的?yi?1?pi?1?，得到实际误差估计值，或直接用来比较误差是否满足要求，或用来对对公式作修正。

6.1.7 常微分方程组、高阶常微分方程 1、常微分方程组初值问题：

常微分方程组初值问题的解法与常微分方程初值问题的解法几乎是一模一样的，唯一区别就在于所有的公式都用向量代替纯量。为此，记向量

?y1(t)??f1(t,y1(t),y2(t),?yn(t))?????y(t)f(t,y(t),y(t),?y(t))????12nY(t)??2?,F(t,Y(t))??2? ???????y(t)??f(t,y(t),y(t),?y(t))?12n?n??n?则可将常微分方程组初值问题

?(t)?f1(t,y1(t),y2(t),?,yn(t))?y1?y?(t)?f(t,y(t),y(t),?,y(t))2212n?????? ??y?(t)?f(t,y(t),y(t),?,y(t))n12n?n??y1(a)?y10,y2(a)?y20,?,yn(a)?yn0,改写成简单的向量形式：

Y(t)?F(t,Y(t)),Y(a)?Y0.

对此向量形式的初值问题的解法，只需将以前一维方程的数值方法改为向量形式便可。

例如，解常微分方程组的Euler方法，写成向量形式：

19

Yi?1?Yi?hF(ti,Yi)，

标准Runge-Kutta公式的向量形式：

1?y?y?(K1?2K2?2K3?K4)i?i?16?K1?hF(ti,Yi)?K1?h K?hF(t?,Y?)?2ii22?Kh?K2?hF(ti?,Yi?2)?22?K2?hF(ti?h,Yi?K3)?2、高阶方程：

对于高阶常微分方程初值问题，一般可以将它改写成常微分方程组，然后用前述的向量方法求解。

例如：3阶常微分方程初值问题

?,y??(a)?y0??， y???(t)?f(t,y(t),y?(t),y??(t)),y(a)?y0,y?(a)?y0令

?y?(t)?u(t)??u?(t)?v(t)?v?(t)?f(t,y(t),u(t),v(t))??y(t)???,Y(t)??u(t)?，

?v(t)????y0????? Y0??y0?y????0?并记

u(t)????v(t) F(t,Y(t))???,?f(t,y(t),u(t),v(t))???则原高阶方程就可以改写成如下形式，从而可按标准方法求解，

Y?(t)?F(t,Y(t)),Y(a)?Y0

6.2边值问题数值方法简介

6.2.1 差分方法 讨论常微分方程：

y??(t)?q(t)y(t)?f(t)a?t?b ，

其中 q(t)?0；对于一般方程：

y??(t)?p(t)y?(t)?q(t)y(t)?f(t)a?t?b

p(t)dt可以通过以下方法转化为上述形式。对方程乘以e?，得：

?p(t)dt???p(t)dt?f(t)e?p(t)dt， ?y?e???q(t)ye??令u?u(t)，使

20

?p(t)dtp(t)dtdudt?e??e?，即 dtdu从而一般方程等价于

p(t)dtp(t)dtd?dydt? ?f(t)e????q(t)ye?dt?dtdu?p(t)dtdt?e?在此方程两端同乘以 ，便得 du2?p(t)dt2?p(t)dtd2y?q(t)e?f(t)e du2利用u?u(t)的反函数（由ut??0知u?u(t)是关于t的单调函数，故有反函数），将 t 转化为 u 的函数，便可得

d2y?Q(u)y?R(u)， 2du此即为所讨论的形式。

第一边值条件 y(a)??,第二边值条件 y?(a)??,y(b)?? y?(b)??

?y?(a)??0y(a)??1第三边值条件 ??0?0,?0?0?0??0?0；

?y(b)??y(b)??01?1、 差分离散法——用数值导数替代导数

等分区间 步长：h?(b?a) ， 节点：ti?a?ihi?0,1,?,n，

n由 y??(ti)?q(ti)y(ti)?f(ti)

y(t?h)?2y(t)?y(t?h)h2(4)?y(?) 及数值微分：y??(t)?212hy(ti?1)?2y(ti)?y(ti?1)h2(4)?q(ti)y(ti)?f(ti)?y(?i) 得

12h2h2(4)y(?i) 注意：R?O(h2)；略去局部截断误差（Local truncation error）：Ri? 12得：在内点 ti(i?1,2,?,n?1) 处的离散逼近方程：

yi?1?2yi?yi?1?qiyi?fi,i?1,2,?,n?1， 2h或

yi?1?(2?qih2)yi?yi?1?fih2i?1,2,?,n?1。

求方程的数值解，y(ti)(i?0,1,?,n) 的近似yi尚差2个方程；由边界条件得到：

对第一边值条件，直接代入即可。

对第三边值条件（因第二边值为其特殊情况）利用数值导数公式：

1h2f?(x0)?[?3f(x0)?4f(x1)?f(x2)]?f???(?0)

2h3

21

1h2f?(x2)?[f(x0)?4f(x1)?3f(x2)]?f???(?2)

2h3略去截断误差部分O(h2)项,可得方程:

?(3?2?0h)y0?4y1?y2??2?1h

yn?2?4yn?1?(3?2?0h)yn??2?1h 至此，已有 n?1 个方程，可以求解。 2、 差分方程求解： 对第一边值条件，方程组：

1??(2?q1h)??y1??f1h2?????????21?(2?q2h)1????y2??f2h? ????????????????21?(2?qn?2h)1???yn?2??fn?2h??????21?(2?qn?1h)????yn?1??fn?1h???如果将上述系数矩阵按追赶法分解，简记为 T?LU

其中

?1???11????????l21???,U? L?? ?????n?21????????ln?11??n?1????显然：

?1??(2?q1h)??2,l2?1??1??1?12

?2??(2?q2h)?l2??1,?li??1,?i??1,i?3,4,?,n?1回顾线性方程组求解过程，可知，此追赶法与选列主元Gauss消去法是一致的，因此算法是稳定的。

但是，对第三边值条件，按上述介绍的方法离散化，将无法解决算法的稳定性问题。因此需要重新离散化。为此，将节点按下述方法确定： b?ahh?,ti?(a?)?ih,i?0,1,?,n,n?1，

n2在内点处（ti,i?1,2,?,n），仍按以前方法离散，而对边界点，则由

1h2y?(t)?[y(t?h)?y(t?h)]?y???(?)2h6 21hy(t)?[y(t?h)?y(t?h)]?y??(?)22将公式中的t分别以a和b代入，以h2代替h，分别代入第3边值条件，再将O(h2)

22

项作为截断误差略去，得

?y0?y1y?y1??00??1,h2?yn?yn?1y?yn?1??0n??1 h2将此按内点方程形式（非对角元为1 ）标准化，便形成方程组：

?1??0h?2???h?1?02?1???????????1h???1??y0???1??0h??2?y1???2?(2?q1h)1f1h???????????????????21?(2?qnh)1??y??fn?2h?h??n???1h1?0???2??yn?1??1??1??0h??h2??1?02??????? ?????? 与前面的分析类似，将系数矩阵按追赶法分解，简记为 T?LU 其中

?1???01????????l11???,U? L?? ?????n1????????ln?11??n?1?????h1?02??1，因此l?1??1,???(2?qh)?l??1, 注意到 ?0??1122?h?001?2由此推出 li??1,?i??1i?2,3,?,n，最后有ln?1??1,?n?1?0。注意，

?n?1?0的情况仅在?0??0?qi?0,i?1,2,?,n时出现。而此时事实上已出现

原方程退化为 y??(t)?f(t)a?t?b 形式，可直接由积分得到函数。除此以外，上述解线性方程组，总是稳定的。

6.2.2 打靶法（Shooting Method） 讨论方程

?y??(t)?f(t,y(t),y?(t))a?t?b （2.1） ??y(a)??,y(b)??将此问题的一个边值 y(b)?? 暂时搁置，而增加一个初值 y?(a)?x，形成初值问题：

?y??(t)?f(t,y(t),y?(t))a?t?b ? (2.2)

?y(a)??,y?(a)?x

23

该初值问题的解由于与新增加的参数x有关，因此将它记作y(t)?y(t,x)；显然，若该初值问题的解就是原边值问题的解，则应满足关系式，或方程：

y(b,x)?? (2.3) 因此，原边值问题的求解便改变成求解上述方程。 若用Newton法求解，则有迭代

y(b,xk)?? (2.4)

?y(b,xk)?x?y(b,xk)值y(b,xk)显然可由初值问题的求解获得。问题是如何求得。为此，将

?x初值问题全部加上参数x，即

xk?1?xk??y??(t,x)?f(t,y(t,x),y?(t,x))a?t?b ???y(a,x)??,y(a,x)?x注意，此处的y?,y??均为对变量t的导数。将此初值问题各项对参量x求导：方程

?y??(t,x)?f(t,y(t,x),y?(t,x))?t?f(t,y(t,x),y?(t,x))?y(t,x)???x?t?x?y?x

?f(t,y(t,x),y?(t,x))?y?(t,x)??y??x由于变量t与参量x无关，因此

?t?0，又对初值条件求导有： ?x?y(a,x)?y?(a,x)?0,?1，

?x?x若记

?y(t,x)?z(t,x)

?x注意到若函数y(t,x) 对变量t,x的各阶导数均连续，则导数与求导顺序无关。

综上所述，便形成关于变量t的初值问题（x仅作为参数）：

z??(t,x)??f(t,y(t,x),y?(t,x))?f(t,y(t,x),y?(t,x))z(t,x)?z?(t,x)a?t?b?y?y? (2.5)

z(a,x)?0,z?(a,x)?1解法：给定参变量xk，解初值问题 (2.2) ，得y(ti,xk),y?(ti,xk)，再解初值

?y(ti,xk)问题(2.5)，得z(ti,xk)?，从而得到解方程(2.3)的Newton迭代(2.4)所需的

?x?y(b,xk)??y(b,xk),。当迭代收敛，即收敛于x，使y(b,x)??，此时解初值问

?x题 (2.2) 得到的y(t)?y(t,x?)便是边值问题 (2.1) 的解。

24

本文来源：https://www.bwwdw.com/article/jgiw.html