常微分方程數值方法
常微分方程數值方法是用以尋找常微分方程(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)