常微分方程數值方法

求常微分方程数值解的方法

常微分方程數值方法是用以尋找常微分方程(ODE)解的數值近似值的方法。其使用也稱作「數值積分」,不過「數值積分」主要是指積分的計算。

微分方程的數值積分求解圖。
  藍:歐拉法
  綠:中點法
  紅:精確解:
步長
時的同一圖。時,中點法比歐拉法收斂得更快。

很多微分方程無法精確求解。但在工程學等領域的實際應用中,通常只需得到數值近似解。本文介紹的算法可用於計算這種近似值,另一種方法是用微積分技術得到解的級數展開表達。

常微分方程出現於物理學化學生物學經濟學等學科中。[1]此外,偏微分方程數值方法中的一部分將偏微分方程轉為常微分方程,然後可用本文所述方法求解。

問題

編輯

一階微分方程是有如下形式的初值問題(IVP):[2]:533–655

  1

其中f是函數 ,初始條件 是已知向量。「一階」是說方程中只出現了y的一階導。

在不失高階系統一般性的前提下,限制(l)式為一階微分方程,因為高階ODE可通過引入額外變量轉換為更大的一階方程組。例如,二階方程 可重寫為2個一階方程: 

本文介紹IVP的數值方法,並指出邊值問題(BVP)需要一套不同的工具:需要在多個點上定義解y的值或成分,因此要用不同方法求解BVP,如打靶法(及其變體),或差分[3]伽遼金法[4]配置法英語Collocation method之類全局方法,都適用於此類問題。

Picard–Lindelöf定理指出,只要f利普希茨連續的,就有唯一解。

方法

編輯

解一階IVP的數值方法可分為兩大類:[5]線性多步法英語Linear multistep method龍格-庫塔法。還可進一步劃為顯式或隱式,例如隱式線性多步法包括亞當斯-莫爾頓法(Adams-Moulton methods)與向後微分公式英語Backward differentiation formula(BDF),隱式龍格-庫塔法[6]:204–215包括對角隱式龍格-庫塔法(DIRK)、[7][8]單對角隱式龍格-庫塔法(SDIRK)[9]與基於高斯求積[10]的高斯–拉道法[11]等等。線性多步法中的顯式方法有亞當斯-巴什福思法,布徹表(Butcher tableau)為下對角的龍格-庫塔法都是顯式方法。根據經驗,剛性微分方程需要用隱式方案,而非剛性問題則可用顯示方案更有效地求解。

所謂一般線性方法英語General linear methods(GLM)是上述兩大類方法的概括。[12]

歐拉法

編輯

從曲線上任意一點出發,沿與曲線相切的直線移動一小段距離,就能得到曲線上鄰近點的近似值。

從微分方程(1)開始,可用差分近似代替導數y′:

  2

重新排列後得到以下公式

 

利用(1)可得

  3

此式的應用通常如下:擇步長h,構造序列 記精確解 的數值估計值為 。受(3)啟發,可用下面的遞歸方法計算估計值:

  4

這就是(前向)歐拉法,是萊昂哈德·歐拉(1768)描述的方法。

歐拉法是顯式方法的一個例子,這是說新值 是根據已知值(如 )定義的。

反向歐拉法

編輯

若不用(2),而用近似值

  5

則得到反向歐拉法

  6

反向歐拉法是隱式方法,這是說需要求解一個方程才能得到新值 。通常用定點迭代牛頓-拉弗森法(的某種修改版)實現之。

隱式方法求解這方程比顯示方法直接代入要花更多時間,選擇方法時必須考慮這一成本。隱式方法(如(6))的優點是求解剛性方程時通常更穩定,便可以使用更大的步長h

一階指數積分法

編輯

指數積分英語Exponential integrator描述了一大類積分器,近來得到了長足發展。[13]它們至少可以追溯到20世紀60年代。

此處不用(1),假設微分方程形式為

  7

或已被局部線性化,圍繞背景狀態產生線性項 與非線性項 

將(7)乘以 ,並在時間區間 內精確積分,便構造得到了指數積分:

 

此積分方程是精確的,但並沒有定義積分。

使 在整個區間內不變,可實現一階指數積分:

  8

