常微分方程数值方法
常微分方程数值方法是用以寻找常微分方程(ODE)解的数值近似值的方法。其使用也称作“数值积分”,不过“数值积分”主要是指积分的计算。
很多微分方程无法精确求解。但在工程学等领域的实际应用中,通常只需得到数值近似解。本文介绍的算法可用于计算这种近似值,另一种方法是用微积分技术得到解的级数展开表达。
常微分方程出现于物理学、化学、生物学、经济学等学科中。[1]此外,偏微分方程数值方法中的一部分将偏微分方程转为常微分方程,然后可用本文所述方法求解。
问题
一阶微分方程是有如下形式的初值问题(IVP):[2]:533–655
其中f是函数,初始条件是已知向量。“一阶”是说方程中只出现了y的一阶导。
在不失高阶系统一般性的前提下,限制(l)式为一阶微分方程,因为高阶ODE可通过引入额外变量转换为更大的一阶方程组。例如,二阶方程可重写为2个一阶方程:。
本文介绍IVP的数值方法,并指出边值问题(BVP)需要一套不同的工具:需要在多个点上定义解y的值或成分,因此要用不同方法求解BVP,如打靶法(及其变体),或差分、[3]伽辽金法、[4]配置法之类全局方法,都适用于此类问题。
Picard–Lindelöf定理指出,只要f是利普希茨连续的,就有唯一解。
方法
解一阶IVP的数值方法可分为两大类:[5]线性多步法与龙格-库塔法。还可进一步划为显式或隐式,例如隐式线性多步法包括亚当斯-莫尔顿法(Adams-Moulton methods)与向后微分公式(BDF),隐式龙格-库塔法[6]:204–215包括对角隐式龙格-库塔法(DIRK)、[7][8]单对角隐式龙格-库塔法(SDIRK)[9]与基于高斯求积[10]的高斯–拉道法[11]等等。线性多步法中的显式方法有亚当斯-巴什福思法,布彻表(Butcher tableau)为下对角的龙格-库塔法都是显式方法。根据经验,刚性微分方程需要用隐式方案,而非刚性问题则可用显示方案更有效地求解。
欧拉法
从曲线上任意一点出发,沿与曲线相切的直线移动一小段距离,就能得到曲线上邻近点的近似值。
重新排列后得到以下公式
利用(1)可得
此式的应用通常如下:择步长h,构造序列记精确解的数值估计值为。受(3)启发,可用下面的递归方法计算估计值:
这就是(前向)欧拉法,是莱昂哈德·欧拉(1768)描述的方法。
欧拉法是显式方法的一个例子,这是说新值是根据已知值(如)定义的。
反向欧拉法
若不用(2),而用近似值
则得到反向欧拉法:
反向欧拉法是隐式方法,这是说需要求解一个方程才能得到新值。通常用定点迭代或牛顿-拉弗森法(的某种修改版)实现之。
隐式方法求解这方程比显示方法直接代入要花更多时间,选择方法时必须考虑这一成本。隐式方法(如(6))的优点是求解刚性方程时通常更稳定,便可以使用更大的步长h。
一阶指数积分法
指数积分描述了一大类积分器,近来得到了长足发展。[13]它们至少可以追溯到20世纪60年代。
此处不用(1),假设微分方程形式为
或已被局部线性化,围绕背景状态产生线性项与非线性项。
将(7)乘以,并在时间区间内精确积分,便构造得到了指数积分:
此积分方程是精确的,但并没有定义积分。
使在整个区间内不变,可实现一阶指数积分:
推广
欧拉法往往不够精确。更准确地说,只有一阶(下面将介绍阶的概念),这就促使数学家寻找高阶方法。
一种方法是,用以及更多之前的值确定,所谓多步法便实现了这种想法。最简单的可能是蛙跳积分法,其是二阶方法,(粗略地说)依赖于2个时间值。
几乎所有实用的多步法都属于线性多步法,形式为
另一种方法是在区间内取更多点。这产生了得名于卡尔·龙格与马丁·威廉·库塔的龙格-库塔法,其中一种4阶方法尤为流行。
高级特征
要用这些方法求解ODE,需要的不仅是时间步长公式。始终相同的步长效率不高,于是开发了可变步长方法。通常,步长的选择应使每步的(局部)误差低于某容差水平,这意味着方法还要计算误差指标(error indicator),即对局部误差的估计。
这一思想的延伸是在不同阶方法之间动态选择(称作可变阶数方法)。基于理查德森外推法[14]的Bulirsch–Stoer算法之类方法[15][16]通常用于构建各种不同阶的方法。
其他理想特征还有:
其他方法
很多方法不在讨论的框架内,如
- 多阶导数法,不仅使用f,还用其导数。这类方法包括Hermite–Obreschkoff法、费尔贝格法与递归计算解y的泰勒级数系数的Parker–Sochacki法[17]、Bychkov–Scherbakov法等。
- 二阶常微分方程组法。前面说过,所有高阶常微分方程都可化为形式如(1)的一阶常微分方程组,这固然没错,但未必是最佳方法。尼斯特伦法可直接用于二阶方程。
- 几何积分法[18][19]专门针对几类特殊的常微分方程(如解哈密顿力学方程的辛积分法),这些方法在数值求解时会尊重这些类的底结构或几何结构。
- 量化状态系统方法是基于状态量化思想的一系列ODE积分方法,在模拟频繁不连续的稀疏系统时非常有效。
实时并行方法
有些IVP要求积分具有很高的时间分辨率和/或很长的时间区间,传统的序列时间步进法无法实时计算(如数值天气预报、等离子体建模与分子动力学中的IVP)。针对这些问题,人们开发了实时并行(PinT)法以便用并行计算缩短运行时间。
早期的PinT法(最早提出于20世纪60年代)[20]最初被研究人员忽视,因为其所需的并行计算架构尚未普及。2000年代初,随着算力的提高,灵活易用的PinT算法——Parareal算法重新吸引了兴趣,它适用于求解各种IVP。百亿亿次级运算的出现使PinT算法获得更多关注,并启动了能用于世界上最强大的超级计算机的算法开发。截至2023年,最流行的方法有Parareal、PFASST、ParaDiag、MGRIT等,[21]
分析
数值分析不仅包括数值方法的设计,还包括其分析。分析中的3个核心概念是:
收敛性
若数值解随着步长h趋近于0而逼近精确解,则称此数值方法具有收敛性(convergent)。更确切地说,要求对利普希茨函数f与每个,ODE (1)
上述所有方法都收敛。
一致性与阶数
设数值方法
方法的局部(截断)误差定义为迭代一步产生的误差。即,假设之前的迭代无误差,则是此方法结果与精确解之间的差
若
则称此方法一致(consistent)。若
则称此方法阶数为p。因此,阶数不为0的方法是一致的。上文介绍的两种欧拉法(4、6)都是1阶的,因此是一致的。实践中使用的大多数方法阶数更高。一致性是收敛的必要条件[来源请求],但不是充分条件;方法要收敛,必须同时具有一致性与零稳性(zero-stable)。
一个相关概念是全局(截断)误差,即达到固定时间t所需所有步骤中持续存在的误差。明确地说,t时刻的全局误差是,其中。第p阶一步法是收敛的,其全局误差是。对多阶方法,这一说法不一定成立。
稳定性与刚性
对某些微分方程,欧拉法、显式龙格-库塔法、多步法(如亚当斯-巴什福思法)之类标准方法会表现出解的不稳定性,其他方法则可能产生稳定的解。方程中的这种“困难行为”(本身不一定复杂)称作“刚性”(stiffness),通常是由于底问题中存在不同时间尺度造成的。[23]例如,机械系统中的碰撞(如碰撞振子中的)发生的时间尺度通常比物体运动的时间尺度小得多,这种差异使状态参数曲线变得非常陡峭。
刚性问题在化学动力学、控制论、固体力学、天气预报、生物学、等离子体物理、电子学等领域中无处不在。克服刚性的一种方法是将微分方程推广到微分包含式,从而允许非光滑性并建立其模型。[24][25]
历史
- 1768 - 莱昂哈德·欧拉发现欧拉法。
- 1824 - 奥古斯丁-路易·柯西证明了欧拉法的收敛性。此证明中,柯西使用了隐性欧拉法。
- 1855 - Francis Bashforth在一封信中首次提到约翰·柯西·亚当斯的多步法。
- 1895 - 卡尔·龙格发现第一个龙格-库塔法。
- 1901 - 马丁·威廉·库塔描述了现在流行的4阶龙格-库塔法。
- 1910 - 路易斯·弗莱·理查德森公布了他的外推法——理查德森外推法。
- 1952 - Charles F. Curtiss与Joseph Oakland Hirschfelder提出刚性方程概念。
- 1963 - Germund Dahlquist引入积分法的A稳定性。
二阶一维边值问题的数值解法
边值问题(BVPs)通常要离散化,得到近似相等的矩阵问题再数值求解。[28]最常用的一维BVP数值求解方法称作有限差分法。[3]这种方法用点值的线性组合构造描述函数导数的有限差分系数,例如一阶导数的二阶中心差分近似为:
二阶导数的二阶中心差分近似为:
两式中,是离散域上相邻x值间的距离。这样就构建了线性系统,然后可用标准矩阵法求解。例如,待解方程:
下一步是将问题离散化,用线性导数近似,如
并解所得线性方程组。将有如下方程
初看之下,这个方程组似乎有困难,因为方程中没有不与变量相乘的项,但事实上这是错误的。i = 1或n − 1时,有一项涉及边值、,由于两值已知,可以简单地代入,就得到了具有非平凡解的非齐次线性方程组。
相关条目
注释
- ^ Chicone, C. (2006). Ordinary differential equations with applications (Vol. 34). Springer Science & Business Media.
- ^ Bradie (2006)
- ^ 3.0 3.1 LeVeque, R. J. (2007). Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems (Vol. 98). SIAM.
- ^ Slimane Adjerid and Mahboub Baccouch (2010) Galerkin methods. Scholarpedia, 5(10):10056.
- ^ Griffiths, D. F., & Higham, D. J. (2010). Numerical methods for ordinary differential equations: initial value problems. Springer Science & Business Media.
- ^ Hairer, Nørsett & Wanner (1993)
- ^ Alexander, R. (1977). Diagonally implicit Runge–Kutta methods for stiff ODE’s. SIAM Journal on Numerical Analysis, 14(6), 1006-1021.
- ^ Cash, J. R. (1979). Diagonally implicit Runge-Kutta formulae with error estimates. IMA Journal of Applied Mathematics, 24(3), 293-301.
- ^ Ferracina, L., & Spijker, M. N. (2008). Strong stability of singly-diagonally-implicit Runge–Kutta methods. Applied Numerical Mathematics, 58(11), 1675-1686.
- ^ Weisstein, Eric W. "Gaussian Quadrature." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/GaussianQuadrature.html (页面存档备份,存于互联网档案馆)
- ^ Everhart, E. (1985). An efficient integrator that uses Gauss-Radau spacings. In International Astronomical Union Colloquium (Vol. 83, pp. 185-202). Cambridge University Press.
- ^ Butcher, J. C. (1987). The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods. Wiley-Interscience.
- ^ Hochbruck (2010,第209–286页) This is a modern and extensive review paper for exponential integrators
- ^ Brezinski, C., & Zaglia, M. R. (2013). Extrapolation methods: theory and practice. Elsevier.
- ^ Monroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.
- ^ Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California.
- ^ Nurminskii, E. A., & Buryi, A. A. (2011). Parker-Sochacki method for solving systems of ordinary differential equations using graphics processors. Numerical Analysis and Applications, 4(3), 223.
- ^ Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations (Vol. 31). Springer Science & Business Media.
- ^ Hairer, E., Lubich, C., & Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numerica, 12, 399-450.
- ^ Nievergelt, Jürg. Parallel methods for integrating ordinary differential equations. Communications of the ACM. 1964, 7 (12): 731–733. S2CID 6361754. doi:10.1145/355588.365137 .
- ^ Parallel-in-Time.org. Parallel-in-Time.org. [2023-11-15]. (原始内容存档于2023-11-15).
- ^ Higham, N. J. (2002). Accuracy and stability of numerical algorithms (Vol. 80). SIAM.
- ^ Miranker, A. (2001). Numerical Methods for Stiff Equations and Singular Perturbation Problems: and singular perturbation problems (Vol. 5). Springer Science & Business Media.
- ^ Markus Kunze; Tassilo Kupper. Non-smooth Dynamical Systems: An Overview. Bernold Fiedler (编). Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems. Springer Science & Business Media. 2001: 431. ISBN 978-3-540-41290-8.
- ^ Thao Dang. Model-Based Testing of Hybrid Systems. Justyna Zander, Ina Schieferdecker and Pieter J. Mosterman (编). Model-Based Testing for Embedded Systems. CRC Press. 2011: 411. ISBN 978-1-4398-1845-9.
- ^ Brezinski, C., & Wuytack, L. (2012). Numerical analysis: Historical developments in the 20th century. Elsevier.
- ^ Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.
- ^ Ascher, U. M., Mattheij, R. M., & Russell, R. D. (1995). Numerical solution of boundary value problems for ordinary differential equations. Society for Industrial and Applied Mathematics.
参考文献
- Bradie, Brian. A Friendly Introduction to Numerical Analysis. Upper Saddle River, New Jersey: Pearson Prentice Hall. 2006. ISBN 978-0-13-013054-9.
- J. C. Butcher, Numerical methods for ordinary differential equations, ISBN 0-471-96758-0
- Ernst Hairer, Syvert Paul Nørsett and Gerhard Wanner, Solving ordinary differential equations I: Nonstiff problems, second edition, Springer Verlag, Berlin, 1993. ISBN 3-540-56670-8.
- Ernst Hairer and Gerhard Wanner, Solving ordinary differential equations II: Stiff and differential-algebraic problems, second edition, Springer Verlag, Berlin, 1996. ISBN 3-540-60452-9.
(This two-volume monograph systematically covers all aspects of the field.) - Hochbruck, Marlis; Ostermann, Alexander. Exponential integrators. Acta Numerica. May 2010, 19: 209–286. Bibcode:2010AcNum..19..209H. CiteSeerX 10.1.1.187.6794 . S2CID 4841957. doi:10.1017/S0962492910000048.
- Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996. ISBN 0-521-55376-8 (hardback), ISBN 0-521-55655-4 (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. ISBN 0-471-92990-5.
(Textbook, slightly more demanding than the book by Iserles.)
外部链接
- Joseph W. Rudmin, Application of the Parker–Sochacki Method to Celestial Mechanics Portuguese Web Archive的存档,存档日期2016-05-16, 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.)
- Pchelintsev, A.N. An accurate numerical method and algorithm for constructing solutions of chaotic systems. Journal of Applied Nonlinear Dynamics. 2020, 9 (2): 207–221. S2CID 225853788. arXiv:2011.10664 . doi:10.5890/JAND.2020.06.004.
- GitHub上的kv页面 (C++ library with rigorous ODE solvers)
- INTLAB (页面存档备份,存于互联网档案馆) (A library made by MATLAB/GNU Octave which includes rigorous ODE solvers)