Numerical Simulation
數值模擬。函數方程式，求得數值解。
Numerical Iteration：數值遞推。建立遞迴公式，逐步遞推求得數值解。
Numerical Error：數值誤差。數值遞推過程，數值解與符號解的差距。
Numerical Convergence：數值收斂。數值遞推過程，數值解與符號解相符。
Numerical Scheme：數值方案。精心設計的遞迴公式，使得數值收斂。
Numerical Iteration
數值遞推。製作遞迴公式，求得微分方程式的數值解。
據我所知，數值模擬目前只有這一套手法，沒有其他手法了。這一套手法基於Euler Method，建立遞迴公式，依序填值。
方才處理單變數函數。現在處理多變數函數，衍生大量細節。
為了做成動畫，時間變化放在左式，空間變化放在右式。
為了做成動畫，左式採用backward difference。
Recurrence
離散化。三種版本差分。
計算過程，建立表格，依序填值，得到數值解。
symbolic solution | numeric solution
------------------------------------| ----------------
f(t, x) | f[⋅][⋅]
f(t₀, x₀) | f[0][0]
f(tₙ, xᵢ) = f(t₀ + n Δt, x₀ + i Δx) | f[n][i]
符號解、數值解，數學符號改寫成遞推風格。
symbolic solution | numeric solution
------------------| -----------------
f⁽˙⁾(⋅) | f̂⁽˙⁾(⋅)
f⁽⁰⁾(x₀) | f̂⁽⁰⁾(x₀)
f⁽ⁿ⁾(xᵢ) | f̂⁽ⁿ⁾(xᵢ)
Boundary
邊界。計算範圍。輪廓形狀。
我們通常只關心一小塊範圍，而不是從負無限大到正無限大。
邊界形狀視實際問題而定。例如矩形、圓形、特殊造型。
其中最棒的形狀是一維區間、二維矩形、三維立方體，恰好符合程式語言的陣列格式。
f(t) 5 ≤ t ≤ 10 time interval
zero-dimensional space
= interval boundary
f(t,x) 5 ≤ t ≤ 10 time interval
-10 ≤ x ≤ +10 space interval
= rectangle boundary
f(t,x,y) 5 ≤ t ≤ 10 time interval
x² + y² ≤ 10 circular space
= cylinder boundary
Initial Value / Boundary Value
數值模擬當中，變數被分成時間變數與空間變數，邊界被分成初始值與邊界值。
Initial Value：初始值。時刻是t₀的時間變數邊界。
Boundary Value：邊界值。空間變數邊界。
數值模擬當中，問題被分成初始值問題與邊界值問題，邊界值被分成不是指定數值與是指定數值。前者較簡單，演算法請見上個章節。後者較困難，演算法更加複雜，此處省略。
Initial Value Problem：初始值問題。邊界值不是指定數值。取決於鄰居數值。
Boundary Value Problem：邊界值問題。邊界值是指定數值。
初始值視實際問題而定。例如零函數、脈衝函數、隨機函數。
(1) Zero Initial Condition 零初始條件（一律是零）
(2) Delta-pulse Initial Condition 脈衝初始條件（一處非零）
(3) Random Initial Condition 隨機初始條件（亂數）
其中最難纏的初始值是隨機函數。必須謹慎內插適當數值，以確保邊界內部數值連續。一種替代方案是【尚待確認】：函數值是隨機弦波疊加。
邊界值（不是指定數值）視實際問題而定。例如無界函數、對稱函數、週期函數。
(1) Open Boundary Condition 開放邊界條件：無限空間
(2) Symmetric Boundary Condition 對稱邊界條件：牆面反彈
(3) Periodic Boundary Condition 週期邊界條件：無限循環
其中最難纏的邊界值是無界函數。必須謹慎外插適當數值，以確保邊界內緣數值不受影響。一種替代方案是Sponge Layer：海綿層置於邊界外緣，海綿層的函數值快速衰減至零。
邊界值（指定數值）視實際問題而定。例如函數值、函數導數值。
(1) Dirichlet Boundary Condition 函數值
(2) Neumann Boundary Condition 函數導數值（的邊界法線分量）
(3) Cauchy Boundary Condition 上述兩者
(4) Robin Boundary Condition 上述兩者的加權總和
這些邊界值都一樣難纏。一種替代方案是Shooting Method：亂槍打鳥、試誤法。邊界值問題簡化為初始值問題。嘗試各種初始值（發射地點），各自進行數值遞推（行進軌道），檢查邊界值是否恰好符合指定數值（恰好射中邊界值）。
實際問題還有各式各樣的初始條件與邊界條件，但是目前尚未形成一套知識體系，無法為大家介紹。例如零散時間零散地點已知數值、外力即時影響數值。等你發表論文。
Numerical Error
數值誤差。符號解與數值解的差距。分為兩個階段。
local error: distance of 1st iteration of
a symbolic solution and a numeric solution.
‖f⁽¹⁾ - f̂⁽¹⁾‖ ≡ distance(f⁽¹⁾(x), f̂⁽¹⁾(x))
≡ distance(f(t₁, ⋅), f[1][⋅])
global error: distance of nth iteration of
a symbolic solution and a numeric solution.
‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ ≡ distance(f⁽ⁿ⁾(x), f̂⁽ⁿ⁾(x))
≡ distance(f(tₙ, ⋅), f[n][⋅])
局部誤差：遞推一步，符號解、數值解兩者差距。
全域誤差：遞推n步，符號解、數值解兩者差距。
順帶一提，累加n處local error並不等於global error。
Error Metric
局部誤差、全域誤差的距離函數，通常設定成平方距離。
沒有空間變數：兩個數字的距離。通常設定成absolute-value norm。
擁有空間變數：兩個函數的距離。通常設定成Hibert norm。
‖f⁽ⁿ⁾ - g⁽ⁿ⁾‖₂ = | f(tₙ) - g(tₙ) |
‖f⁽ⁿ⁾ - g⁽ⁿ⁾‖₂ = √ ∫ (f(tₙ,x) - g(tₙ,x))² dx
‖f⁽ⁿ⁾ - g⁽ⁿ⁾‖₂ = √ ∫∫ (f(tₙ,x,y) - g(tₙ,x,y))² dxdy
兩個離散化函數的距離：如法炮製。黎曼積分改成黎曼和。
我們習慣在邊界取樣，導致黎曼和在左右邊界各多出½Δx寬度。無傷大雅。如果你很在意，你可以扣掉它。
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖₂ = √ sum (f̂⁽ⁿ⁾(xᵢ) - ĝ⁽ⁿ⁾(xᵢ))² Δx
= √ sum (f[i] - g[i])² Δx
一個函數與一個離散化函數的距離：大家採用兩個離散化函數的距離，事情比較簡單。
‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖₂ = √ sum (f⁽ⁿ⁾(xᵢ) - f̂⁽ⁿ⁾(xᵢ))² Δx
Numerical Convergence
數值收斂。數值遞推過程，全域誤差趨近零，答案精準。
中文譯作「收斂」，原意應是「合流」。
Fundamental Theorem of Numerical Analysis
「數值分析基本定理」。
if consistent, then convergent = stable
拆成兩句話。本文主角是第二句話。
(1) consistent & convergent => stable
(2) consistent & stable => convergent
這三個詞彙的中文翻譯：
adjective 形容詞丨 noun 名詞
------------一一一丨-------------一一一
consistent 一致的丨 consistency 一致性
stable 穩定的丨 stability 穩定性
convergent 收斂的丨 convergence 收斂性
這三個詞彙的原始意義：
consistency: 1. theorems have no contradiction.
2. an equation have a unique solution.
3. symbol and numeral are approximately equal.
stability: two sequence are bounded by a certain distance
that is depending on the initial distance
|aₙ - bₙ| < ε when |a₀ - b₀| < δ
where ε and δ are constants
convergence: a sequence leads to a constant.
aₙ → C when n → ∞
where C is a constant
一致性：定理無矛盾。等式有唯一解。符號與數值相符。表裡如一。
穩定性：兩串數列，對應項差距範圍受限。牽繩遛狗。
收斂性：一串數列，數字趨近定值。閒逛到家。
這三個詞彙在數值模擬當中的用法：
（我沒有照抄文獻裡的句子。我不確定我的理解是否正確。）
differential equation: d/dt f(t,x) = F(f(t,x))
symbolic solution: f(t,x)
consistency: difference of the functional equation and
its discretization approaches zero.
L(f(t,x)) - L̂(f(t,x)) → 0 when Δt→0, Δx→0
where L(f(t,x)) = (d/dt f(t,x)) - F(f(t,x))
big-O notation is used somehow.
L(f(t,x)) = L̂(f(t,x)) + O((Δt)ᵖ + (Δx)𐞥)
where p>0, q>0
stability: distance of corresponding iteration of
two numeric solutions with
different initial conditions is bounded.
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ K ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ when n≥0
where K is a constant and 0 < K < 1
with consistency,
distance of corresponding iteration of
a symbolic solution and a numeric solution with
same initial condition is bounded.
‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ ≤ K ‖f⁽⁰⁾ - f̂⁽⁰⁾‖ when n≥0
where K is a constant and 0 < K < 1
convergence: distance of corresponding iteration of
a symbolic solution and a numeric solution with
same initial condition leads to zero.
‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ → 0 when n→∞
一致性：函數方程式、其離散化（遞迴公式），兩者趨近相等。
穩定性：兩組數值解（初始條件不同），各回合距離受限。
一致性且穩定性：符號解、數值解，各回合距離受限。
收斂性：符號解、數值解，各回合趨近相等。
一致性：截斷誤差趨近零。
穩定性：任兩數值解，各回合距離受限。
一致性且穩定性：全域誤差受限。
收斂性：全域誤差趨近零。
大部分文獻又進一步修改定義，將f(t,x)改成d/dt f(t,x)再改成微分方程式的右式F(f(t,x))。我推測原因是避免在Δt→0的情況下乘以除以Δt，並且確保時間變化空間變化可以分離成左式右式。
differential equation: d/dt f(t,x) = F(f(t,x))
consistency: F(f(t,x)) - F̂(f(t,x)) → 0 when Δx→0
stability: ‖F(f̂⁽ⁿ⁾) - F(ĝ⁽ⁿ⁾)‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ when n≥0
with consistency,
‖F(f⁽ⁿ⁾) - F(f̂⁽ⁿ⁾)‖ ≤ ‖f⁽⁰⁾ - f̂⁽⁰⁾‖ when n≥0
convergence: ‖F(f⁽ⁿ⁾) - F(f̂⁽ⁿ⁾)‖ → 0 when n→∞
將f(t,x)改成F(f(t,x))，此時穩定性宛如Lipschitz function。數值分析基本定理宛如Banach Fixed-point Theorem。「平緩函數套用不動點遞推法保證收斂」。
順帶一提，不動點定理則是源自monotone convergence theorem。「if monotone, then convergent = bounded」。
此處版本，遞迴公式不受時間影響，每回合一律使用相同泛函數F̂(f(t,x))。遞迴公式受時間影響，每回合都要微調泛函數F̂ₜ(f(t,x))，事情更加棘手，我沒有研究。穩定性也推廣成AN-stability和BN-stability。
證明我還沒有學會。等我學完泛函分析不知何年何月。Banach space、uniform boundedness principle、dominated convergence theorem。
總之，當一致性和穩定性同時成立，就能保證數值收斂。
接下來分別介紹一致性和穩定性。
Consistency
「一致性」。符號轉數值，符號與數值相符。分為兩點討論。
一、數字大小相符。
二、運算結果相符。加減乘除微積求值代入。
Consistency Error（Numerical Error）
「一致性誤差」。符號轉數值導致的誤差。分為兩種類型。
round-off error：四捨五入誤差。實數四則運算從符號轉數值導致的誤差。
truncation error：截斷誤差。函數微分運算從符號轉數值導致的誤差。
實數從符號轉數值，一般採用浮點數。四捨五入誤差，一般就是浮點數運算誤差。請見本站文件「Bit」。
函數微分運算從符號轉數值，一般採用泰勒近似。截斷誤差，一般就是泰勒近似誤差。請見本站文件「Approximation」。
Floating-point Arithmetic
「浮點數運算」。離散化實數之四則運算。
實數資料結構採用浮點數，實施四則運算。
本節不是本文重點。細節請自行上網搜尋講義。
Round-off Error
「四捨五入誤差」。浮點數運算結果將消滅多餘的低位數。
簡易統計方式：單次運算誤差ｘ運算次數＝總誤差。
機器精度ϵ，時間複雜度O(N)，總誤差O(ϵN)。
本節不是本文重點。細節請自行上網搜尋講義。
Finite Difference
「有限差分」。離散化函數之微分運算。
函數資料結構採用Regular Discretization，實施微分運算。
微分運算的泰勒近似，保留越多項，遞迴公式越精準。
微分運算的泰勒近似可以多取幾項，函數點可以多取幾個。
Truncation Error
「截斷誤差」。泰勒近似結果將消滅多餘的高次方項。
Consistency of Numerical Analysis
數值模擬當中的一致性：原本函數（微分方程式）、離散化函數（遞迴公式），兩者相符。
數值模擬當中的一致性誤差：原本函數（微分方程式）、離散化函數（遞迴公式），兩者差異。
一、函數值相符：考慮四捨五入誤差。
二、運算結果相符：加減乘除，考慮四捨五入誤差。微分積分，還得考慮截斷誤差。
數學家的說法：當取樣間隔Δx趨近零，原本函數與離散化函數的運算趨近相同。注意到，數學家只有討論「趨近零則趨近相同」、並未討論「越微小則越相同」。不太務實就是了。
數學家的說法：當取樣間隔Δx趨近零，截斷誤差趨近零。數學家只有考慮截斷誤差、並未考慮四捨五入誤差，但是指數穩定性可以消滅四捨五入誤差。
舉例來說，下述微分方程式：
d
—— f(t) = f(t)
dt
微分方程式重新整理成泛函數的模樣，才能相減得到差距。
differential equation:
d/dt f(t) = f(t)
d/dt f(t) - f(t) = 0
discretization:
(f(t) - f(t - Δt)) / Δt = f(t)
(f(t) - f(t - Δt)) / Δt - f(t) = 0
functional:
let L(f(t)) = d/dt f(t) - f(t)
let L̂(f(t)) = (f(t) - f(t - Δt)) / Δt - f(t)
consistency:
L(f(t)) = lim L̂(f(t))
Δt→0
lim { L(f(t)) - L̂(f(t)) } = 0
Δt→0
truncation error:
L(f(t)) - L̂(f(t))
大家習慣用大O符號呈現截斷誤差。
順帶一提，截斷誤差同時考慮時間變數與空間變數的步進，局部誤差只考慮時間變數的步進。有些文獻誤將截斷誤差與局部誤差視為相同，誤將O((Δt)ᵖ + (Δx)𐞥)與‖f⁽¹⁾ - f̂⁽¹⁾‖視為相同。
時間變數的步進，最簡單的方式是採用一階泰勒近似，p = 1。有些文獻因而省略時間變數的步進，並將截斷誤差與局部誤差視為相同，將O((Δx)𐞥)與‖f⁽¹⁾ - f̂⁽¹⁾‖視為相同。
Stability
「穩定性」。兩個數列距離受限。
穩定性細分為許多類型。
Lyapunov stability：內部穩定性。距離小於等於定值，定值取決於最初距離。
asymptotic stability：漸進穩定性。距離收斂至零。
exponential stability：指數穩定性。距離指數衰減。
https://en.wikipedia.org/wiki/Lyapunov_stability
Stability of Numerical Analysis
數值模擬當中的穩定性：數值解距離受限。
穩定性細分為許多類型。數值分析基本定理，採用收縮穩定性。
Lipschitz stability：平緩穩定性。距離小於等於最初距離的特定倍率。
contractive stability：收縮穩定性。距離小於等於最初距離。
https://en.wikipedia.org/wiki/Contraction_mapping
實務上容易檢查一致性。檢查遞迴公式是否符合泰勒近似即可。實務上很難檢查穩定性。實務上大家將收縮穩定性換成指數穩定性、換成更強更嚴格的性質。
(1) 線性遞迴公式，改寫成矩陣，特徵值絕對值皆小於等於1。
(2) 線性遞迴公式，改寫成特徵方程式，根絕對值皆小於等於1。
指數穩定性又細分為許多類型。數學家做證明會用到。
weak stability：弱穩定性。特徵值／根絕對值等於1有兩個以上。
strong stability：強穩定性。特徵值／根絕對值等於1至多一個。
zero-stability：零穩定性。沒有多重特徵值／重根。
https://sites.math.rutgers.edu/~falk/math573/lecture21.pdf
補充說明一下。Banach Fixed-point Theorem當中，收縮穩定性只能是小於。數值分析基本定理當中，收縮穩定性可以追加等於。第零回合，令數值解趨近符號解、令初始距離趨近零；往後回合，即便距離不變也無妨。一開始就夾緊，往後保持夾緊，依然收斂。
Contractive Stability（B-stability）
「收縮穩定性」。數值解距離小於等於最初距離。
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ contractive stability
當遞迴公式不因時刻而變，收縮穩定性變成：數值解每回合距離變小或不變。
‖f̂⁽ⁿ⁺¹⁾ - ĝ⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ time-invariant
當遞迴公式恰是線性遞迴函數，兩種數值解相加減仍是數值解。收縮穩定性變成：數值解每回合範數變小或不變。
‖f̂⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾‖ linearity
Exponential Stability（A-stability）
「指數穩定性」。數值解距離指數衰減。
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ exp(-c n) ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ exponential stability
c is a constant and 0 ≤ exp(-c n) ≤ 1
當遞迴公式不因時刻而變，指數穩定性變成：數值解每回合距離乘上固定倍率，其小於等於1。
‖f̂⁽ⁿ⁺¹⁾ - ĝ⁽ⁿ⁺¹⁾‖ ≤ k ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ time-invariant
k is a constant and 0 ≤ k ≤ 1
當遞迴公式恰是線性遞迴函數，兩種數值解相加減仍是數值解。指數穩定性變成：數值解每回合範數乘上固定倍率，其小於等於1。
‖f̂⁽ⁿ⁺¹⁾‖ ≤ k ‖f̂⁽ⁿ⁾‖ linearity
k is a constant and 0 ≤ k ≤ 1
Stability Analysis with Matrix
收縮穩定性：數值解距離小於等於最初距離。
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ contractive stability
當遞迴公式恰是線性遞迴函數，得以簡化成矩陣運算。
‖Aⁿf̂⁽⁰⁾ - Aⁿĝ⁽⁰⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ linear scheme
當遞迴公式恰是線性遞迴函數，兩種數值解相加減仍是數值解。收縮穩定性變成：數值解每回合範數變小或不變。
‖Aⁿf̂⁽⁰⁾‖ ≤ ‖f̂⁽⁰⁾‖ linearity
另一方面，根據範數的性質，得到不等式。
‖Aⁿf̂⁽⁰⁾‖ ≤ ‖Aⁿ‖ ‖f̂⁽⁰⁾‖ ≤ ‖A‖ⁿ ‖f̂⁽⁰⁾‖ sub-multiplicative
實務上很難檢查收縮穩定性。一種替代方案是指數穩定性：令矩陣範數受限，導致矩陣範數小於等於1。
let ‖A‖ⁿ ‖f̂⁽⁰⁾‖ ≤ ‖f̂⁽⁰⁾‖ for all n exponential stability
thus ‖A‖ⁿ ≤ 1 for all n exponential stability
thus ‖A‖ ≤ 1 exponential stability
矩陣範數採用譜範數：矩陣特徵值絕對值最大者小於等於1。矩陣特徵值皆介於±1（若是實數）。
‖A‖₂ = max(abs(eigval(A))) ≤ 1 spectral norm
0 ≤ abs(eigval(A)) ≤ 1 inside the unit circle
-1 ≤ eigval(A) ≤ 1 when eigval(A) are real
想要檢查指數穩定性，只需將遞迴公式改寫成矩陣，然後檢查所有特徵值絕對值皆小於等於1。
0 ≤ abs(eigval(A)) ≤ 1 => exponential stability
Contractive Stability ⇔ Exponential Stability
當遞迴公式恰是線性遞迴函數，收縮穩定性等價於指數穩定性。
for linear scheme,
contractive stability <=> exponential stability
援引三角不等式，即可證明。
‖f̂⁽ⁿ⁺¹⁾‖ = ‖Af̂⁽ⁿ⁾‖ ≤ ‖A‖ ‖f̂⁽ⁿ⁾‖ sub-multiplicative
‖f̂⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾‖ contractive stability
‖A‖ ≤ 1 exponential stability
常見的範數類型是L₂-norm。矩陣範數採用譜範數：特徵值絕對值最大者。向量範數採用平方範數：元素平方和開根號。
‖f̂⁽ⁿ⁺¹⁾‖₂ = ‖Af̂⁽ⁿ⁾‖₂ ≤ ‖A‖₂ ‖f̂⁽ⁿ⁾‖₂ sub-multiplicative
‖f̂⁽ⁿ⁺¹⁾‖₂ ≤ ‖f̂⁽ⁿ⁾‖₂ contractive stability
‖A‖₂ ≤ 1 exponential stability
不斷乘以一個小於等於1的倍率，想當然變小或不變。倍率從數字換成矩陣，依然如此。就這樣。
順帶一提，指數穩定性得以消滅或抑制四捨五入誤差。數值解變小或不變，四捨五入誤差也隨之變小或不變。
Stability Analysis with Characteristic Equation
矩陣特徵值改成特徵方程式的根。總之都是頻譜。
linear difference equation
solution f(t) = c1 λ1^t + ... cn λn^t
想要檢查指數穩定性，只需將遞迴公式改寫成特徵方程式，然後檢查所有根絕對值皆小於等於1。
Algebraic Stability
數值分析基本定理，要求Δt與Δx趨近零。然而浮點數位數有限，不可能趨近零。實務上必須檢查特定大小的Δt與Δx是否滿足穩定性。嚴謹起見，數學家另造一詞，稱作代數穩定性。簡單起見，大家還是稱作穩定性。
代數穩定性細分許多類型。數學家做證明會用到。
https://sites.math.rutgers.edu/~falk/math573/lecture21.pdf
http://minitorn.tlu.ee/~jaagup/uk/dynsys/ds2/num/Absolute/Absolute.html
absolute stability：絕對穩定性。特定大小的Δt與Δx，滿足指數穩定性。
relative stability：相對穩定性。當Δt與Δx趨近零，可以滿足指數穩定性。
絕對穩定性又細分為許多類型。數學家做證明會用到。
https://dl.acm.org/doi/pdf/10.1145/322139.322147
https://en.wikipedia.org/wiki/L-stability
A-stability：特定大小的Δt與Δx，滿足指數穩定性。
L-stability：承上，當Δt與Δx趨近無限大，指數趨近負無限大。
Interval of Absolute Stability
當遞迴公式恰是線性遞迴函數，特徵方程式的根皆介於±1，才會滿足指數穩定性。以此推導Δt與Δx的適當範圍。
舉例來說，下述微分方程式：
d
—— f(t) = f(t)
dt
套用Backward Euler Method，選擇適當的Δt，就會穩定。
(f(t) - f(t - Δt)) / Δt = f(t) 1st-order backward derivative
(f[i] - f[i-1] ) / Δt = f[i] implicit method
f[i] = f[i-1] / (1 - Δt) recurrence
f[i] = a f[i-1] where a = 1 / (1 - Δt)
f[i] is stable iff |a| ≤ 1
f[i] is stable iff Δt ≤ 0 and Δt ≥ 2
f[i] is stable when Δt ≥ 2
套用Forward Euler Method，無論如何不可能穩定。
(f(t + Δt) - f(t)) / Δt = f(t) 1st-order forward derivative
(f[i+1] - f[i]) / Δt = f[i] explicit method
f[i+1] = f[i] (1 + Δt) recurrence
f[i+1] = a f[i] where a = (1 + Δt)
f[i+1] is stable iff |a| ≤ 1
f[i+1] is stable iff Δt ≥ -2 and Δt ≤ 0
f[i+1] is unstable since Δt > 0
Function of Absolute Stability
數學家將a改寫成函數，叫做Function of Absolute Stability，簡稱Stability Function。
舉例來說，下述微分方程式：
d
—— f(t) = f(t)
dt
套用Backward Euler Method，將數值a換成函數ϕ(Δt)。
f[i+1] = a f[i] where a = 1 / (1 - Δt)
f[i+1] = ϕ(Δt) f[i] where ϕ(Δt) = 1 / (1 - Δt)
數學家甚至追加變數。讓微分方程式係數是變數，讓遞迴公式係數是變數，方便判斷各種係數的代數穩定性，方便調整係數。
舉例來說，下述微分方程式：
d
—— f(t) = k f(t) k is a coefficient
dt
套用Backward Euler Method，將數值a換成函數ϕ(k Δt)。
f[i+1] = a f[i] where a = 1 / (1 - k Δt)
f[i+1] = ϕ(k Δt) f[i] where ϕ(k Δt) = 1 / (1 - k Δt)
k可以是複數，形成複數函數ϕ(z)。函數輸入是z = k Δt。
ϕ(z) = 1 / (1 - z)
數學家描述等高線，習慣使用符號ϕ。我也採用符號ϕ。
Region of Absolute Stability
有些教科書喜歡畫一種圖，叫做Region of Absolute Stability，簡稱Stable Region。此圖毫無實際用處，只是想讓讀者轉換心情。
舉例來說，下述微分方程式：
d
—— f(t) = k f(t) k is a coefficient
dt
套用Backward Euler Method。穩定區域|a| ≤ 1即是|ϕ(k Δt)| ≤ 1畫在複數平面。座標軸是(k Δt)的實部虛部。得到左圖。
套用Forward Euler Method，得到右圖。
【圖片尚待修正】
k可能是複數，於是畫在複數平面、視作二維向量。Δt只能是實數，於是視作向量縮放倍率。穩定區域|a| ≤ 1與射線k的交點，觀察縮放倍率，得到適當的Δt。
數學家利用ϕ(z)定義代數穩定性。此例當中：
A-stability：原點放射線可以碰到穩定區域、且不是碰到原點。
L-stability：承上，當ϕ(z)趨近零，z趨近負無限大。
Characteristic Equation
觀察遞迴公式係數，判斷指數穩定性。
舉例來說，下述微分方程式：
d
—— f(t) = k f(t) k is a coefficient
dt
解析解：
exp(k t)介於±1，(k t)是負數或零，才會滿足指數穩定性。
f(t) = a exp(k t) + c(t) 初始值咧
套用Backward Euler Method，指數穩定性變成：特徵方程式之分式是負數或零（實部為負或零）。
套用Forward Euler Method，無論如何不可能穩定。
方程式 d/dt f(t) = F(f(t))
遞迴公式 (ak f(t+k) + ... + a0 f(t)) = Δt (bk F(f(t+k)) + ... + b0 F(f(t)))
特徵方程式 a(λ) - Δt b(λ) F(λ) = 0
微分方程式 d/dt f(t) = k f(t)
解析解 x = a exp(k t) + c(t)
遞迴公式 (ak f(t+k) + ... + a0 f(t)) = k Δt (bk f(t+k) + ... + b0 f(t))
特徵方程式 a(λ) = k Δt b(λ)
特徵方程式 k Δt = a(λ)/b(λ)
特徵多項式的根 λ
指數穩定性 |λ| ≤ 1
指數穩定性 Re[k Δt] ≤ 0
指數穩定性 Re[a(λ)/b(λ)] ≤ 0
https://www.math.unipd.it/~alvise/AN_2017/LETTURE/DAHLQUIST_STAB.pdf
http://cds.cern.ch/record/1069162
https://www.jstor.org/stable/2156562
https://www.jstor.org/stable/2156895
https://encyclopediaofmath.org/wiki/Non-linear_stability_of_numerical_methods
複數根λ，形成複數函數q(λ)。函數輸出是q(λ) = k Δt。
q(λ) = a(λ)/b(λ) = k Δt
數學家利用q(λ)定義代數穩定性。此例當中：
A-stability：q(λ)完全位於左半複數平面。
L-stability：承上，當q(λ)趨近無限大，λ趨近零。
數學家描述除法商數，習慣使用符號q。我也採用符號q。
Stability Analysis with Fourier Transform
使用傅立葉轉換。總之都是頻譜。
Stiffness
微分方程式的造型太過特殊，選擇很小很小的Δt與Δx才能穩定，需要很多很多的遞推回合才能到達指定時刻，導致很久很久的執行時間才能跑完數值模擬。數學家稱作「剛硬性」。硬篤、歹舞。
剛硬性沒有明確數學定義。只是一種感覺。
Numerical Scheme
數值方案。面對各式各樣的微分方程式，想方設法調整遞迴公式，導致數值收斂。
具體來說，滿足一致性、滿足穩定性。
consistent numerical scheme：設定遞迴公式的係數，滿足一致性。
stable numerical scheme：設定數值解的數學性質，滿足穩定性。
關於一致性，分別討論時間差分與空間差分。
spatial finite difference：對空間變數偏微分的離散化
temporal finite difference：對時間變數偏微分的離散化
linear numerical scheme：遞迴公式是線性遞迴函數
線性微分方程式＋線性時間差分＋線性空間差分
nonlinear numerical scheme：遞迴公式是非線性遞迴函數
關於穩定性，分別討論各式各樣的數值解特性。
(1) contractive scheme：數值解距離遞減／數值解範數遞減
(2) monotone scheme：數值解相對大小不變／數值解是單調函數
(3) monotonic upstream-centered scheme：數值解是單峰函數
為了應付非線性微分方程式，為了應付時間變化空間變化無法分離成左式右式，花招百出。
Spatial Finite Difference
(1) 決定數值。
(1a) implicit method：取新值。
(1b) explicit method：取舊值。
(1c) implicit–explicit method：有些取新值、有些取舊值。
(1d) Crank–Nicolson method：新值加舊值然後除以二。雞肋。
(2) operator splitting：分割步驟。
(3) semi-implicit method (symplectic method)：更動步驟。
implicit。取新值。需要移項整理，甚至需要一次方程組求解。
d/dx f(x) = f(x) + 1 ---> (f[i] - f[i-1]) / Δx = f[i] + 1
^^^^
explicit。取舊值。
d/dx f(x) = f(x) + 1 ---> (f[i] - f[i-1]) / Δx = f[i-1] + 1
^^^^^^
operator splitting。項次太多，逐步更新。
d/dx f(x) = 2f(x) + 3g(x) ---> 1. (f[i] - f[i-1]) / Δx = 2 f[i]
2. (f[i] - f[i-1]) / Δx = 3 g[i]
semi-implicit。式子太多，設定順序。
d/dx f(x) = 2f(x) + 3g(x) ---> 1. (f[i] - f[i-1]) / Δx = 3 g[i]
2. (f[i] - f[i-1]) / Δx = 2 f[i]
Stencil
Temporal Finite Difference
(1) 取左方斜率、取右方斜率。
(2) Multistep Method：取多處斜率。
(3) Multistage Method：求多次斜率。
(4) Multiderivative Method：取多階微分。
Adams method: Maclaurin series
https://mathworld.wolfram.com/AdamsMethod.html
Richardson Extrapolation
https://en.wikipedia.org/wiki/Richardson_extrapolation
取左方斜率，無論是一處或多處，無論是遞增或遞推，無論是顯式或隱式，只要強行展開數學式子，多階段就會變成單階段，缺乏討論意義。取右方斜率，才有討論意義。
重新歸納為左右兩種。宛如差分。
Euler Method：取當地斜率。參考先前當前函數點。
Multistep Method：取左方斜率。參考先前幾個函數點。優點是快。
Multistage Method：取右方斜率。預測之後幾個函數點。優點是準。
左右都取：同時獲得兩者缺點，沒有人這樣做。
Butcher Tableau
Stable Numerical Scheme
針對穩定性：
Lax Equivalence Theorem：
IVP初始值問題、一次微分方程式、線性方案、指數穩定性
Dahlquist Equivalence Theorem：
IVP初始值問題、非一次微分方程式、線性方案、Multistep Method、零穩定性
針對代數穩定性：
Courant–Friedrichs–Lewy Condition：
IVP初始值問題、物理學之連續方程式、線性方案、CFL condition
Strong-Stability Preserving Property：
IVP初始值問題、物理學之連續方程式、Euler Method家族系列、CFL condition
參考資料。
Hyperbolic partial differential equations
Peter D. Lax
線性遞迴公式，數值解，冪級數不超過一
Difference methods for initial-value problems
Robert D. Richtmyer, K. W. Morton.
線性遞迴公式，數值解，頻譜強度不超過一
應該就是Von Neumann Stability Analysis
Numerical Solution of Partial Differential Equations: An Introduction
K. W. Morton, D. F. Mayers
線性遞迴公式，數值解距離，冪級數收斂
An Introduction to Numerical Analysis
Endre Süli, David F. Mayers
一般遞迴公式，數值解距離，冪級數收斂
Numerical solution of partial differential equations: finite difference methods
Gordon D. Smith
線性遞迴公式，數值解，最大的奇異值不超過一
Constractive Scheme
線性方案當中，收縮穩定性等同指數穩定性。
for linear scheme,
contractive stability <=> exponential stability
contractive stability 數值解範數每回合變小或不變
exponential stability 數值解距離每回合呈指數衰減
解析解不斷變大、或者反覆變大變小，不適用指數穩定性。線性遞迴公式的解析解應該不會反覆變大變小。【尚待確認】
彈簧振動、琴弦振動、鼓膜振動的解析解，可以寫成複數的實部。儘管實部乍看反覆變大變小，但是複數其實保持變小或不變，滿足指數穩定性。【尚待確認】
Algebraic Stability Property
「代數穩定性性質」。
套用隱式遞迴公式，指數穩定性變成：特徵方程式分式是負數或零（實部為負或零）。
套用顯式遞迴公式，無論如何不可能穩定。
我要找一下Dahlquist的書看看有沒有證明。
Courant–Friedrichs–Lewy Condition
物理學之連續方程式，一個線性方案，一旦滿足指數穩定性，也會連帶滿足某些特殊性質。此處要談的特殊條件是「CFL條件」。
for linear scheme of continuity equation,
contractive stability => CFL condition
連續方程式請見後面章節。連續方程式的某些衍生版本（例如波動方程式、熱傳導方程式）符合這個結論，但是無法保證各種衍生版本都符合這個結論。這個結論不是推理而得，而是觀察而得。因此這個結論不被認為是數學定理。
CFL條件的數學式子：波動方程式是u Δt / Δx ≤ 1。熱傳導方程式是2 u ∆t / (∆x)² ≤ 1。視情況做調整。其中u是速度。
證明請見講義：
【待補範例】
CFL條件的幾何意義：移動距離u Δt不得超過取樣距離Δx，否則空間差分會抓錯樣本。
空間差分必須抓到相鄰樣本。當移動距離太大，則空間差分誤抓到從遠處跨越鄰居而來的樣本。
CFL 條件是必要條件：滿足指數穩定性，導致滿足CFL條件。CFL條件不是充分條件：滿足CFL條件，可能穩定、可能不穩定。
for linear scheme of continuity equation,
contractive stability ⇍ CFL condition
CFL條件的主要功用是快篩試劑。逆否命題：不滿足CFL條件，導致不滿足指數穩定性。
Strong-Stability Preserving Property
如果一個數值方案，時間差分採用Forward Euler Method，一旦滿足CFL Condition，並且導致指數穩定性（或者特定的穩定性，例如單調穩定性），那麼該數值方案，時間差分採用Forward Euler Method家族系列演算法，一旦滿足CFL Condition乘上適當倍率，也會導致指數穩定性（或者特定的穩定性）。
「保強穩定性性質」原理簡單。Foward Euler Method只取一處斜率。Foward Euler Method家族系列，則是取多處斜率的加權平均數。視作多個Foward Euler Method各自加權，而CFL Condition也跟著加權。
證明請見專著《Strong Stability Preserving Runge–Kutta and Multistep Time Discretizations》。
【待補範例】
「保強穩定性性質」名稱很雲。此概念既非強穩定性、亦非保Ｘ變換。此概念是將時間差分從一階泰勒展開推廣成Euler Method家族系列，在採用CFL Condition的情況下。
「保強穩定性性質」內容很雲。CFL Condition只是必要條件，此性質卻以充分條件作為前提。前提是一件未必成立的事情。
順帶一提，仔細修改時間差分的係數，讓適當倍率是1，就能去除「乘上適當倍率」六個字。
Monotone Scheme
「單調方案」。數值解每回合維持高低關係。
順帶一提，數學家命名方案，有時候按照遞迴公式的格式，有時候按照數值解的性質。線性方案屬於前者。單調方案屬於後者。
實務上很難檢查單調性。實務上大家將單調性換成保單調性、換成更強更嚴格的性質。
根據Harten's theorem，單調性可以換成更強更嚴格的性質。
for linear scheme,
monotone => TVD => monotonicity preserving
monotonicity preserving 數值解每回合維持單調函數
total variation diminishing 數值解每回合總差值減少或不變
monotone 兩個數值解每回合維持高低關係
根據Godunov's theorem，單調方案只可能是一階空間差分。
「保單調方案」用來消除Runge Phenomenon。單調函數不會上下抖動。
然而大部分的微分方程式暨初始條件，解通常不是單調函數。保單調方案通常無法派上用場。
一種方式是空間分段，空間小段滿足單調函數。
Monotonic Upstream-centered Scheme
「單調逆流中心方案」。
https://en.wikipedia.org/wiki/Godunov's_scheme
https://en.wikipedia.org/wiki/MUSCL_scheme
https://en.wikipedia.org/wiki/Shock-capturing_method