推廣

編輯

歐拉法往往不夠精確。更準確地說,只有一階(下面將介紹階的概念),這就促使數學家尋找高階方法。

一種方法是,用 以及更多之前的值確定 ,所謂多步法便實現了這種想法。最簡單的可能是蛙跳積分法,其是二階方法,(粗略地說)依賴於2個時間值。

幾乎所有實用的多步法都屬於線性多步法英語Linear multistep method,形式為  

另一種方法是在區間 內取更多點。這產生了得名於卡爾·龍格馬丁·威廉·庫塔龍格-庫塔法,其中一種4階方法尤為流行。

高級特徵

編輯

要用這些方法求解ODE,需要的不僅是時間步長公式。始終相同的步長效率不高,於是開發了可變步長方法。通常,步長的選擇應使每步的(局部)誤差低於某容差水平,這意味着方法還要計算誤差指標(error indicator),即對局部誤差的估計。

這一思想的延伸是在不同階方法之間動態選擇(稱作可變階數方法)。基於理查德森外推法[14]Bulirsch–Stoer算法之類方法[15][16]通常用於構建各種不同階的方法。

其他理想特徵還有:

  • 輸出稠密:不僅對點 ,還對整個積分區間進行低成本數值逼近;
  • 事件定位:找到事件發生的位置,例如某函數取0的時間。通常用求根算法
  • 支持並行計算
  • 用於時間積分時,具有時間可逆性。

其他方法

編輯

很多方法不在討論的框架內,如

實時並行方法

編輯

有些IVP要求積分具有很高的時間分辨率和/或很長的時間區間,傳統的序列時間步進法無法實時計算(如數值天氣預報、等離子體建模與分子動力學中的IVP)。針對這些問題,人們開發了實時並行(PinT)法以便用並行計算縮短運行時間。

早期的PinT法(最早提出於20世紀60年代)[20]最初被研究人員忽視,因為其所需的並行計算架構尚未普及。2000年代初,隨着算力的提高,靈活易用的PinT算法——Parareal算法英語Parareal重新吸引了興趣,它適用於求解各種IVP。百億億次級運算英語Exascale computing的出現使PinT算法獲得更多關注,並啟動了能用於世界上最強大的超級計算機的算法開發。截至2023年,最流行的方法有Parareal、PFASST、ParaDiag、MGRIT等,[21]

分析

編輯

數值分析不僅包括數值方法的設計,還包括其分析。分析中的3個核心概念是:

  • 收斂性:方法是否逼近解;
  • 階數:近似解的程度;
  • 數值穩定性:誤差是否能得到抑制。[22]

收斂性

編輯

若數值解隨着步長h趨近於0而逼近精確解,則稱此數值方法具有收斂性(convergent)。更確切地說,要求對利普希茨函數f與每個 ,ODE (1)

 

上述所有方法都收斂。

一致性與階數

編輯

設數值方法

 

方法的局部(截斷)誤差定義為迭代一步產生的誤差。即,假設之前的迭代無誤差,則是此方法結果與精確解之間的差

 

 

則稱此方法一致(consistent)。若

 

則稱此方法階數為p。因此,階數不為0的方法是一致的。上文介紹的兩種歐拉法(4、6)都是1階的,因此是一致的。實踐中使用的大多數方法階數更高。一致性是收斂的必要條件[來源請求],但不是充分條件;方法要收斂,必須同時具有一致性與零穩性(zero-stable)。

一個相關概念是全局(截斷)誤差,即達到固定時間t所需所有步驟中持續存在的誤差。明確地說,t時刻的全局誤差是 ,其中 。第p階一步法是收斂的,其全局誤差是 。對多階方法,這一說法不一定成立。

穩定性與剛性

編輯

對某些微分方程,歐拉法、顯式龍格-庫塔法、多步法(如亞當斯-巴什福思法)之類標準方法會表現出解的不穩定性,其他方法則可能產生穩定的解。方程中的這種「困難行為」(本身不一定複雜)稱作「剛性」(stiffness),通常是由於底問題中存在不同時間尺度造成的。[23]例如,機械系統中的碰撞(如碰撞振子中的)發生的時間尺度通常比物體運動的時間尺度小得多,這種差異使狀態參數曲線變得非常陡峭。

