常微分方程初值问题数值解法

Rn?1?y(xn?1)?[yn?1??ajyn?j?h?bjf(xn?j,yn?j)]

j?0j??1k?1k?1为k步公式(2.7)在xn?1处的局部截断误差. 当

Rn?1?O(hp?1)

时称式(2.7)是p阶的.

应用方程y'(x)?f(x,y(x))可知局部截断误差也可写成为

Rn?1?y(xn?1)?[yn?1??ajyn?j?h?bjy'(xn?j)]

j?0j??1k?1k?1定义2.4 如果线性k步方法(2.7)至少是1阶的, 则称是相容的; 如果线性k步法(2.7)是p阶的, 则称是p阶相容的. 2.4.1 Adams外插法

将微分方程

dx?f(x,y),x?[a,b], dy的两端从xn到xn?1进行积分, 得到

xn?1y(xn?1)?y(xi)??我们用插值多项式代替右端的被积函数.

xnf(x,y(x))dx (2.8)

Adams外插法选取k个点xn,xn?1,?,?,xn?k?1作为插值基点构造f(x,y)的k-1阶多项式Adams外插法的计算公式为

yn?1?yn?h?aj?jfm

j?0k?1其中aj满足如下代数递推式:

111aj?aj?1?aj?2???a0?1,

23j?1j?0,1,?

根据此递推公式, 可逐个的计算aj?j?0,1,??, 表2.1给出了aj的部分数值:

表2.1

j 0 1 1 1/2 2 5/12 3 3/8 4 251/720 5 95/288 ? aj ? 2.4.2 Adams内插法

根据插值理论知道, 插值节点的选择直接影响着插值公式的精度, 同样次数的内插公式的精度要比外插公式的高.

仍假定已按某种方法求得问题(2.1)的解y(x)在xi(i?0,1,?n)处的数值yi, 并选取插值节点xn?k?, p是正整数, 用Lagrange型插值多项式Ln,k?1(x)构造,,xn?p1,?,xn?(p)y'?f(x,y)可以导出解初值问题(2.1)的Adams内插公式:

xn?1y(xn?1)?y(xn)?当p?0时上式就退化为内插公式.

xn?L(p)n,k?1(x)dx (2.9)

内插公式(2.9)除了包含f在xn?k?1,?,xn处的已知值外, 还包含了在点xn,?,xn?p, 处的未知值.因此内插公式(2.9)只给出了未知量yn,?,yn?p的关系式, 实际计算时仍需要

解联立方程组. p?1的内插公式是最适用的, Ln,k?1(x)采用Newton向后插值公式得到Adams内插公式

(1)yn?1?yn?h?a*j?jfm

j?0k其中系数aj定义为

0????a*j?(?1)j???d?

?1?j?*j?0,1,?

其中aj满足如下代数递推式:

*1*1*1*?1,j?0 aj?aj?1?aj?2???a0??0,j?023j?1?*根据此递推公式, 可逐个的计算a*j?j?0,1,??, 表2.2给出了a*j的部分数值:

表2.2

j 0 1 1 -1/2 2 -1/12 3 -1/24 4 -19/720 5 -3/160 ? ? a*j 2.5 算法的收敛性和稳定性 2.5.1 相容性

初值问题(2.1):

dx?f(x,y),x?[a,b], dyy(x0)?y0

的显式单步法的一般形式为

yn?1?yn?h?(xn,yn,h), n?0,1,?,N?1 (2.10)

其中h?b?a, xn?a?nh. 这样我们用差分方程初值问题(2.10)的解yn作为问题(2.1)的N解y(xn)在x?xn处的近似值, 即y(xn)?yn.因此, 只有在问题(2.1)的解使得

y(x?h)?y(x)??(x,y(x),h)

h逼近

y'?f(x,y(x))?0

时, 才有可能使(2.10)的解逼近于问题(2.1)的解. 从而, 我们期望对任一固定的x?[a,b], 都有

lim[h?0y(x?h)?y(x)??(x,y(x),h)]?0

h假设增量函数?(x,y,h)关于h连续, 则有

?(x,y,h)?f(x,y)

定义 2.5 若关系式

?(x,y,h)?f(x,y)

成立, 则称单步法(2.10)与微分方程初值问题(2.1)相容.

容易验证, Euler法与改进Euler法均满足相容性条件. 事实上, 对Euler法, 增量函数为

?(x,y,h)?f(x,y)

自然满足相容性条件, 改进Euler法的增量函数为

?(x,y,h)?[f(x,y)?f(x?h,y?hf(x,y))]

因为f(x,y), 从而有

12?(x,y,h)?[f(x,y)?f(x,y)]?f(x,y)

所以Euler法与改进Euler法均与初值问题(2.1)相容. 一般的, 如果显示单步法有p阶精度(p?0), 则其局部误差为

12y(x?h)?[y(x)?h?(x,y(x),h)]?O(hp?1)

联系客服:779662525#qq.com(#替换为@)