常微分方程的数值解法

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

§6.0 引言

§6.1 算法构造的主要途径 §6.2 Runge-Kutta Method算法 §6.3 线性多步法

§6.4 线性多步法的一般形式 §6.5 算法的稳定性、收敛性

§6.0 引 言

1． 主要考虑如下的一阶常微分方程初值问题的求解:

?dy?dx?f?x,y????y?x0??y0??

微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。 2. 例如微分方程:

xy'-2y=4x ；初始条件: y(1)=-3。

于是可得一阶常微分方程的初始问题

???y??2y??x4?y(1)??3。

显然函数y(x)=x2-4x满足以上条件,因而是该初始问题 的微分方程的解。

3. 但是,只有一些特殊类型的微分方程问题能

够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。因此，只能依赖于数值方法去获得微分方程的数值解。 4． 微分方程的数值解：设微分方程问题的解y(x)的存在区间是[a,b]，初始点x0=a，将[a,b]进行划分得一系列节点x0 , x1 ,...,xn，其中a= x0< x1<…< xn =b。

y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk的近似值y(xk)，即y≈y(xk)，这样y0 , y1 ,...,yn称为微分方程的数值解。如图所示:

a b

x0

x1

x2

...

§6.1 算法构造的主要途径

1 欧拉公式 1.1 构造的思想：

?dy??f?x,y??dx?y?x0??y0微分方程初值问题: ? 利用差商代替一阶导数，即

y(x1)?y(x0)dy?dxx?x0h，

y(x1)?y(x0)?f(x0,y0)h则 。

于是，可求出y(x1)的近似值y1，

同样地，可利用x1处的微分方程可得：

y2?y1?hf(x1,y1)

一般地，利用在xn处的微分方程可得：

y1?y0?hf(x0,y0)

yn?1?yn?hf(xn,yn),n?0,1,2,...

此式称为欧拉公式。 1.2 几何意义：

对于微分方程y'=2(x+1),其通解是y=(x+1)2+c,是一个曲线族，当给定初值条件y(0)=2，其特解为y=(x+1)2+1。 如图所示：

由y(0)=2，过该曲线上一点(0,2)作曲线的切线，

dy?f?0,2?其斜率dxx?0，

∴切线为： y?2?f?0,2?(x?0), 因此可计算出y1,

如此，可根据：y?y1?f?x1,y1?(x?x1)，…

故欧拉法又称欧拉折线法。

1.3 算例：

?dy2y???2?dxx?y?1???1例：? 解：h=0.2 , xi=1+ih

?2?(?1)??1.0?2???1? y1= y0+hf(x0，y0)=-1+0.2×?

y2= y1+hf(x1，

?2?(?1)??1.2?2???0.9333?y1)=-1+0.2×?

y3= y2+hf(x2，y2)=

?2?(??0.9333)??2???0.8?1.4?-0.9333+0.2×? ?2?(?0.8)??1.6?2???0.6?y4= y3+hf(x3，y3)=0.8+0.2×?

y5= y4+hf(x4，

?2?(?0.6)??1.8?2???0.3333?y4)=0.6+0.2×?

y6= y5+hf(x5，

?2?(?0.3333)??2??0?2? y5)=0.3333+0.2×?精确解为：y=x2-2x

xk 1.2

y(xk) -0.96

-1

yk ek 0.04

1.4 1.6 1.8 2.0 2.2

-0.84 -0.64 -0.36 0 0.44

-0.9333 -0.8 -0.6 -0.3333 0

0.0933 0.16 0.24 0.3333 0.44

可以看出误差随着计算在积累。 1.4 Euler法的特点和误差

迭代格式：

yn?1?yn?hf?xn,yn?n?0,1,2,?,N?1 特点：

(1) 单步方法； (2) 显式格式；

2?h(3) 局部截断误差为??。

局部截断误差：

当yn?y?xn?时，由y?xn?按照欧拉方法计算来的yk?1的

误差称为局部截断误差。即，

y?xn?1??[y?xn??hf?xn,y?xn??]

是局部截断误差。 如：

12y?xn?1??y?xn??hy??xn??hy?????2

欧拉法得：

yn?1?yn?hf?xn,yn?

2?h因此，局部截断误差是??。

2 改进Euler法 2.1方法构造

dy?f(x,y)在微分方程初值问题dx,对其从xk到

xk?1进行定积分得:

将右端的定积分用梯形公式来进行近似计算。

xky(xk?1)?y(xk)??xk?1f(x,y(x))dxhy(xk?1)?y(xk)??f(xk,y(xk))?f(xk?1,y(xk?1))?2

