常微分方程数值方法

微分方程$y^{\prime }=y,y(0)=1$的数值积分求解图。 Template:Legend Template:Legend Template:Legend 步长$h=1.0$。

$h=0.25$时的同一图。$h\to 0$时，中点法比欧拉法收敛得更快。

常微分方程数值方法是用以寻找常微分方程（ODE）解的数值近似值的方法。其使用也称作“数值积分”，不過「数值积分」主要是指积分的计算。

很多微分方程无法精确求解。但在工程学等领域的实际应用中，通常只需得到数值近似解。本文介绍的算法可用于计算这种近似值，另一种方法是用微积分技术得到解的级数展开表达。

常微分方程出现于物理学、化学、生物学、经济学等学科中。此外，偏微分方程数值方法中的一部分将偏微分方程转为常微分方程，然后可用本文所述方法求解。

问题

一阶微分方程是有如下形式的初值问题（IVP）：Template:Rp Template:NumBlk 其中f是函数$f:[t_0,\mathrm{\infty })\times {ℝ}^d\to {ℝ}^d$，初始条件$y_0\in {ℝ}^d$是已知向量。“一阶”是说方程中只出现了y的一阶导。

在不失高阶系统一般性的前提下，限制(l)式為一阶微分方程，因为高阶ODE可通过引入额外变量转换为更大的一阶方程组。例如，二阶方程$y^{″}=-y$可重写为2个一阶方程：$y^{\prime }=z,z^{\prime }=-y$。

本文介绍IVP的数值方法，并指出边值问题（BVP）需要一套不同的工具：需要在多个点上定义解y的值或成分，因此要用不同方法求解BVP，如打靶法（及其变体），或差分、伽辽金法、Template:Le之类全局方法，都适用于此类问题。

Picard–Lindelöf定理指出，只要f是利普希茨连续的，就有唯一解。

方法

解一阶IVP的数值方法可分为两大类：Template:Le与龙格-库塔法。还可进一步划为显式或隐式，例如隐式线性多步法包括亚当斯-莫尔顿法（Adams-Moulton methods）与Template:Le（BDF），隐式龙格-库塔法Template:Rp包括对角隐式龙格-库塔法（DIRK）、单对角隐式龙格-库塔法（SDIRK）与基于高斯求积的高斯–拉道法等等。线性多步法中的显式方法有亚当斯-巴什福思法，布彻表（Butcher tableau）为下对角的龙格-库塔法都是显式方法。根据经验，刚性微分方程需要用隐式方案，而非刚性问题则可用显示方案更有效地求解。

所谓Template:Le（GLM）是上述两大类方法的概括。

欧拉法

Template:Details 从曲线上任意一点出发，沿与曲线相切的直线移动一小段距离，就能得到曲线上邻近点的近似值。

从微分方程(Template:EquationNote)开始，可用差分近似代替导数y′： Template:NumBlk 重新排列后得到以下公式

$y(t+h)\approx y(t)+h{y}^{\prime }(t)$

利用(Template:EquationNote)可得 Template:NumBlk 此式的应用通常如下：择步长h，构造序列$t_0,t_1=t_0+h,t_2=t_0+2h,...$记精确解$y(t_n)$的数值估计值为$y_n$。受(Template:EquationNote)启发，可用下面的递归方法计算估计值： Template:NumBlk 这就是（前向）欧拉法，是莱昂哈德·欧拉（1768）描述的方法。

欧拉法是显式方法的一个例子，这是说新值$y_{n+1}$是根据已知值（如$y_n$）定义的。

反向欧拉法

Template:Details 若不用(Template:EquationNote)，而用近似值 Template:NumBlk 则得到反向欧拉法： Template:NumBlk 反向欧拉法是隐式方法，这是说需要求解一个方程才能得到新值$y_{n+1}$。通常用定点迭代或牛顿-拉弗森法（的某种修改版）实现之。

隐式方法求解这方程比显示方法直接代入要花更多时间，选择方法时必须考虑这一成本。隐式方法（如(Template:EquationNote)）的优点是求解刚性方程时通常更稳定，便可以使用更大的步长h。

一阶指数积分法

Template:Details Template:Le描述了一大类积分器，近来得到了长足发展。它们至少可以追溯到20世纪60年代。

此處不用(Template:EquationNote)，假设微分方程形式为 Template:NumBlk 或已被局部线性化，围绕背景状态产生线性项$-Ay$与非线性项$𝒩(y)$。

将(Template:EquationNote)乘以$e^{At}$，并在时间区间$[t_n,t_{n+1}=t_n+h]$内精确积分，便构造得到了指数积分：

$y_{n+1}=e^{-Ah}{y}_n+\underset{0}{\overset{h}{\int }}{e}^{-(h-\tau )A}𝒩\left(y(t_n+\tau )\right)d\tau .$

此积分方程是精确的，但并没有定义积分。

使$𝒩(y(t_n+\tau ))$在整个区间内不变，可实现一阶指数积分： Template:NumBlk

推广

欧拉法往往不够精确。更准确地说，只有一阶（下面将介绍阶的概念），这就促使数学家寻找高阶方法。

