常微分方程數值方法

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

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

外部連結

編輯