用yk?1和yk来分别代替y(xk?1)和y(xk)得计算格式:

hyk?1?yk?(f(xk,yk)?f(xk?1,yk?1))2

这就是改进欧拉方法。 2.2 显式格式和隐式格式 在欧拉式中

yk?1?yk?hf(xk,yk)

每一步计算已知yk,直接用格式可以计算出yk?1,此类格式称为显式格式。

而在改进欧拉方法中

hyk?1?yk?(f(xk,yk)?f(xk?1,yk?1))2

在每一步计算中yk?1是未知,待求的,未知量在f(x,y)中这是一个方程,如f是非线形或超越函

数,此方程是无法直接解出来(要依靠迭代法才能解出)。这类格式称为隐式格式。 2.3 算例

?y??y?x?例: ?y(0)?2用改进欧拉方法求解h?0.2。

hyk?1?yk?(yk?xk?yk?1?xk?1)2解：解得: hh1?yk?1?2yk?2(xk?xk?1)hh1?1? 22注意：由于f(x,y)?y?x,是线形函数可以从隐

式格式中 解出yk?1。

xy(x)?e?x?1 问题的精确解是

xk y(xk) yk ek 欧拉方法 误差 0.2 2.421403 2.422222 0.000813 0.4 2.891825 2.893827 0.00200 0.02140 0.05183 0.6 3.422119 3.425789 0.00367 … … … … 0.09411 … 2.0 10.38906 10.43878 0.04872 3 预测－校正方法

1.1973 由于改进Euler法是隐式格式，无法从格式中直接求出yn?1必须要解方程。下面用预测－校正方法来求隐式格式中的yn?1。 预测值: yn?1?yn?hf?xn,yn? 校正值:

hyn?1?yn??f?xn,yn??f?xn?1,yn?1???? 2此式相当于对隐式格式求yk?1时采用迭代的方法,用欧拉格式得到的yk?1作为初始值迭代公式迭代

一次而已,此公式代入后得:

hyk?1?yk?(f(xk,yk)?f(xk?h,yk?hf(xk,yk)))2

??yp?yk?hf(xk,yk)??yc?yk?hf(xk?1,yp)??yk?1?1(yp?yc)2如改写成平均的形式为:?

§6.2 龙格-库塔法Runge-Kutta

Method

1 龙格-库塔法的思想

?y'?f(x,y)?1.1 考虑微分方程的初值问题:?y(x0)?y0 根据微分中值定理有：

y(xk?1)?y(xk)?y'(xk??h)h，其中0<θ<1。

于是

y'(xk??h)?f(xk??h,y(xk??h))

即

y(xk?1)?y(xk)?hf(xk??h,y(xk??h))

我们称f(xk??h,y(xk??h))为y(x)在区间[xk, xk+1]上的平均斜率，记作K，其中，θ是存在但是未知的。 因此，如何对平均斜率K进行近似计算，相应地就得到一种近似公式，或称为微分方程的一种计算格式。

1.2 例如： 用f(xk,yk)作为平均斜率K的近似值就得到欧拉格式：

yk?1?yk?hf(xk,yk)。