一种方法是，用$y_n$以及更多之前的值确定$y_{n+1}$，所谓多步法便实现了这种想法。最简单的可能是蛙跳积分法，其是二阶方法，（粗略地说）依赖于2个时间值。

几乎所有实用的多步法都属于Template:Le，形式为 $$\begin{array}{cc}& {\alpha }_k{y}_{n+k}+{\alpha }_{k-1}{y}_{n+k-1}+\cdots +{\alpha }_0{y}_n\\ & =h[{\beta }_{k}f(t_{n+k},y_{n+k})+{\beta }_{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +{\beta }_{0}f(t_n,y_n)].\end{array}$$

另一种方法是在区间$[t_n,t_{n+1}]$内取更多点。这产生了得名于卡尔·龙格与马丁·威廉·库塔的龙格-库塔法，其中一种4阶方法尤为流行。

高级特征

要用这些方法求解ODE，需要的不仅是时间步长公式。始终相同的步长效率不高，于是开发了可变步长方法。通常，步长的选择应使每步的（局部）误差低于某容差水平，这意味着方法还要计算误差指标（error indicator），即对局部误差的估计。

这一思想的延伸是在不同阶方法之间动态选择（称作可变阶数方法）。基于理查德森外推法的Bulirsch–Stoer算法之类方法通常用于构建各种不同阶的方法。

其他理想特征还有：

• 输出稠密：不仅对点$t_0,t_1,t_2,\dots$，还对整个积分区间进行低成本数值逼近；
• 事件定位：找到事件发生的位置，例如某函数取0的时间。通常用求根算法；
• 支持并行计算；
• 用于时间积分时，具有时间可逆性。

其他方法

很多方法不在讨论的框架内，如

• 多阶导数法，不仅使用f，还用其导数。这类方法包括Hermite–Obreschkoff法、Template:Le与递归计算解y的泰勒级数系数的Template:Le、Bychkov–Scherbakov法等。
• 二阶常微分方程组法。前面说过，所有高阶常微分方程都可化为形式如(1)的一阶常微分方程组，这固然没错，但未必是最佳方法。Template:Le可直接用于二阶方程。
• 几何积分法专门针对几类特殊的常微分方程（如解哈密顿力学方程的Template:Le法），这些方法在数值求解时会尊重这些类的底结构或几何结构。
• Template:Le是基于状态量化思想的一系列ODE积分方法，在模拟频繁不连续的稀疏系统时非常有效。

实时并行方法

有些IVP要求积分具有很高的时间分辨率和/或很长的时间区间，传统的序列时间步进法无法实时计算（如数值天气预报、等离子体建模与分子动力学中的IVP）。针对这些问题，人们开发了实时并行（PinT）法以便用并行计算缩短运行时间。

早期的PinT法（最早提出于20世纪60年代）最初被研究人员忽视，因为其所需的并行计算架构尚未普及。2000年代初，随着算力的提高，灵活易用的PinT算法——Template:Le重新吸引了兴趣，它适用于求解各种IVP。Template:Le的出现使PinT算法获得更多关注，并启动了能用于世界上最强大的超级计算机的算法开发。截至2023年，最流行的方法有Parareal、PFASST、ParaDiag、MGRIT等，

分析

数值分析不仅包括数值方法的设计，还包括其分析。分析中的3个核心概念是：

• 收敛性：方法是否逼近解；
• 阶数：近似解的程度；
• 数值稳定性：误差是否能得到抑制。

收敛性

Template:Main 若数值解随着步长h趋近于0而逼近精确解，则称此数值方法具有收敛性（convergent）。更确切地说，要求对利普希茨函数f与每个$t^{*}>0$，ODE (1)

$\underset{h\to 0^{+}}{\lim }\underset{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }{\max }\Vert y_{n,h}-y(t_n)\Vert =0.$

上述所有方法都收敛。

一致性与阶数

Template:Details 设数值方法

$y_{n+k}=\mathrm{\Psi }(t_{n+k};y_n,y_{n+1},\dots ,y_{n+k-1};h).$

方法的局部（截断）误差定义为迭代一步产生的误差。即，假设之前的迭代无误差，则是此方法结果与精确解之间的差

${\delta }_{n+k}^h=\mathrm{\Psi }(t_{n+k};y(t_n),y(t_{n+1}),\dots ,y(t_{n+k-1});h)-y(t_{n+k}).$

若

$\underset{h\to 0}{\lim }\frac{{\delta }_{n+k}^h}{h}=0.$

则称此方法一致（consistent）。若

${\delta }_{n+k}^h=O(h^{p+1})\text{as }h\to 0.$

则称此方法阶数为p。因此，阶数不为0的方法是一致的。上文介绍的两种欧拉法(4、6)都是1阶的，因此是一致的。实践中使用的大多数方法阶数更高。一致性是收敛的必要条件Template:Citation needed，但不是充分条件；方法要收敛，必须同时具有一致性与零稳性（zero-stable）。

一个相关概念是全局（截断）误差，即达到固定时间t所需所有步骤中持续存在的误差。明确地说，t时刻的全局误差是$y_N-y(t)$，其中$N=(t-t_0)/h$。第p阶一步法是收敛的，其全局误差是$O(h^p)$。对多阶方法，这一说法不一定成立。

稳定性与刚性

Template:Details 对某些微分方程，欧拉法、显式龙格-库塔法、多步法（如亚当斯-巴什福思法）之类标准方法会表现出解的不稳定性，其他方法则可能产生稳定的解。方程中的这种“困难行为”（本身不一定复杂）称作“刚性”（stiffness），通常是由于底问题中存在不同时间尺度造成的。例如，机械系统中的碰撞（如碰撞振子中的）发生的时间尺度通常比物体运动的时间尺度小得多，这种差异使状态参数曲线变得非常陡峭。

刚性问题在化学动力学、控制论、固体力学、天气预报、生物学、等离子体物理、电子学等领域中无处不在。克服刚性的一种方法是将微分方程推广到微分包含式，从而允许非光滑性并建立其模型。

历史

下面是该领域一些重要进展的时间线。

• 1768 - 莱昂哈德·欧拉发现欧拉法。
• 1824 - 奥古斯丁-路易·柯西证明了欧拉法的收敛性。此证明中，柯西使用了隐性欧拉法。
• 1855 - Francis Bashforth在一封信中首次提到约翰·柯西·亚当斯的多步法。
• 1895 - 卡尔·龙格发现第一个龙格-库塔法。
• 1901 - 马丁·威廉·库塔描述了现在流行的4阶龙格-库塔法。
• 1910 - 路易斯·弗莱·理查德森公布了他的外推法——理查德森外推法。
• 1952 - Charles F. Curtiss与Joseph Oakland Hirschfelder提出刚性方程概念。
• 1963 - Germund Dahlquist引入积分法的A稳定性。

二阶一维边值问题的数值解法

边值问题（BVPs）通常要离散化，得到近似相等的矩阵问题再数值求解。最常用的一维BVP数值求解方法称作有限差分法。这种方法用点值的线性组合构造描述函数导数的有限差分系数，例如一阶导数的二阶中心差分近似为：

$\frac{u_{i+1}-u_{i-1}}{2h}=u^{\prime }(x_i)+𝒪(h^2),$

二阶导数的二阶中心差分近似为：

$\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}=u^{″}(x_i)+𝒪(h^2).$