剛性問題在化學動力學控制論固體力學天氣預報生物學等離子體物理、電子學等領域中無處不在。克服剛性的一種方法是將微分方程推廣到微分包含式,從而允許非光滑性並建立其模型。[24][25]

歷史

編輯

下面是該領域一些重要進展的時間線。[26][27]

二階一維邊值問題的數值解法

編輯

邊值問題(BVPs)通常要離散化,得到近似相等的矩陣問題再數值求解。[28]最常用的一維BVP數值求解方法稱作有限差分法[3]這種方法用點值的線性組合構造描述函數導數的有限差分係數,例如一階導數的二階中心差分近似為:

 

二階導數的二階中心差分近似為:

 

兩式中, 是離散域上相鄰x值間的距離。這樣就構建了線性系統,然後可用標準矩陣法求解。例如,待解方程:

 

下一步是將問題離散化,用線性導數近似,如

 

並解所得線性方程組。將有如下方程

 

初看之下,這個方程組似乎有困難,因為方程中沒有不與變量相乘的項,但事實上這是錯誤的。i = 1或n − 1時,有一項涉及邊值  ,由於兩值已知,可以簡單地代入,就得到了具有非平凡解的非齊次線性方程組

相關條目

編輯

注釋

編輯
  1. ^ Chicone, C. (2006). Ordinary differential equations with applications (Vol. 34). Springer Science & Business Media.
  2. ^ Bradie (2006)
  3. ^ 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.
  4. ^ Slimane Adjerid and Mahboub Baccouch (2010) Galerkin methods. Scholarpedia, 5(10):10056.
  5. ^ Griffiths, D. F., & Higham, D. J. (2010). Numerical methods for ordinary differential equations: initial value problems. Springer Science & Business Media.
  6. ^ Hairer, Nørsett & Wanner (1993)
  7. ^ Alexander, R. (1977). Diagonally implicit Runge–Kutta methods for stiff ODE’s. SIAM Journal on Numerical Analysis, 14(6), 1006-1021.
  8. ^ Cash, J. R. (1979). Diagonally implicit Runge-Kutta formulae with error estimates. IMA Journal of Applied Mathematics, 24(3), 293-301.
  9. ^ Ferracina, L., & Spijker, M. N. (2008). Strong stability of singly-diagonally-implicit Runge–Kutta methods. Applied Numerical Mathematics, 58(11), 1675-1686.
  10. ^ Weisstein, Eric W. "Gaussian Quadrature." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/GaussianQuadrature.html頁面存檔備份,存於網際網路檔案館
  11. ^ Everhart, E. (1985). An efficient integrator that uses Gauss-Radau spacings. In International Astronomical Union Colloquium (Vol. 83, pp. 185-202). Cambridge University Press.
  12. ^ Butcher, J. C. (1987). The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods. Wiley-Interscience.
  13. ^ Hochbruck (2010,第209–286頁) This is a modern and extensive review paper for exponential integrators
  14. ^ Brezinski, C., & Zaglia, M. R. (2013). Extrapolation methods: theory and practice. Elsevier.
  15. ^ Monroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.
  16. ^ Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California.
  17. ^ 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.
  18. ^ Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations (Vol. 31). Springer Science & Business Media.
  19. ^ Hairer, E., Lubich, C., & Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numerica, 12, 399-450.
  20. ^ 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 . 
  21. ^ Parallel-in-Time.org. Parallel-in-Time.org. [2023-11-15]. (原始內容存檔於2023-11-15). 
  22. ^ Higham, N. J. (2002). Accuracy and stability of numerical algorithms (Vol. 80). SIAM.
  23. ^ Miranker, A. (2001). Numerical Methods for Stiff Equations and Singular Perturbation Problems: and singular perturbation problems (Vol. 5). Springer Science & Business Media.
  24. ^ 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. 
  25. ^ 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. 
  26. ^ Brezinski, C., & Wuytack, L. (2012). Numerical analysis: Historical developments in the 20th century. Elsevier.
  27. ^ Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.
  28. ^ 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.)

外部連結

編輯