茶年 · 中文维基百科

常微分方程数值方法

Wikipedia · 2026-08-01 维基百科 👁 2 阅读

常微分方程数值方法是用以寻找常微分方程（ODE）解的数值近似值的方法。其使用也称作“数值积分”，不過数值积分一般是指用數值方法計算定積分，不是求解常微分方程。

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

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

问题
一阶微分方程是有如下形式的初值问题（IVP）：

其中f是函数f:[t_0, \infty) \times \R^d \to \R^d，初始条件y_0 \in \R^d 是已知向量。“一阶”是说方程中只出现了y的一阶导。

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

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

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

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

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

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

从微分方程()开始，可用差分近似代替导数y′：
{{NumBlk|:| y'(t) \approx \frac{y(t+h) - y(t)}{h}, |}}
重新排列后得到以下公式
: y(t+h) \approx y(t) + hy'(t)
利用()可得

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

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

反向欧拉法
若不用()，而用近似值
{{NumBlk|:| y'(t) \approx \frac{y(t) - y(t-h)}{h}, |}}
则得到反向欧拉法：
{{NumBlk|:| y_{n+1} = y_n + hf(t_{n+1},y_{n+1}). |}}
反向欧拉法是隐式方法，这是说需要求解一个方程才能得到新值y_{n+1}。通常用定点迭代或牛顿-拉弗森法（的某种修改版）实现之。

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

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

此處不用()，假设微分方程形式为
{{NumBlk|:|y'(t) = -A\, y+ \mathcal{N}(y), |}}
或已被局部线性化，围绕背景状态产生线性项-Ay与非线性项\mathcal{N}(y)。

将()乘以e^{A t}，并在时间区间[t_n, t_{n+1} = t_n + h]内精确积分，便构造得到了指数积分：
: y_{n+1} = e^{-A h } y_n + \int_{0}^{h} e^{ -(h-\tau) A } \mathcal{N}\left( y\left( t_n+\tau \right) \right)\, d\tau.
此积分方程是精确的，但并没有定义积分。

使\mathcal{N}( y( t_n+\tau ) )在整个区间内不变，可实现一阶指数积分：
{{NumBlk|:|y_{n+1} = e^{-Ah}y_n + A^{-1}(1-e^{-Ah}) \mathcal{N}( y( t_n ) )\ .|}}

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

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

几乎所有实用的多步法都属于线性多步法，形式为
\begin{align}
&{} \alpha_k y_{n+k} + \alpha_{k-1} y_{n+k-1} + \cdots + \alpha_0 y_n \\
&{} \quad = h \left[ \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) \right].
\end{align}

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

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

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

其他理想特征还有：

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

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

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

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

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

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

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

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

: \lim_{h\to0^+} \max_{n=0,1,\dots,\lfloor t^*/h\rfloor} \left\| y_{n,h} - y(t_n) \right\| = 0.

上述所有方法都收敛。

一致性与阶数
设数值方法

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

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

: \delta^h_{n+k} = \Psi \left( t_{n+k}; y(t_n), y(t_{n+1}), \dots, y(t_{n+k-1}); h \right) - y(t_{n+k}).
若
: \lim_{h\to 0} \frac{\delta^h_{n+k}}{h} = 0.
则称此方法一致（consistent）。若
: \delta^h_{n+k} = O(h^{p+1}) \quad\mbox{as } h\to0.
则称此方法阶数为p。因此，阶数不为0的方法是一致的。上文介绍的两种欧拉法(4、6)都是1阶的，因此是一致的。实践中使用的大多数方法阶数更高。一致性是收敛的必要条件，但不是充分条件；方法要收敛，必须同时具有一致性与零稳性（zero-stable）。

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

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

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

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

  • 1768 - 莱昂哈德·欧拉发现欧拉法。
  • 1824 - 奥古斯丁-路易·柯西证明了欧拉法的收敛性。此证明中，柯西使用了隐性欧拉法。
  • 1855 - 弗朗西斯·巴什福斯在一封信中首次提到约翰·柯西·亚当斯的多步法。
  • 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'(x_i) + \mathcal{O}(h^2),

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

: \frac{u_{i+1}- 2 u_i + u_{i-1}}{h^2} = u*(x_i) + \mathcal{O}(h^2).

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

: \begin{align}
&{} \frac{d^2 u}{dx^2} -u =0,\\
&{} u(0)=0, \\
&{} u(1)=1.
\end{align}

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

: u*_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, \quad \forall i={1,2,3,...,n-1}.

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

相關條目
*
*
*
*Modelica语言及软件

注释
参考文献
*
J. C. Butcher, Numerical methods for ordinary differential equations*,
Ernst Hairer, Syvert Paul Nørsett and Gerhard Wanner, Solving ordinary differential equations I: Nonstiff problems,* second edition, Springer Verlag, Berlin, 1993. .
Ernst Hairer and Gerhard Wanner, Solving ordinary differential equations II: Stiff and differential-algebraic problems,* second edition, Springer Verlag, Berlin, 1996. .
(This two-volume monograph systematically covers all aspects of the field.)
*
Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations,* Cambridge University Press, 1996. (hardback), (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. .
(Textbook, slightly more demanding than the book by Iserles.)

外部链接

*

  • (C++ library with rigorous ODE solvers)

*[http://www.ti3.tu-harburg.de/intlab/ INTLAB] (A library made by MATLAB/GNU Octave which includes rigorous ODE solvers)

#数值微分方程 #常微分方程

评论 (0)

  • 还没有评论，来抢沙发吧。