使用多步法精确求解微分方程的二阶数值积分算法

另一种求解式（5.1）的方法是用一个k次多项式来逼近非线性函数x（t），

式中，系数α₀，α₁，…，α_k为常数。可以证明，任何函数都能在有限区间[t₀，tN]上用一个足够高次的多项式进行逼近（误差小于事先给定的ε）。引入多步法后，式（5.1）的求解就与多项式逼近联系起来了。在多步法中，x_n+1依赖于前面若干节点处的x_n，x_n-1，…及相对应的f（x_n，t_n），f（x_n-1，t_n-1），…；而单步法（比如Runge-Kutta法）只依赖于前一步的信息。一般地，

为了将数值积分算法与多项式逼近联系起来，必须建立两套系数之间的关系。一个k次多项式由k+1个系数（α₀，…，α_k）唯一确定。而上述数值积分算法有2p+3个系数，因此p和k之间必须满足如下关系：

2p+3≥k+1 （5.9）

数值积分算法的阶数等于以t为变量的多项式的最高次数k，对此多项式，数值解与精确解完全重合。这些系数可以通过选择一系列基函数[ϕ₁（t）ϕ₂（t）…ϕ_k（t）]来确定，基函数的表达式为

ϕ_j（t）=t^jj=0，1，…，k

将这些基函数代入到多步法式（5.8）中，得到

式中，j=0，1，…，k。

用上述方法可以推导出几个一阶数值积分算法。考虑p=0和k=1的情形，这种情形满足式（5.9）的限制，因此可以使用上述方法来确定多步法的系数，使得对次数为1的多项式，多步法公式是精确成立的。当k=1时，基函数集合为

ϕ₀（t）=1 （5.10）

ϕ₁（t）=t （5.11）

基函数的导数为

多步法方程为

x_n+1=a₀x_n+b_-1h_n+1f（x_n+1，t_n+1）+b₀h_n+1f（x_n，t_n） （5.14）

将基函数代入多步法式（5.14），得到如下两个方程

将式（5.10）和式（5.11）代入式（5.15）和式（5.16），得到

1=a₀（1）+b_-1h_n+1（0）+b₀h_n+1（0） （5.17）

t_n+1=a0tn+b_-1hn+1（1）+b0h_n+1（1） （5.18）

由式（5.17），可得a₀=1。由t_n+1-t_n=h_n+1和式（5.18），可得

b_-1+b₀=1 （5.19）

阶数p和次数k的选择导致了具有3个未知数的2个方程，因此，其中一个未知数可以取任意值。当选择a₀=1，b_-1=0，b₀=1时，可以得到Euler法：

x_n+1=x_n+h_n+1f（x_n，t_n）

而当选择a₀=1，b_-1=1，b₀=0时，可以得到另一种数值积分算法：

x_n+1=x_n+h_n+1f（x_n+1，t_n+1） （5.20）

这种特殊的数值积分算法通常被称为后退Euler法。注意，在这种算法中，系数b_-1不为零，因此x_n+1的表达式隐式地依赖于函数f（x_n+1，t_n+1）。对于b_-1≠0的数值积分算法，通常被称为隐式算法，否则就称为显式算法。因为f（x_n+1，t_n+1）隐式地（且通常是非线性地）依赖于x_n+1，所以隐式算法一般需要在每个时间节点上进行迭代求解。

现在考察p=0，k=2的情形，此时2p+3=k+1，因此所有系数可以被唯一确定。采用与前面一样的基函数取法，且取ϕ₂（t）=t²，，可得如下3个方程：

1=a₀（1）+b_-1h_n+1（0）+b₀h_n+1（0） （5.21）

t_n+1=a₀tn^+b_-1h_n+1（1）+b₀h_n+1（1） （5.22）

t²_n+1=a₀t²_n+h_n+1（b_-1（2t_n+1）+b₀（2t_n）） （5.23）

如果t_n=0，那么有t_n+1=h_n+1。由式（5.21）～式（5.23），可得a₀=1，，，故有

