常微分方程初值问题数值解法 下载本文

(1) 单步法-----计算y(x)在x?xn?1处的值仅取决于xn处的应变量及其导数值.

(2) 多步法-----计算y(x)在x?xn?1处的值需要应变量及其导数在xn?1之前的多个网个节点出的值.

2.2 Euler方法 2.2.1 Euler公式

若将函数y(x)在点xn处的导数y'(xn)用两点式代替, 即

y'(xn)?y(xn?1)?y(xn),

h再用yn近似地代替y(xn), 则初值问题(2.1)变为

?yn?1?yn?hf(xn,yn) (2.2) ??y0?y(x0),n?0,1,2,?(2.2)式就是著名的欧拉(Euler)公式.

2.2.2 梯形公式

欧拉法形式简单精度低, 为了提高精度, 对方程y?f(x,y)的两端在区间

'[xn,xn?1]上积分得

y(xn?1)?y(xn)??改用梯形方法计算其积分项, 即

xn?1xnf[x,y(x)]dx (2.3)

?xn?1xnf[x,y(x)]dx?xn?1?xn[f(xn,y(xn))?f(xn?1,y(xn?1))] 2代入式(2.3), 并用yn近似代替式中y(xn)即可得到梯形公式

hyn?1?yn?[f(xn,yn)?f(xn?1,yn?1)] (2.4)

2由于数值积分的梯形公式比矩形公式精度高, 所以梯形公式(2.4)比欧拉公式(2.2)的精度高一个数值方法.

式(2.4)的右端含有未知的yn?1, 它是关于yn?1的函数方程, 这类方法称为隐式方法.

2.2.3 改进的欧拉公式

梯形公式实际计算时要进行多次迭代, 因而计算量较大. 在实用上, 对于梯形公式(2.4)只迭代一次, 即先用欧拉公式算出yn?1的预估值yn?1, 再用梯形公式(2.4)进行一次迭代得到校正值yn?1, 即

?yn?1?yn?hf(xn,yn)? (2.5) ?h?yn?1?yn?[f(xn,yn)?f(xn?1,yn?1)]?22.2.4 欧拉法的局部截断误差

衡量求解公式好坏的一个主要标准是求解公式的精度, 因此引入局部截断误差和阶数概念.

定义2.1 在yn准确的前提下, 即yn?y(xn)时, 用数值方法计算yn?1的误差

Rn?1?y(xn?1)?yn?1, 称为该数值方法计算yn?1时的局部截断误差.

对于欧拉公式, 假定yn?y(xn), 则有

yn?1?y(xn)?h[f(xn,y(xn))]?y(xn)?hy'(xn)

而将真解y(x)在xn处按二阶泰勒展开式有

h2'''yn?1?y(xn)?hy(xn)?y(?),??(xn,xn?1)

2!'因此有

h2''y(xn?1)?yn?1?y(?)

2!定义2.2 若数值方法的局部截断误差为O(hp?1), 则称这种数值方法的阶数是P. 步

长(h<1)越小, P越高, 则局部截断误差越小, 计算精度越高.

2.3 Runge-Kutta方法

2.3.1 Runge-Kutta方法的基本思想

欧拉公式可改写成

?yn?1?yn?hK1 ?K?f(x,y)nn?1则yn?1的表达式与y(xn?1)的泰勒展开式的前两项完全相同, 即局部截断误差为O(h2).

改进的欧拉公式又可改写成

h?y?y?(K1?K2)n?n?12?. K1?f(xn,yn)??K?f(x,y?hk)n?1n1?2?上述两组公式在形式上有一个共同点: 都是用f(x,y)在某些点上值的线形组合得出

y(xn?1)的近似值yn?1, 而且增加计算f(x,y)的次数, 可提高截断误差的阶. 如欧拉公式,

每步计算一次f(x,y)的值, 为一阶方法. 改进的欧拉公式需计算两次f(x,y)的值, 它是二阶方法. 它的局部截断误差为O(h3).

于是可考虑用函数f(x,y)在若干点上的函数值的线形组合来构造近似公式, 构造时要求近似公式在(xn,yn)处的泰勒展开式与解y(x)在xn处的泰勒展开式的前面几项重合, 从而使近似公式达到所需要的阶数. 既避免求偏导, 又提高了计算方法精度的阶数. 或者说, 在[xn,xn?1]这一步内多预报几个点的斜率值, 然后将其加权平均作为平均斜率. 则可构造出更高精度的计算格式, 这就是Runge-Kutta方法的基本思想.

2.3.2 Runge-Kutta方法的构造

一般地, Runge-Kutta方法设近似公式为 (下面的公式修改了)

p?yn?1?yn?h?ciKi?i?1??K1?f(xn,yn) ?i?2,3,?,p? ??i?1?Ki?f(xn?aih,yn?h?bijKj)?j?1?其中ai, bij, ci都是参数, 确定它们的原则是使近似公式在ym处的泰勒展开式与y(x)在

xn处的泰勒展开式的前面的项尽可能多的重合, 这样就使近似公式有尽可能高的精度. 以

此我们可以通过一个复杂的计算过程得出常用的的三阶和四阶Runge-Kutta公式

h?y?y?(K1?4K2?K3)n?n?16?K1?f(xn,yn)? ??K?f(x?h,y?hK)2nn1?22?K?f(x?h,y?hK?2hK)nn12?3和

h?y?y?(K1?2K2?2K3?K4)n?n?16?K1?f(xn,yn)??hh? (2.6) K2?f(xn?,yn?K1)?22??hhK?f(x?,y?K2)3nn?22??K4?f(xn?h,yn?hK3)?式(2.6)称为经典Runge-Kutta方法.

2.4 线性多步法

在逐步推进的求解过程中, 计算yn?1之前事实上已经求出了一系列的近似值y0, y1,

?, yn. 如果充分利用前面多步的信息来预测yn?1, 则可期望获得较高的精度, 这就是

构造多步法的基本思想.

线性k步方法的一般公式为

yn?1??ajyn?j?h?bjf(xn?j,yn?j) (2.7)

j?0j??1k?1k?1其中aj, bj均为与n无关的常数, ak?1?bk?1?0. 当b?1?0时为显格式; 当b?1?0时为隐格式. 特别当k?1,a0?b0?1,b?1?0时为Euler公式;当k?1,a0?1,b0?b?1?为梯形公式.

1时2定义 2.3 称