常微分方程初值问题数值解法(卷名：数学)
numerical method for initial value problem of ordinary differential equations
根据给定的初始条件，确定常微分方程惟一解的问题叫常微分方程初值问题。大多数实际问题难以求得解析解，必须将微分问题离散化，用数值方法求其近似解。
一阶常微分方程的初值问题的提法是，求出函数 y(x)，使满足条件
(1)利用数值方法解问题 (1)时，通常假定解存在且惟一，解函数y（x）及右端函数 ƒ(x，y)具有所需的光滑程度。数值解法的基本思想是:先取自变量一系列离散点，把微分问题(1)离散化，求出离散问题的数值解，并以此作为微分问题解y(x)的近似。例如取步长h>0，以h剖分区间[α，b]，令xi=α+ih，把微分方程离散化成一个差分方程。以y(x)表微分方程初值问题的解，以yi表差分问题的解，
就是近似解的误差，称为全局误差。因此，设计各种离散化模型，求出近似解，估计误差以及研究数值方法的稳定性和收敛性等构成了数值解法的基本内容。离散化方法 常用的有三种：
① 基于数值微分的方法 将方程(1)左端的导数用某个一阶数值微分公式代替，例如在xn点以(yn+1-yn)/h代替yń即得到欧拉向前公式
(2)若在xn+1点以(yn+1-yn)/h代替
则得到欧拉向后公式
(3)取(2)、(3)的平均，可导出二阶精度的梯形公式
(4)② 基于泰勒展开的方法 设计一个算法，假定公式中含有某些待定常数，在函数光滑的假定下，将其按泰勒展开并与微分方程解 y(xn+h)的展式中h的同幂次项相比较，按照给定的精度阶得到待定常数应满足的一些方程，通过这些方程确定待定常数，即可得到所要的差分公式。由此法可导出龙格－库塔公式。
设计算公式有下列形式
(5)
αi、βij为待定常数。取定N值，可按上述泰勒展开的方法确定它们。最常用的显式4阶龙格-库塔公式为
(6)式中
③ 基于函数数值积分的方法 将微分方程的解y(x)代入方程(1)，在子区间[xn-ih，xn+jh]上积分得到公式
(7)积分号下是x的函数，若用某些结点上ƒ的值的数值积分公式近似这个积分，便得到各种差分公式，特别地，若取i=0，j=1，并用ƒ在结点xn，xn-1，xn-2，…上的插值代替(7)式中的被积函数，便得亚当斯外推公式。4阶亚当斯外推公式为
(8)若取ƒ在结点xn+1，xn，xn-1，…上的插值代替(7)式的被积函数，则得亚当斯内插公式。4阶亚当斯内插公式为
(9)解法可按计算 yn+1时用多少个结点上的值分为单步法和多步法，又可以按yn+1出现的形式分为显式法和隐式法。
单步法是指已知结点xn上yn的值便可计算yn+1的值的解法，如(2)、(3)、(4)。单步法是可以自己起步的，即可从方程的初值y0一步步算出y1，y2，…的值。
多步法是指已知yn，yn-1，…，yn-k+1(k≥2)的值才能计算yn+1的值的解法，又称k步法。例如，(8)是四步法，(9)是三步法。多步法不能自己起步，即给了初值y0以后，还要用其他解法(如单步法)，算出y1，y2，…，yk-1后，才能使用多步法，继续往下计算。多步法公式若对yi和ƒi都是线性的，则称作线性多步法，k步线性多步法的一般形式为
显式法的公式中，未知的yn+1明显地被表示，即公式中除yn+1一项外，其他的项中不再含有yn+1，如公式(2)。
隐式法的公式不显含yn+1，求未知的yn+1时一般需要解方程，如公式(3)或(4)。通常用各种迭代方法解隐式差分方程，也可采用较简单的预估-校正方法，如使用梯形公式(4)时，可先用显式公式(2)求得yn+1的预估值，代入式(4)的函数ƒn+1中， 再求得yn+1的值。此法又称改进的欧拉折线法。
数值解法满足相容的、收敛的、数值稳定的条件时，才有实用价值。为此要研究以下的一些问题。
相容性 将微分方程离散化所带来的误差叫截断误差。当h→0时，截断误差趋于零，则称离散化后的方程与微分方程具有相容性，表示离散化后的方程是微分方程的近似。若截断误差的主要项为Chp+1，则称截断误差的阶是p+1，而称该解法是p阶的。p越大表示离散化后的方程与微分方程近似程度越高。
收敛性 是指当h→0时，全局误差εi→0，即离散问题的解yn收敛于微分问题的解y(x)，这是离散解可用的理论基础。p阶的解法，即是当h→0时，εi以hp的速度收敛。
误差估计 对全局误差εi的估计，是应用数值解法时最关心的问题。先验估计通常只能给出误差的阶，即误差的主要项中步长h的幂次。一般采用事后估计，即在计算的过程中估计误差，例如用理查森外推法估计误差。外推法也是提高解的精确度的有效方法。
数值稳定性 是指计算过程中，某一步上产生的误差一步一步地传递下去，是衰减、不增或有界，使得传递下来的误差不致于影响数值解的精度，至少是不会湮没数值解。数值稳定性是常微分方程数值积分时必须考虑的问题。
1956年G.达赫尔斯特证明:存在2k阶k步线性多步法，但数值稳定的k步线性多步法，当k为偶数时，其阶不能超过k+2，当k为奇数时，其阶不能超过k+1。称为限制性定理。
判别一个数值方法的稳定性时，微分方程
(10)有较广的代表性，这里λ=α+iβ，α<0，
许多数值稳定性的定义都以这个方程为基础。通常称它为测验方程。
稳定区域 是指将一个数值积分方法应用于测验方程(10)，在λh复平面上使方法数值稳定的区域。欧拉公式（2）的稳定条件为｜1+λh｜<1，其稳定区域是以[-2，0]为直径的圆的内部（图1）。龙格-库塔法(6)的稳定区域由条件
(11)来确定，实际上所有4阶显式龙格-库塔法的稳定区域都是由条件(11)确定，是图2图形的内部。欧拉向后公式 (3)的稳定区域是以[0，2]为直径的圆的外部（图3）。稳定区域是稳定条件的几何表示，其作用在于解线性常系数常微分方程组
(12)时，若步长 h取得使所有λh（λ是A的本征值）都落在稳定区域内，用这个步长积分时数值稳定，理论上稳定区域可以用来选择可用的步长h。由图1、图2可知，对于α<0的λ，只要h取得足够小，λh 就可落在稳定区域内。对于图3，h没有限制。稳定区域包含 λh平面的整个左半平面的方法叫做 A稳定的。如公式(3)、(4)是A稳定的。刚性方程组、常微分方程组的初值问题为：
式中
皆为n维向量。假定微分方程右端函数ƒ的雅可比矩阵
的特征值为 
且 Reλj<0(j=1，2，…，n)，当比值
时，用通常显式公式计算，│hλj│(j=1，2，…，n)不能超过某个量，也就是h必须很小，从而大大增加计算时间。这类问题称为刚性问题。因为 A稳定的方法从数值稳定来说对步长h没有限制，适用于刚性方程组，如欧拉向后公式(3)，梯形公式(4)，隐式龙格-库塔公式等都是有效的方法。参考书目
P. Henrici，Discrete variable Methods in Ordinary Differential Equations， John Wiley & Sons， New York， 1962.
C. W.吉尔著，费景高、刘德贵、高永春译：《常微分方程初值问题的数值解法》，科学出版社，北京，1978。(C. W. Gear，Numerical Initial Value Problems in Ordinary Differential Equations，Prentice-Hall， Englewood Cliffs， New Jersey，1971.)