这个二阶数值积分算法被称为梯形法，它也是一种隐式算法。上述公式之所以被称为梯形法，是因为式（5.24）的右边第二项可以被理解为一个梯形的面积。由于在计算x_n+1时使用了t_n和t_n+1时的信息，因此梯形法可以被看作为两步法。

例5.1采用不同的固定步长，基于Euler法、后退Euler法、梯形法和Runge-Kutta法数值求解如下微分方程。

解5.1这个二阶微分方程必须首先转化为常微分方程组的形式，令x₁=x，，有

通过观察，可知该方程组的解析解为

x₁（t）=cost （5.28）

x₂（t）=-sint （5.29）通常难以找到常微分方程的精确解，但在本例中，该精确解将被用来与数值解作比较。(www.xing528.com)

（1）向前Euler法

使用向前Euler法解上述常微分方程组，可得

x_1，n+1=x_1，n+hf₁（x_1，n，x_2，n） （5.30）

=x_1，n+hx_2，n （5.31）

x_2，n+1=x_2，n^+hf₂（x_1，n，x_2，n） （5.32）

=x_2，n-hx_1，n （5.33）

用矩阵形式表示

（2）后退Euler法

使用后退Euler法解上述常微分方程组，可得

x_1，n+1=x_1，n+hf₁（x_1，n+1，x_2，n+1） （5.35）

=x_1，n+hx₂，n+1 （5.36）

x_2，n+1=x_2，n+hf₂（x_1，n+1，x_2，n+1） （5.37）

=x_2，n-hx_1，n+1 （5.38）

用矩阵形式表示

在求解式（5.39）时，矩阵的逆实际上不是显式给出的，而是采用LU分解的方法来求解此方程的。

（3）梯形法

使用梯形法解上述常微分方程组，可得

用矩阵形式表示

（4）Runge-Kutta法

使用Runge-Kutta法解上述常微分方程组，可得

k₁₁=x_2，nk₂₁=-x1，n

k₁₄=x_2，n+hk₁₃k₂₄=-x_1，n-hk₂₃

及

图5.1 例5.1的各种数值解

采用不同的数值积分算法求解式（5.25）的结果如图5.1所示，图中也给出了精确解cost。注意，梯形法和Runge-Kutta法所得的结果与精确解几乎不可区分。由于向前Euler法和后退Euler法是一阶算法，其精度不如高阶的梯形法和Runge-Kutta法。注意，向前Euler法所得的结果比精确解的幅值稍大，且幅值随着时间的推移越来越大。相反地，后退Euler法所得的结果比精确解的幅值稍小，且幅值随着时间的推移越来越小。两者都是由算法的局部截断误差引起的。向前Euler法倾向于产生随时间增加的数值解（欠阻尼），而采用后退Euler法得到的数值解倾向于过阻尼。因此，在使用这些一阶数值积分算法时必须谨慎。

图5.2给出上述几种数值积分算法的全局误差随时间变化的曲线。注意，向前Euler法和后退Euler法其误差具有相同的幅值但符号相反，这个特性将在本章后面做进一步讨论。放大后的梯形算法和Runge-Kutta算法误差曲线重新画于图5.3中，尽管梯形法是一个二阶的多项式逼近算法，而Runge-Kutta法是一个四阶的Taylor级数展开算法，两者之间的误差仍然是可比的。5.3节将进一步探讨各种数值积分算法误差评估表达式的推导问题。

图5.2 例5.1各种数值解的误差

图5.3 例5.1梯形法和Runge-Kutta法的误差

当使用诸如梯形法的隐式算法来求解非线性微分方程组时，每个时步上都必须采用迭代算法进行求解。例如，考察如下的非线性微分方程组：

采用梯形法对上述方程进行数值积分，得到如下的离散化方程：

由于此非线性表达式中隐含了x_n+1，必须用数值算法进行求解：

式中，k是Newton-Raphson迭代指数，I是单位矩阵，x_n是上一步所得的收敛值。