顯式方法(explicit method)和隱式方法(implicit methods)是数值分析中計算以時間為自變數的常微分方程和偏微分方程的數值近似法,也是偏微分方程中计算机模拟會使用的方法。顯式方法會用系統目前的狀態來計算下一個時間的狀態,隱式方法會將系統目前狀態和下一個時間的狀態以方程式的方式表示,下一個時間的狀態為未知數,求解方程式來得到下一個時間的狀態。考慮數學的型式,若Y(t)是目前系統狀態,Y(t+\Delta t)是下一個時間的狀態(\Delta t是很小的時間間隔),在顯式方法下,下一個時間的狀態為
: Y(t+\Delta t) = F(Y(t))\,
在隱式方法下,下一個時間的狀態可用以下方程來表示
: G\Big(Y(t), Y(t+\Delta t)\Big)=0 \qquad (1)\,
需求解方程才能得到下一個時間的狀態。
隱式方法需要針對方程式求解,需要額外的計算,也比較不容易實現。顯式方法比較容易實現,但許多問題屬於刚性方程,為了要使其誤差限制在一定範圍內(和数值稳定性有關),需要非常小的\Delta t。這類的問題,若要得到相同的精度,可以用隱式方法,選取較大的時間間隔,就算將隱式方法需要在每一步針對方程式(1)的求解考慮在內,隱式方法仍可以在較少的運算量下得到結果。因此,要用隱式方法或是顯式方法求解需視問題而定。隱式方法無法用在每一種微分運算子上,因此有時會用所謂的運算子分離法(operator splitting method),將微分算子改寫為兩種互補算子的結合
:Y(t+\Delta t) = F(Y(t+\Delta t))+G(Y(t)),\,
一個是顯式的,另一個則是隱式的。
一般的應用會讓隱式項是線性的,而顯式項可以為非線性。這種組合稱為「隱式-顯式方法」(Implicit-Explicit Method),簡稱IMEX。
用前向歐拉方法和後向歐拉方法的說明
考慮常微分方程
: \frac{dy}{dt} = -y^2, \ t\in [0, a]\quad \quad (2)
初始條件是y(0)=1。考慮格點t_k=a\frac{k}{n},0 ≤ k ≤ n,意思是,時間間隔是\Delta t=a/n,,且y_k=y(t_k)。用最簡單的顯式和隱式方法將此方程式离散化,分別是「前向歐拉方法」及「後向歐拉方法」,並且比較其差異。
;前向歐拉方法:
[[File:Result of applying integration schemes.png|thumb|用不同的積分法所得的結果 y'=-y^2, \; t\in[0, 5], \; y_0=1 ,\Delta t = 5/10.]]
前向欧拉方法
:\left(\frac{dy}{dt}\right)_k \approx \frac{y_{k+1}-y_k}{\Delta t} = - y_k^2
可得
: y_{k+1}=y_k-\Delta t y_k^2 \quad \quad \quad(3)\,
對所有k=0, 1, \dots, n.,這是y_{k+1}的顯式公式。
;後向歐拉方法:
用後向歐拉方法
:\frac{y_{k+1}-y_k}{\Delta t} = - y_{k+1}^2
可以得到y_{k+1}的隱式方程
: y_{k+1}+\Delta t y_{k+1}^2=y_k
比較上式和公式(3),公式(3)的y_{k+1}可以直接求得,而此處是方程式中的未知數,需要求解。
這是一元二次方程,有一個正根和一個負根,因為其初值為正,選擇其正根,則下一步的y為
: y_{k+1}=\frac{-1+\sqrt{1+4\Delta t y_k}}{2 \Delta t}. \quad \quad (4)
大部份的隱式方程中,要求解的方程會比一元二次方程複雜的多,也有可能不存在解析解,因此需要用其他求根算法(例如牛顿法)來求得數值解。
;克兰克-尼科尔森方法:
用克兰克-尼科尔森方法
:\frac{y_{k+1}-y_k}{\Delta t} = -\frac{1}{2}y_{k+1}^2 -\frac{1}{2}y_{k}^2
可以求得y_{k+1}的隱式方程
: y_{k+1}+\frac{1}{2}\Delta t y_{k+1}^2=y_k - \frac{1}{2}\Delta t y_{k}^2
這可以用求根算法(例如牛顿法)來求得y_{k+1}的數值解。
克兰克-尼科尔森方法可以視為是通用的IMEX(Implicit-Explicit,隱式-顯式)架構。
;前向-後向歐拉方法:
為了應用IMEX架構,考慮另一個微分方程:
: \frac{dy}{dt} = y-y^2, \ t\in [0, a]\quad \quad (5)
可以得到
: \left(\frac{dy}{dt}\right)_k \approx y_{k+1}-y_{k}^2, \ t\in [0, a]
因此
: y_{k+1}=\frac{y_k(1-y_k\Delta t)}{1-\Delta t} \quad\quad(6)
針對k=0, 1, \dots, n
相關條目
*
- SIMPLE算法,壓力耦合方程組的半隱式方法
來源
评论 (0)