1(f(xk,y(xk)?f(xk?1,yk?1))用2作为平均斜率K的近似值就得到比欧拉格式高一阶精度的格式，即，改进欧拉格式的预测-校正方法:

1.3 启发(Motivative)

hyk?1?yk?(f(xk,yk)?f(xk?1,yk?1))2

能否在二维平面中x∈[xk, xk+1], y∈[yk, yk+1]上

多找一些f(x,y)点，在这些点上作函数值的平均数并以此作为平均斜率K的近似值。由于自由度

p?1O(h)的p能够提高，从而达到提的增加，使得

高精度的目标，这就是龙格-库塔法的基本思想。 1.4 R-K的一般形式

在?(x,y)|x?[xk, xk?1], y?[yk, yk?1]?二维平面区域上取N个点，得到公式:

，

其中Kj是二元函数f(x,y)在这N个点上的值。

K1?f(xk,yk),K2?f(xk??2h,yk??21hK1),K3?f(xk??3h,yk?h(?31K1??32K2)),??Kj?f(xk??jh,yk?h(?j,1K1????j,j?1Kj?1)) ??其中?j,?j,?j是待定常数。

j?1yk?1?yk?h??jKjN2 二阶龙格-库塔法

?K1?f?xk,yk???K2?f?xk??2h,yk??21hK1??2.1 计算格式: ?yk?1?yk?h??1K1??2K2? 应用Taylor公式求参数：?2,?21,?1,?2。

将f(x,y)在(xk , yk)上展开

K1?f(xk,yk)?y?(xk)K2?f(xk??2h,yk??21hK1)?f(xk,yk)?h?2fx??h?21K1fy??O(h2)

又

12y(xk?1)?y(xk)?hy'(xk)?hy''(xk)?O(h3)2

因

y'?fdy''(xk)?f(xk,y(xk))?fx??fy??y'?fx??fy??f dx 以上两式代入y(xk?1)中，又由yk?1?yk?h(?k1??k2)得

h23??y?hf?(fx?fyf)?O(h)22???y?h??f??(f??hf?h?ff?O(h)?22x21y?1?

比较两边得

??1??2?1?f:?1???fx:?2?2?2?1??fyf:?21?2??2 ?四个未知量,三个方程,有无穷多组解，但局部截

3O(h),是二阶精度的。 断误差均为

1???2?12, 例如: (1) 取?2?1,则?21?1,

就是改进欧拉公式的预测-校正方法(或二阶龙格-库塔法)。

??K1?f?xk,yk???K2?f?xk?h,yk?hK1??h?yk?1?yk?(K1?K2)2? 11?2??21?2，则2，则?1?0,?2?1, (2) 取

得二阶中点法：

?K1?f?xk,yk??11????K2?f?xk?h,yk?hK1?22????yk?1?yk?hK2?

3 三阶龙格－库塔方法

3.1 计算格式: ?K1?f(xk,yk)?K?f(x?ah,y?h?K)?2k2k211??K3?f(xk?a3h,yk?h(?31K1??32K2))??yk?1?yk?h(?1K1??2K2??3K3)

参数有?1,?2,?3,?2,?3,?21,?31,?32共有八个参数。 将f和y?xk?1?均在(xk,yk)点Tayor展开,有 K1?f?xk,yk??y??xk?

12212222???a2?21hffxy????21hffyy???O(h3)K2?f?a2hfx???21hffy??a2hfxx22122??K3?f??3hfx??(?31K1??32K2)hffy???3hfxx2122?????31K1??32K2?fyy?????h3???3??31K1??32K2?hfxy 2利用K1?f,K2按展开式代入(高于h2的不计)得:

K3?f??3hfx??(?31??32)hffy???2?32hfx?fy???21?32hf?fy??22212212222?????????h3??3hfxx??3??31??32?hffxy???31??32?hffyy 22y?xk?1??y?xk??hy??xk??y?(xk)?f

dy???xk??f?x,y?xk???fx??ffy?dx2dd2???2ffxy???ffyy???fx?fy??f?fy??y????xk???y???x????fx??ffy???fxxdxdx 以上展开式代入yk?1?yk?h??1K1??2K2??3K3?中比较得:

h1??y?xk??y????xk????h4? 262??f ?1??2??3?1 1?2?2??3?3?2 1?21?2???31??32??3?2 12121?2?2??3?3?226 (1) (2) (3) (4) fx ffy? ?? fxx?? ffxy1(5) ?2?2?21??3??31??32??3?3 121122?? ?21?2???31??32??3?(6) ffyy226 1fx?fy? (7) ?3?2?32?6 21(8) ?3?21?32?f?fy?? 6 上述八个方程有六个独立方程。

??1??2??3?1???2?2??3?3?12???21??2???31??32??3?122??2?2??3?3??3?1??2?32?6?

八个未知量,六个方程有无穷多组解,

4?(h),具有三阶精度。 但其截断误差均为

y?xk?1??y?xk???xk?1?xk?1xkP?x?dx?5?h截断??部分，用等距步长，上面积分很简单，得到的方法就是显式四阶Adams方法。

hyk?1?yk?[?9fk?3?37fk?2?59fk?1?55fk]24

5?h以上可以看到该方法的局部截断误差是??因

?????xkP?x?dx???h5?xk?1xkxk?1l0?x?dxfk?2k?2xk???l?x?dx?fl?x?dx?f???l?x?dx?fxk?1xk1xk?1xk3

k?1k?3

而是四阶精度的。

2x???y?y?y??y?0??1例如：?

解：取h?0.2，首先用四阶Runge－Kutta

y2?1.3416803，方法来起步，计算出y1?1.1832293，

y3?1.4832838,下面不必用Runge－Kutta方法，而开始用四阶Adams方法。 (1)、求y4

2x02x1f0?y0??1f1?y1??0.8451714y0y1

2x32x2f2?y2??0.745413f3?y3??0.674268y2y3 hy4?y3?[?9f0?37f1?59f2?55f3]?1.611423124

(2)、求y5

2x4f4?y4??0.6185115y4只要补算

hy5?y4?[?9f1?37f2?59f3?55f4]?1.729840324

(3)、求y6

2x5f5?y5??0.5736642y5只要补算 hy6?y5?[?9f2?37f3?59f4?55f5]?1.880661624

现列表看用Adams方法求出的误差，精解为?2x?1 y?xk? xk yk ek 0.8 1.6114231 1.6124515 1.02?10-3 1.0 1.2 1.7298403 1.8406616 1.7320508 1.8439089 2.21?10-3 3.24?10-3 2.3 隐式Adams方法

用xk?2,xk?1,xk,xk?1作为插值结点，由于xk?1也是插值结点，必带来fk?1从而导致是隐式格式。 P?x??l?2?x?fk?2?l?1?x?fk?1?l0?x?fk?l1?x?fk?1 用插值多项式P?x?来代替积分中的f?x,y?x??得：

5?h截掉??得近似公式：

yk?1?yk?A?2fk?2?A?1fk?1?A0fk?A1fk?1

y?xk?1??y?xk???xk?1xkP?x?dx???h5?

Ai??xk?1xkli?x?dx得：

h5199A?2?,A?1??h,A0?h,A1?h24242424，

从而得四阶隐式Adams方法。

hyk?1?yk?(fk?2?5fk?1?19fk?9fk?1)24

因fk?1?f?xk?1,yk?1?，而yk?1是未知的，故这是隐式格式。

隐式格式的解法―用预测-校正法：

用显式格式作为预测值，再用隐式格式来校正。

预测值：

hyk?1?yk?(?9fk?3?37fk?2?59fk?1?55fk)24

校正值：

hyk?1?yk?(fk?2?5fk?1?19fk?9f?xk?1,yk?1?)24

§6.4 线性多步法的一般形式

?dy??f?x,y??dx?y?x0??y0微分方程初值问题: ? 1 k步线性法的一般结构

??yii?0kn?i?h??ifn?ii?0k

?i,?i是常数 ——线性多步法

?k?0 , ?k?1 ——k步法 ?0,?0不同时为零 ——k步法 ?k?0 ,隐式结构;否则为显式结构。

2 k步p阶线性法的推导

?y?由:

i?0kin?i

?h??ifn?ii?0k

定义算子：

'????yx；h??yx?ih?h?y????i?x?ih??????i?i?0k

由Taylor展开右端有

y?x?ih??y?x??ihy?x?'2!2ih?'''?''''y?x?ih??y?x??ihy?x??y?x???2!

ih???2y''?x???代入上式：

'p?p????yx；h?Cyx?Chyx???Chy?x??? p????0??1??其中，

?C0??0??1????k?k?C1???i??i??i??i?0?????11ppp?1p?1C???2????k????2????k?k????p12k?12p!?p?1?!? 因此，若y?x?有p?1次可微，令C0?C1???Cp?0,Cp?1?0,

??y?h?f??则

iii?0kx?ih?Cp?1hp?1y?p?1??x?

?y?x?；h??对应的线性多步法为p阶k因而算子??步方法。

??yii?0kn?i?h??ifn?ii?0k3 算例

例：考虑2步法 k?2,?2?1,记?0?a

?C0?a??1?1?0??C1??1?2???0??1??2??0?1?C2???1?4????1?2?2??02!??11?C3???1?8????1?4?2??03!2!? ??1??1?a,???0??1?1?5a?,12???2??1?3?1?a?????1?5?a?2?12于是： ?

故一般的二步法为：

yn?2??1?a?yn?1?aynh??5?a?fn?2?8?1?a?fn?1??1?5a?fn???? 12同时可以算出：

1C4???1?a?4!1C5??17?13a?? 3?5!当a??1时，C4?0,上述公式为三阶的。

当a??1时，C4?0,C5?0上述公式为四阶的。

□

§6.5 收敛性与稳定性

1 定义[收敛性]：

若某算法对于任意固定的 x = xn= x0 + n h，当 h?0 ( 同时 n ? ?) 时,有 yn? y( xn)，则称该算法是收敛的。 1.2 算例

?y???y?例：对于初值问题?y(0)?y0 考察欧拉显式格式的收敛性。

?xy(x)?ye解：该问题的精确解为 0欧拉公式为 yn?1?yn?hf(xn,yn)

即：

yn?1?yn?h?yn?(1??h)yn

对任意固定的 x = xn = n h ，有

yn?(1??h)yn?1?(1??h)2yn?2???(1??h)ny0?y0(1??h)xn/h?y0[(1??h)1/?h]?xn

?xnlim(1??h)1/?h?e?ye?y(xn) 由h?0 0 □

2 稳定性

x)?y?(x)??30y(??1例：考察初值问题 ?y(0) 在区间[0,

0.5]上的解。分别用欧拉显式、隐式格式和改进的欧拉格式计算数值解。

节点 xi 0.0 0.1 0.2 0.3 0.4 0.5

欧拉显式 1.0000 ?2.0000 4.0000 ?8.0000 1.6000?101 ?3.2000?101

欧拉隐式 1.0000 2.5000?10?1 6.2500?10?2 1.5625?10?2 3.9063?10?3 9.7656?10?4

改进欧拉法 1.0000 2.5000 6.2500 1.5626?101 3.9063?101 9.7656?101

精确解 1.0000 4.9787?10?2 2.4788?10?3 1.2341?10?4 6.1442?10?6 3.0590?10?7

误差包括截断误差(算法理论误差)y?xi??yi和 舍入误差两个部分。后者由计算机字长等决定，属于稳定性问题。

若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的

讨论数值方法的稳定性，通常只考虑试验方程

y???y

当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差 ?0?y0?y 0 ，则若此误差以后逐步衰减，就称该算法相对于h??h绝对稳定，h??h的全体构成绝对稳定区域。有时，我们称算法A 比算法B 稳定，就是指 A 的绝对稳定区域比 B 的绝对稳定区域大。 2.1 算例：考察显式欧拉法算法的稳定性

解：

yn?1?yn?h?yn?(1?h?)yn?(1?h?)n?1y0

由?0?y0?y0 得：

n?1yn?1?(1??h)y0

于是：

?n?1?yn?1?yn?1?(1??h)n?1?0

因此，要保证初始误差?0 以后逐步衰减， 必须满足：|1??h|?1

Im -2 - 1 0 Re 例：考察隐式欧拉法

yn?1?yn?h?yn?1

?1??1???n?1???0yn?1??yn???1?h?? ?1??h?

可见绝对稳定区域为：|1??h|?1

Imn?1 0 1 2 Re

ODE数值方法的小结

?dy??f?x,y??dx?y?x0??y0微分方程初值问题: ? 1 欧拉公式

2?h单步方法；显式格式；局部截断误差为??

yn?1?yn?hf(xn,yn),n?0,1,2,...

2 改进Euler法

3?h?? 单步方法；隐式格式；局部截断误差是

hyn?1?yn?(f(xn,yn)?f(xn?1,yn?1))2

3 预测－校正方法

预测值: yn?1?yn?hf?xn,yn?

hyn?1?yn??f?xn,yn??f?xn?1,yn?1???? 2 校正值: 4 龙格-库塔法

??K1?f?xn,yn???K2?f?xn?h,yn?hK1???yn?1?yn?h(K1?K2)2 ·二阶龙格-库塔法：?

· 二阶中点法： ?K1?f?xn,yn??11????K2?f?xn?h,yn?hK1?22????yn?1?yn?hK2? ·三阶龙格－库塔方法：

?K1?f?xn,yn???K?f?x?1h,y?hK?2n1??n2?2?????K3?f?xn?h,yn?hK1?2hK2??121??yn?1?yn?h?K1?K2?K3??636??? ??K1?f?xn,yn???K?f?x?1h,y?hK??n3n31??2????22??K?fx?h,y?hKn2??3?n33????3??1?yn?1?yn?h?K1?K3?4? ?4·三阶Heun法： ?·四阶龙格－库塔方法：

?K1?f?xn,yn???K?f?x?1h,y?hK?n1??n2?22???hh????K3?f?xn?,yn?K2?22????K4?f?xn?h,yn?hK3??h??yn?1?yn?6?K1?2K2?2K3?K4??

5 多步法及其它

·显式四阶Adams方法：

hyn?1?yn?[?9fn?3?37fn?2?59fn?1?55fn]24

·隐式四阶Adams方法：

hyn?1?yn?(fn?2?5fn?1?19fn?9fn?1)24

·隐式格式的解法―用预测-校正法： 预测值：

hyn?1?yn?[?9fn?3?37fn?2?59fn?1?55fn]24

校正值：

hyn?1?yn?(fn?2?5fn?1?19fn?9f?xn?1,yn?1?)24

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