两式中，$h=x_i-x_{i-1}$是离散域上相邻x值间的距离。这样就构建了线性系统，然后可用标准矩阵法求解。例如，待解方程：

$\begin{array}{cc}& \frac{d^{2}u}{d{x}^2}-u=0,\\ & u(0)=0,\\ & u(1)=1.\end{array}$

下一步是将问题离散化，用线性导数近似，如

$u^{\prime }{\text{'}}_i=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$

并解所得线性方程组。将有如下方程

$\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}-u_i=0,\mathrm{\forall }i=1,2,3,...,n-1.$

初看之下，这个方程组似乎有困难，因为方程中没有不与变量相乘的项，但事实上这是错误的。i = 1或n − 1时，有一项涉及边值$u(0)=u_0$、$u(1)=u_n$，由于两值已知，可以简单地代入，就得到了具有非平凡解的非齐次线性方程组。

相關條目

• Template:Le
• Template:Le
• Template:Le
• Modelica语言及Template:Le软件

参考文献

• Template:Cite book
• J. C. Butcher, Numerical methods for ordinary differential equations, Template:Isbn
• Ernst Hairer, Syvert Paul Nørsett and Gerhard Wanner, Solving ordinary differential equations I: Nonstiff problems, second edition, Springer Verlag, Berlin, 1993. Template:Isbn.
• Ernst Hairer and Gerhard Wanner, Solving ordinary differential equations II: Stiff and differential-algebraic problems, second edition, Springer Verlag, Berlin, 1996. Template:Isbn.
  (This two-volume monograph systematically covers all aspects of the field.)
• Template:Cite journal
• Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996. Template:Isbn (hardback), Template:Isbn (paperback).
  (Textbook, targeting advanced undergraduate and postgraduate students in mathematics, which also discusses numerical partial differential equations.)
• John Denholm Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991. Template:Isbn.
  (Textbook, slightly more demanding than the book by Iserles.)

外部链接

• Joseph W. Rudmin, Application of the Parker–Sochacki Method to Celestial Mechanics Template:Webarchive, 1998.
• Dominique Tournès, L'intégration approchée des équations différentielles ordinaires (1671-1914), thèse de doctorat de l'université Paris 7 - Denis Diderot, juin 1996. Réimp. Villeneuve d'Ascq : Presses universitaires du Septentrion, 1997, 468 p. (Extensive online material on ODE numerical analysis history, for English-language material on the history of ODE numerical analysis, see, for example, the paper books by Chabert and Goldstine quoted by him.)
• Template:Cite journal
• Template:GitHub (C++ library with rigorous ODE solvers)
• INTLAB Template:Wayback (A library made by MATLAB/GNU Octave which includes rigorous ODE solvers)