數值穩定性
編輯
數值穩定性與數值誤差密切相關。使用有限差分方法進行計算時，若任意時間步的誤差不會導致其後計算結果的發散，則可稱該有限差分法是數值穩定的。如果誤差隨着進一步計算降低最終消失，該算法被認為穩定；若誤差在進計算中保持為常量，則認為該算法「中性穩定」。但如果誤差隨着進一步計算增長，結果發散，則數值方法不穩定。數值方法的穩定性可以通過馮諾依曼穩定性分析得到驗證。穩定性一般不易分析，特別是針對非線性偏微分方程。
馮諾依曼穩定性方法只適用於滿足 Lax–Richtmyer 條件 （Lax 等價定理 ） 的某些特殊差分法： 偏微分方程系統須線性，常係數，滿足周期性邊界條件，只有兩個獨立變量，差分法中最多含兩層時間步[4]
方法描述
編輯
馮諾依曼誤差分析將誤差分解為傅立葉級數 。為了描述此過程，考慮一維熱傳導方程
∂
u
∂
t
=
α
∂
2
u
∂
x
2
{\displaystyle {\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}}
空間網格間隔為
L
{\displaystyle L}
FTCS （Forward-Time Central-Space，時間步前向歐拉法，空間步三點中心差分) 離散處理，
(
1
)
u
j
n
+
1
=
u
j
n
+
r
(
u
j
+
1
n
−
2
u
j
n
+
u
j
−
1
n
)
{\displaystyle \quad (1)\qquad u_{j}^{n+1}=u_{j}^{n}+r\left(u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}\right)}
其中
r
=
α
Δ
t
Δ
x
2
{\displaystyle r={\frac {\alpha \,\Delta t}{\Delta x^{2}}}}
u
j
n
{\displaystyle u_{j}^{n}}
u
(
x
,
t
)
{\displaystyle u(x,t)}
定義捨入誤差
ϵ
j
n
=
N
j
n
−
u
j
n
{\displaystyle \epsilon _{j}^{n}=N_{j}^{n}-u_{j}^{n}}
u
j
n
{\displaystyle u_{j}^{n}}
N
j
n
{\displaystyle N_{j}^{n}}
u
j
n
{\displaystyle u_{j}^{n}}
ϵ
j
n
{\displaystyle \epsilon _{j}^{n}}
[5]
(
2
)
ϵ
j
n
+
1
=
ϵ
j
n
+
r
(
ϵ
j
+
1
n
−
2
ϵ
j
n
+
ϵ
j
−
1
n
)
{\displaystyle \quad (2)\qquad \epsilon _{j}^{n+1}=\epsilon _{j}^{n}+r\left(\epsilon _{j+1}^{n}-2\epsilon _{j}^{n}+\epsilon _{j-1}^{n}\right)}
此式將確定誤差的遞推關係。方程 (1) 和 (2) 中，誤差和數值解隨時間具有一致的變化趨勢。對於含周期性邊界條件的線性微分方程，間隔
L
{\displaystyle L}
(
3
)
ϵ
(
x
)
=
∑
m
=
1
M
A
m
e
i
k
m
x
{\displaystyle \quad (3)\qquad \epsilon (x)=\sum _{m=1}^{M}A_{m}e^{ik_{m}x}}
其中波數
k
m
=
π
m
L
{\displaystyle k_{m}={\frac {\pi m}{L}}}
m
=
1
,
2
,
…
,
M
{\displaystyle m=1,2,\ldots ,M}
M
=
L
/
Δ
x
{\displaystyle M=L/\Delta x}
A
m
{\displaystyle A_{m}}
(
4
)
ϵ
(
x
,
t
)
=
∑
m
=
1
M
e
a
t
e
i
k
m
x
{\displaystyle \quad (4)\qquad \epsilon (x,t)=\sum _{m=1}^{M}e^{at}e^{ik_{m}x}}
其中
a
{\displaystyle a}
由於誤差所滿足的差分方程是線性的（級數每一項的行為與整個級數一致），只估計一項的誤差變化便足以估計整體趨勢：
(
5
)
ϵ
m
(
x
,
t
)
=
e
a
t
e
i
k
m
x
.
{\displaystyle \quad (5)\qquad \epsilon _{m}(x,t)=e^{at}e^{ik_{m}x}.}
為找出誤差隨時間步的變化, 將方程 (5) 式應用於離散後的誤差表達式上
ϵ
j
n
=
e
a
t
e
i
k
m
x
ϵ
j
n
+
1
=
e
a
(
t
+
Δ
t
)
e
i
k
m
x
ϵ
j
+
1
n
=
e
a
t
e
i
k
m
(
x
+
Δ
x
)
ϵ
j
−
1
n
=
e
a
t
e
i
k
m
(
x
−
Δ
x
)
,
{\displaystyle {\begin{aligned}\epsilon _{j}^{n}&=e^{at}e^{ik_{m}x}\\\epsilon _{j}^{n+1}&=e^{a(t+\Delta t)}e^{ik_{m}x}\\\epsilon _{j+1}^{n}&=e^{at}e^{ik_{m}(x+\Delta x)}\\\epsilon _{j-1}^{n}&=e^{at}e^{ik_{m}(x-\Delta x)},\end{aligned}}}
再代入到 (2) 式中，求解方程後可得
(
6
)
e
a
Δ
t
=
1
+
α
Δ
t
Δ
x
2
(
e
i
k
m
Δ
x
+
e
−
i
k
m
Δ
x
−
2
)
.
{\displaystyle \quad (6)\qquad e^{a\Delta t}=1+{\frac {\alpha \Delta t}{\Delta x^{2}}}\left(e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}-2\right).}
使用已知的指數三角關係式
cos
(
k
m
Δ
x
)
=
e
i
k
m
Δ
x
+
e
−
i
k
m
Δ
x
2
{\displaystyle \qquad \cos(k_{m}\Delta x)={\frac {e^{ik_{m}\Delta x}+e^{-ik_{m}\Delta x}}{2}}}
sin
2
k
m
Δ
x
2
=
1
−
cos
(
k
m
Δ
x
)
2
{\displaystyle \sin ^{2}{\frac {k_{m}\Delta x}{2}}={\frac {1-\cos(k_{m}\Delta x)}{2}}}
可以將方程 (6) 變作
(
7
)
e
a
Δ
t
=
1
−
4
α
Δ
t
Δ
x
2
sin
2
(
k
m
Δ
x
/
2
)
.
{\displaystyle \quad (7)\qquad e^{a\Delta t}=1-{\frac {4\alpha \Delta t}{\Delta x^{2}}}\sin ^{2}(k_{m}\Delta x/2).}
定義漲幅因子
G
≡
ϵ
j
n
+
1
ϵ
j
n
,
{\displaystyle G\equiv {\frac {\epsilon _{j}^{n+1}}{\epsilon _{j}^{n}}},}
則誤差有限的充要條件為
|
G
|
≤
1
{\displaystyle \vert G\vert \leq 1}
(
8
)
G
=
e
a
(
t
+
Δ
t
)
e
i
k
m
x
e
a
t
e
i
k
m
x
=
e
a
Δ
t
,
{\displaystyle \quad (8)\qquad G={\frac {e^{a(t+\Delta t)}e^{ik_{m}x}}{e^{at}e^{ik_{m}x}}}=e^{a\Delta t},}
聯立 (7) 和 (8) 兩式，易得穩定性條件為
(
9
)
|
1
−
4
α
Δ
t
Δ
x
2
sin
2
(
k
m
Δ
x
/
2
)
|
≤
1
{\displaystyle \quad (9)\qquad \left\vert 1-{\frac {4\alpha \Delta t}{\Delta x^{2}}}\sin ^{2}(k_{m}\Delta x/2)\right\vert \leq 1}
即
(
10
)
α
Δ
t
Δ
x
2
≤
1
2
.
{\displaystyle \quad (10)\qquad {\frac {\alpha \Delta t}{\Delta x^{2}}}\leq {\frac {1}{2}}.}
(10) 即為該算法的穩定性條件。 對於 FTCS 求解一維熱傳導方程，給定
Δ
x
{\displaystyle \Delta x}
Δ
t
{\displaystyle \Delta t}
參考資料
編輯
^ Analysis of Numerical Methods by E. Isaacson, H. B. Keller . [2011-05-20 ] . （原始內容存檔 於2011-05-21）. ^
Crank, J.; Nicolson, P., A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of Heat Conduction Type, Proc. Camb. Phil. Soc., 1947, 43 : 50–67, doi :10.1007/BF02127704
^
Charney, J. G.; Fjørtoft, R.; von Neumann, J., Numerical Integration of the Barotropic Vorticity Equation , Tellus, 1950, 2 : 237–254
^
Smith, G. D., Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd ed.: 67–68, 1985
^ Anderson, J. D., Jr. Computational Fluid Dynamics: The Basics with Applications. McGraw Hill . 1994.
