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 is an 1D interval
space is a 0D point
=> boundary is two endpoints
f(t,x) 5 ≤ t ≤ 10 time is an 1D interval
-10 ≤ x ≤ +10 space is an 1D interval
=> boundary is a rectangle
f(t,x,y) 5 ≤ t ≤ 10 time is an 1D interval
x² + y² ≤ 10 space is a 2D disk
=> boundary is a cylinder
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
數值收斂。數值遞推過程，全域誤差趨近零，答案精準。
中文譯作「收斂」，原意應是「合流」。
數值收斂的數學定理的進化過程：
(1) problem: series convergence
algorithm: Cauchy condensation test
convergence: monotone convergence theorem
(if monotone, then bounded = convergent)
(2) problem: equation solving
algorithm: fixed point iteration
convergence: Banach fixed-point theorem
(if contractive, then convergent)
(3) problem: differential equation solving
algorithm: Euler method
convergence: fundamental theorem of numerical analysis
(if consistent, then stable = convergent)
問題：級數收斂
演算法：柯西收縮測試
收斂性：單調收斂定理（單調時，受限即收斂。）
問題：方程式求解
演算法：不動點遞推法
收斂性：巴拿赫不動點定理（收縮則收歛。）
問題：微分方程式求解
演算法：歐拉法
收斂性：數值分析基本定理（一致時，穩定即收斂。）
fundamental theorem of numerical analysis
「數值分析基本定理」。
if consistent, then stable = convergent
拆成兩句話。第一句話是本文主角。第二句話顯而易見。
(1) consistent & stable => consistent & convergent
(2) consistent & convergent => consistent & stable (trivial)
這三個詞彙的中文翻譯：
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 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
L(f(t,x)) = 0
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̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ when n≥0
with consistency,
distance of corresponding iteration of
a symbolic solution and a numeric solution with
same initial condition is bounded.
‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ ≤ ‖f⁽⁰⁾ - f̂⁽⁰⁾‖ when n≥0
convergence: distance of corresponding iteration of
two numeric solutions with
different initial conditions leads to zero.
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ → 0 when n→∞
with consistency,
distance of corresponding iteration of
a symbolic solution and a numeric solution with
same initial condition leads to zero.
‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ → 0 when n≥0
一致性：函數方程式、其離散化（遞迴公式），兩者趨近相等。
穩定性：兩組數值解（初始條件不同），各回合距離受限。
一致性且穩定性：符號解、數值解，各回合距離受限。
收斂性：兩組數值解（初始條件不同），最終趨近相等。
一致性且收斂性：符號解、數值解，各回合趨近相等。
一致性：截斷誤差趨近零。
一致性且穩定性：全域誤差受限。
一致性且收斂性：全域誤差趨近零。
大部分文獻又進一步修改定義，將f(t,x)改成d/dt f(t,x)再改成微分方程式的右式F(f(t,x))。我推測原因是避免在Δt→0的情況下乘以除以Δt，並且確保時間變化空間變化可以分離成左式右式。
如此一來，穩定性宛如平緩函數Lipschitz function，數值分析基本定理宛如不動點定理Banach fixed-point theorem。平緩函數套用不動點遞推法收斂至不動點。請見本站文件「fixed point 」。
differential equation
d/dt f(t,x) = F(f(t,x))
symbolic solution
f(t,x)
consistency
F(f(t,x)) - F̂(f(t,x)) → 0 when Δt→0, Δx→0
stability with consistency
‖F(f⁽ⁿ⁾) - F(f̂⁽ⁿ⁾)‖ ≤ ‖f⁽⁰⁾ - f̂⁽⁰⁾‖ when n≥0
convergence with consistency
‖F(f⁽ⁿ⁾) - F(f̂⁽ⁿ⁾)‖ → 0 when n≥0
時間變化空間變化可以分離成左式右式，遞迴公式不受時間影響，每回合一律使用相同泛函數F̂(f(t,x))，事情比較簡單。時間變化空間變化無法分離成左式右式，遞迴公式受時間影響，每回合都要微調泛函數F̂ₜ(f(t,x))，事情更加棘手，我沒有研究。
數值分析基本定理的證明，我還沒有學會。等我學完泛函分析不知何年何月。關鍵字也許是Banach space、uniform boundedness principle、dominated convergence theorem。
https://math.stackexchange.com/questions/4460567/
https://www.colorado.edu/amath/sites/default/files/attached-files/lax_equivalence_theorem.pdf
https://users.soe.ucsc.edu/~hongwang/AM213B/Notes/Lecture02.pdf
https://users.soe.ucsc.edu/~hongwang/AM213B/Notes/Lecture03.pdf
補充說明一下。不動點定理當中，平緩函數只能是小於。數值分析基本定理當中，穩定性可以追加等於。第零回合，令數值解趨近符號解、令初始距離趨近零；往後回合，即便距離不變也無妨。一開始就夾緊，往後保持夾緊，依然收斂。
補充說明一下。數值分析基本定理當中，只有設置初始條件，沒有設置邊界條件。符號解和數值解，左右延伸至無限遠。
總之，當一致性和穩定性同時成立，就能保證數值收斂。
接下來分別介紹一致性和穩定性。
consistency與consistency error
「一致性」。符號轉數值，符號與數值相符。分為兩點討論。
一、數字大小相符。
二、運算結果相符。加減乘除微積求值代入。
「一致性誤差」。符號轉數值導致的誤差。分為兩點討論。
round-off error：四捨五入誤差。實數四則運算從符號轉數值導致的誤差。
truncation error：截斷誤差。函數微分運算從符號轉數值導致的誤差。
實數四則運算從符號轉數值，一般採用浮點數運算。四捨五入誤差，一般是指浮點數運算誤差。請見本站文件「bit 」。
函數微分運算從符號轉數值，一般採用泰勒近似。截斷誤差，一般是指泰勒近似誤差。請見本站文件「approximation 」。
floating-point arithmetic與round-off error
「浮點數運算」。離散化實數之四則運算。
實數資料結構採用浮點數，實施四則運算。
「四捨五入誤差」。浮點數運算總是消滅多餘的低位數。
本節不是本文重點。細節請自行上網搜尋講義。
finite difference與truncation error
「有限差分」。離散化實函數之微分運算。
實函數資料結構採用regular discretization，或者其他資料結構，實施微分運算。
「截斷誤差」。泰勒近似總是消滅多餘的高次方項。
一次微分的差分公式。細分成左邊、右邊、中央三種版本。
Taylor approximation
(1) f(x+Δx) = f(x) + f′(x)Δx + ½f″(x)(Δx)² + ⅙f‴(x)(Δx)³ + ...
(2) f(x-Δx) = f(x) - f′(x)Δx + ½f″(x)(Δx)² - ⅙f‴(x)(Δx)³ + ...
first-order backward difference
f(x-Δx) = f(x) - f′(x)Δx + O((Δx)²)
f(x) - f(x-Δx) = f′(x)Δx + O((Δx)²)
(f(x) - f(x-Δx)) / Δx = f′(x) + O(Δx)
f′(x) = (f(x) - f(x-Δx)) / Δx + O(Δx)
f′[i] = (f[i] - f[i-1]) / Δx + O(Δx)
first-order forward difference
f(x+Δx) = f(x) + f′(x)Δx + O((Δx)²)
f(x+Δx) - f(x) = f′(x)Δx + O((Δx)²)
(f(x+Δx) - f(x)) / Δx = f′(x) + O(Δx)
f′(x) = (f(x+Δx) - f(x)) / Δx + O(Δx)
f′[i] = (f[i+1] - f[i]) / Δx + O(Δx)
second-order central difference
⎰ f(x+Δx) = f(x) + f′(x)Δx + ½f″(x)(Δx)² + O((Δx)³)
⎱ f(x-Δx) = f(x) - f′(x)Δx + ½f″(x)(Δx)² + O((Δx)³)
f(x+Δx) - f(x-Δx) = 2 f′(x)Δx + O((Δx)³)
(f(x+Δx) - f(x-Δx)) / (2 Δx) = f′(x) + O((Δx)²)
f′(x) = (f(x+Δx) - f(x-Δx)) / (2 Δx) + O((Δx)²)
f′[x] = (f[i+1] - f[i-1]) / (2 Δx) + O((Δx)²)
大家習慣使用大O符號呈現截斷誤差。大家習慣用Δx次方值作為截斷誤差的數量級。n階差分的n就是截斷誤差的數量級。
Δx趨近於零、Δx遠小於一。Δx次方越高，截斷誤差越小。泰勒展開低次方項保留越多，截斷誤差越小。
想讓截斷誤差趨近零，一種方式是採用平滑函數。如此一來，泰勒展開高次方項通通存在、而且通通趨近零。
想讓截斷誤差趨近零，也可以採用其他函數。找到截斷誤差的數學式子，證明其趨近零。比較麻煩就是了。
負面詞彙「誤差error」、正面詞彙「精度accuracy」，兩者概念相同，都有人用。
二次微分的差分公式。
second-order central difference of second derivative
⎰ f(x+Δx) = f(x) + f′(x)Δx + ½f″(x)(Δx)² + ⅙f‴(x)(Δx)³ + O((Δx)⁴)
⎱ f(x-Δx) = f(x) - f′(x)Δx + ½f″(x)(Δx)² - ⅙f‴(x)(Δx)³ + O((Δx)⁴)
f(x+Δx) + f(x-Δx) = 2 f(x) + f″(x)(Δx)² + O((Δx)⁴)
f(x+Δx) - 2 f(x) + f(x-Δx) = f″(x)(Δx)² + O((Δx)⁴)
(f(x+Δx) - 2 f(x) + f(x-Δx)) / (Δx)² = f″(x) + O((Δx)²)
f″(x) = (f(x+Δx) - 2 f(x) + f(x-Δx)) / (Δx)² + O((Δx)²)
f″[i] = (f[i+1] - 2 f[i] + f[i-1]) / (Δx)² + O((Δx)²)
一次微分的差分公式，增加點數以增加精度。我沒有研究。
4th-order central difference
f′[x] = (- f[i+2] + 8 f[i+1] - 8 f[i-1] + f[i-2]) / (12 Δx)
+ O((Δx)⁴)
https://en.wikipedia.org/wiki/Finite_difference_coefficient
一些特殊技巧。我沒有研究。
compact finite difference：用更少的點數達到更高的精度。
degraded finite difference：達到更低的精度。自虐。
https://en.wikipedia.org/wiki/Compact_finite_difference
https://people.bath.ac.uk/ensdasr/COMPACT/dasr.compact.pdf
想要迅速獲得中央差分公式，除了泰勒展開以外，也可以使用Lagrange interpolation然後微分。想要簡化差分公式並且維持精度，可以使用Padé approximation。我沒有研究。
consistency of numerical simulation
數值模擬當中的一致性：當取樣間隔Δx趨近零，原本函數（微分方程式）、離散化函數（遞迴公式），兩者相符。
注意到，數學家只有討論「Δ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 + O((Δt)ᵖ) = f(t)
(f(t) - f(t - Δt)) / Δt - f(t) = O((Δt)ᵖ)
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((Δt)ᵖ)
數學家只有考慮截斷誤差，並未考慮四捨五入誤差。馬馬虎虎就是了。
微分方程式的各個微分運算的截斷誤差，組成了局部誤差。截斷誤差趨近零，那麼局部誤差也趨近零。想讓截斷誤差趨近零，一種方式是符號解f(t)採用平滑函數。
順帶一提，截斷誤差同時考慮時間變數與空間變數的步進，局部誤差只考慮時間變數的步進。有些文獻誤將截斷誤差與局部誤差視為相同，誤將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 simulation
數值模擬當中的穩定性：數值解距離受限。
穩定性細分為許多類型。
bounded stability：受限穩定性。距離小於等於定值。
Lipschitz stability：平緩穩定性。距離小於等於最初距離的特定倍率。
contractive stability：收縮穩定性。距離小於等於最初距離。
https://en.wikipedia.org/wiki/Contraction_mapping
數值分析基本定理，理論上是受限穩定性，實務上是收縮穩定性。受限穩定性很難判斷，大家換成更強更嚴格的收縮穩定性。
順帶一提，收縮穩定性得以消滅或抑制數值誤差。數值解變小或不變，數值誤差也隨之變小或不變。
contractive stability
「收縮穩定性」。數值解距離小於等於最初距離。
‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ contractive stability
當遞迴公式不因時刻而變，收縮穩定性變成：數值解每回合距離變小或不變。
‖f̂⁽ⁿ⁺¹⁾ - ĝ⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ time-invariant
當遞迴公式恰是線性遞迴函數，兩種數值解相加減仍是數值解。收縮穩定性變成：數值解每回合範數變小或不變。
‖f̂⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾‖ linearity
exponential 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
當遞迴公式恰是線性遞迴函數，收縮穩定性改寫成矩陣運算。
實際範例請見本站文件「linear recurrence 」。
‖f̂⁽ⁿ⁺¹⁾‖ = ‖Af̂⁽ⁿ⁾‖ ≤ ‖f̂⁽ⁿ⁾‖ contractive stability
當遞迴公式恰是線性遞迴函數，指數穩定性改寫成矩陣運算。
‖f̂⁽ⁿ⁺¹⁾‖ = ‖Af̂⁽ⁿ⁾‖ ≤ k ‖f̂⁽ⁿ⁾‖ exponential stability
k is a constant and 0 ≤ k ≤ 1
另一方面，根據範數的性質，得到不等式。
‖Af̂⁽ⁿ⁾‖ ≤ ‖A‖ ‖f̂⁽ⁿ⁾‖ sub-multiplicative
根據不等式，收縮穩定性等價於矩陣範數小於等於1。
矩陣範數視作倍率k。收縮穩定性等價於指數穩定性！
不斷乘以一個小於等於1的倍率，顯然指數衰減。倍率從數字換成矩陣，依然指數衰減。就這樣。
‖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。矩陣特徵值皆介於±1（若是實數）。
‖A‖₂ = max(abs(eigval(A))) ≤ 1 spectral norm
0 ≤ abs(eigval(A)) ≤ 1 inside the unit circle
-1 ≤ eigval(A) ≤ 1 if eigval(A) are real
想要檢查指數穩定性，只需將線性遞迴公式改寫成矩陣，然後檢查所有特徵值絕對值皆小於等於1。
0 ≤ abs(eigval(A)) ≤ 1 exponential stability
順帶一提，指數穩定性顯然滿足收斂性，畢竟指數衰減至零。大可不必煩惱數值分析基本定理。
absolute stability
數值分析基本定理，要求Δt與Δx趨近零。然而浮點數位數有限，況且陣列大小有限，因此Δt與Δx不可能趨近零。實務上必須檢查特定大小的Δt與Δx是否滿足穩定性。嚴謹起見，數學家另造一詞，稱作絕對穩定性。簡單起見，大家還是稱作穩定性。
穩定性細分許多類型。
relative stability：相對穩定性。當Δt與Δx趨近零，滿足穩定性。
absolute stability：絕對穩定性。特定大小的Δt與Δx，滿足穩定性。
stiff stability：剛硬穩定性。滿足絕對穩定性，而且，
當Δt與Δx趨近無限大，
則收斂速度趨近無限大（指數趨近負無限大）。
stiffness
「剛硬性」。微分方程式的造型太過特殊，選擇很小很小的Δ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) backward difference
(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)
a numerical scheme is stable iff |a| ≤ 1
this numerical scheme is stable iff Δt ≤ 0 and Δt ≥ 2
this numerical scheme is stable when Δt ≥ 2
套用forward Euler method，無論如何不可能穩定。
(f(t + Δt) - f(t)) / Δt = f(t) forward difference
(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)
a numerical scheme is stable iff |a| ≤ 1
this numerical scheme is stable iff Δt ≥ -2 and Δt ≤ 0
this numerical scheme 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，得到左圖。
套用forward Euler method，得到右圖。
穩定區域|a| ≤ 1即是|ϕ(k Δt)| ≤ 1畫在複數平面。座標軸是(k Δt)的實部虛部。
【圖片尚待修正】
k可能是複數，於是畫在複數平面、視作二維向量。Δt只能是實數，於是視作向量縮放倍率。穩定區域|a| ≤ 1與射線k的交點，觀察縮放倍率，得到適當的Δt。
數學家利用ϕ(z)定義穩定性。此例當中：
relative stability：穩定區域包含原點。
absolute stability：原點放射線可以碰到穩定區域、且不是碰到原點。
stiff stability：承上，當ϕ(z)趨近零，z趨近負無限大。
stability analysis with Fourier transform
專著《Computational Aerodynamics》。
舉例來說，下述微分方程式：
d d
—— f(t,x) + k —— f(t,x) = 0 k is a coefficient
dt dx
遞迴公式：
f[i+1][j] = f[i][j] - ½ C (f[i][j+1] - f[i][j-1])
where C = k Δt / Δx is Courant number
forward difference of t
central difference of x
傅立葉轉換：
𝓕(f[i+1]) = 𝓕(f[i]) - ½ C (𝓕(f[i]) exp(+𝑖wΔx) - 𝓕(f[i]) exp(-𝑖wΔx))
𝓕(f[i+1]) = (1 - ½ C (exp(+𝑖wΔx) + exp(-𝑖wΔx)) 𝓕(f[i])
無論如何不可能穩定。
g(w) = 1 - ½ C (exp(𝑖wΔx) - exp(-𝑖wΔx)) amplification factor
g(w) = 1 - C 𝑖 sin(wΔx) amplification factor
a numerical scheme is stable iff |g(w)| ≤ 1
this numerical scheme is unstable since |g(w)| > 1
修改遞迴公式，採用Lax–Friedrichs method，就有機會穩定。
f[i+1][j] = ½ (f[i][j+1] + f[i][j-1])
- ½ C (f[i][j+1] - f[i][j-1])
g(w) = ½ (exp(𝑖wΔx) + exp(-𝑖wΔx)) - ½ C (exp(𝑖wΔx) - exp(-𝑖wΔx))
g(w) = cos(wΔx) - C 𝑖 sin(wΔx)
let |C| ≤ 1 is CFL condition
thus |g(w)|² = cos²(wΔx) + C² sin²(wΔx)
≤ cos²(wΔx) + sin²(wΔx)
= 1
a numerical scheme is stable iff |g(w)| ≤ 1
this numerical scheme is stable under CFL condition |C| ≤ 1
Courant–Friedrichs–Lewy condition
物理學之連續方程式，一個線性方案，一旦滿足指數穩定性，也會連帶滿足某些特殊性質。此處要談的特殊性質是「CFL條件」。
for linear numerical scheme of continuity equation,
contractive stability => CFL condition
連續方程式請見後面章節。連續方程式的某些衍生版本（例如波動方程式、熱傳導方程式）符合這個結論，但是無法保證各種衍生版本都符合這個結論。這個結論不是推理而得，而是觀察而得。因此這個結論不被認為是數學定理。
CFL條件的數學式子：波動方程式是v Δt / Δx ≤ 1。熱傳導方程式是2 v ∆t / (∆x)² ≤ 1。其中v是速度。
證明請見講義：
【待補範例】
CFL條件的幾何意義：移動距離v Δt不得超過取樣距離Δx，否則空間離散化會抓錯樣本。
空間離散化必須抓到相鄰樣本。當移動距離太大，則空間離散化誤抓到從遠處跨越鄰居而來的樣本。
CFL 條件是必要條件：滿足收縮穩定性，導致滿足CFL條件。CFL條件不是充分條件：滿足CFL條件，可能穩定、可能不穩定。
for linear numerical scheme of continuity equation,
contractive stability ⇍ CFL condition
CFL條件的主要功用是快篩試劑。逆否命題：不滿足CFL條件，導致不滿足收縮穩定性。
https://scicomp.stackexchange.com/questions/2927/
https://scicomp.stackexchange.com/questions/25398
https://people.maths.ox.ac.uk/trefethen/pdetext.html
A-stability
專著《Numerical Methods for Ordinary Differential Equations》 。
先前只有考慮遞迴公式暨數值解的穩定性，現在還要考慮微分方程式暨符號解的穩定性。
微分方程式暨符號解的數學性質、遞迴公式暨數值解的數學性質，兩邊必須一致，包含穩定性。
B-stability：B穩定性。B是作者Butcher。
一旦微分方程式暨符號解滿足單調遞減性（等價於收縮穩定性），
那麼遞迴公式暨數值解必須滿足收縮穩定性。
A-stability：A穩定性。A是角度angle。
一旦微分方程式暨符號解滿足穩定性，
那麼遞迴公式暨數值解必須滿足絕對穩定性。
L-stability：L穩定性。L是左邊left。
一旦微分方程式暨符號解滿足穩定性，
那麼遞迴公式暨數值解必須滿足剛硬穩定性。
甚至深入觀察經典的數值方法的穩定性。
AN-stability：滿足B穩定性的情況下，
multistage method的數學性質。
G-stability：滿足A穩定性的情況下，
multistep method的數學性質。
algebraic stability：滿足A穩定性的情況下，
multistage method的數學性質。
https://encyclopediaofmath.org/wiki/Non-linear_stability_of_numerical_methods
https://www.jstor.org/stable/2156562
https://www.jstor.org/stable/2156895
http://cds.cern.ch/record/1069162
function of absolute stability
舉例來說，下述微分方程式：
d
—— f(t) = k f(t) k is a coefficient
dt
符號解：exp(k t)收斂至零或受限，(k t)是負數或零，才會滿足穩定性。
f(t) = f(0) exp(k t)
數值解：特定的數值方案，得到特定的ϕ(z)。
f[i+1] = ϕ(z) f[i] let z = k Δt
數學家利用ϕ(z)定義穩定性。此例當中：
absolute stability: |ϕ(z)| ≤ 1
A-stability：|ϕ(z)| ≤ 1 for all Re[z] ≤ 0
A-stability：{z : Re[z] ≤ 0} ⊆ {z : |ϕ(z)| ≤ 1}
L-stability：A-stability and ϕ(z)→0 as z→∞
針對A穩定性，t和Δt都是非負數，Re[k t] ≤ 0和Re[k] ≤ 0和Re[k Δt] ≤ 0意義相同。針對L穩定性，Δt → ∞和(k Δt) → ∞意義相同。為了讓數學式子變漂亮，數學家一律採用z = k Δt。
stability analysis with characteristic equation
關鍵字linear difference equation。
遞迴公式：
f[i] = c₁ f[i-1] + ... + cₘ f[i-m]
特徵方程式：
λᵐ = c₁ λᵐ⁻¹ + ... + cₘ λᵐ⁻ᵐ
移項整理，左式是特徵多項式：
λᵐ - c₁ λᵐ⁻¹ - ... - cₘ λᵐ⁻ᵐ = 0
特徵多項式分解：
(x - λ₁)...(x - λₘ) = 0
遞迴公式改寫成一般公式：
f[i] = c₁ λ₁ⁱ + ... + cₘ λₘⁱ
想要檢查指數穩定性，只需將遞迴公式改寫成特徵方程式，然後檢查所有根絕對值皆小於等於1。
characteristic polynomial of absolute stability
針對multistep method。數學家將遞迴公式改寫成特徵方程式，叫做characteristic polynomial of absolute stability，簡稱stability polynomial。
舉例來說，下述微分方程式：
微分方程式 d/dt f(t) = k f(t)
符號解 x = f(0) exp(k t)
多步法遞迴公式 (a₀ f[i] + ... + aₘ f[i-m])
= k Δt (b₀ f[i] + ... + bₘ f[i-m])
where a₀ = 1
特徵方程式 a(λ) = k Δt b(λ)
特徵方程式 a(λ) - k Δt b(λ) = 0
特徵方程式 k Δt = a(λ) / b(λ)
特徵多項式 π(λ, z) = a(λ) - z b(λ) where z = k Δt
特徵多項式分解 π(λ, z) = (x - λ₁)...(x - λₖ)
特徵多項式的根 r(z) = {λ : π(λ, z) = 0} = {λ₁, ..., λₖ}
數值解穩定性 |r(z)| ≤ 1
符號解穩定性 Re[z] ≤ 0
數學家利用r(z)定義穩定性。此例當中：
absolute stability：|r(z)| ≤ 1
A-stability：|r(z)| ≤ 1 for all Re[z] ≤ 0
L-stability：A-stability and r(z) → 0 as z → ∞
absolute stability：特徵多項式的根全部小於等於一。
A-stability：特徵多項式的根完全位於複數平面左半部。
L-stability：承上，當z趨近負無限大，則特徵多項式的根趨近零。
https://users.soe.ucsc.edu/~hongwang/AM213B/Notes/Lecture06.pdf
https://www.youtube.com/watch?v=fdbZpKyhpU0
https://www.math.unipd.it/~alvise/AN_2017/LETTURE/DAHLQUIST_STAB.pdf
https://link.springer.com/article/10.1007/BF01934126
此穩定性分析衍生許多定理。有興趣的讀者請自行研究。
Dahlquist's stability theorem
Dahlquist's order barrier theorem
Dahlquist's equivalence theorem
專論《Dahlquist's classical papers on stability theory》。
stability analysis
穩定性分析。使用數學工具推導穩定性。業已介紹，茲此總結。
穩定性有三種層次。總之都是收縮。
relative stability：考慮遞迴公式暨數值解的穩定性，當Δt與Δx趨近零。
absolute stability：考慮遞迴公式暨數值解的穩定性，針對特定大小的Δt與Δx。
A-stability：同時考慮微分方程式暨符號解的穩定性。
https://sites.math.rutgers.edu/~falk/math573/lecture21.pdf
http://minitorn.tlu.ee/~jaagup/uk/dynsys/ds2/num/Absolute/Absolute.html
https://www.google.com.tw/books?id=wbCjEAAAQBAJ&pg=PA404
穩定性分析有三種手法。總之都是頻譜。
(1) stability analysis with matrix
線性遞迴公式，改寫成矩陣，特徵值絕對值皆小於等於1。
(2) stability analysis with Fourier transform
(von Neumann stability analysis)
遞迴公式，改寫成弦波，振幅絕對值皆小於等於1。
(3) stability analysis with characteristic equation
(Dahlquist stability analysis)
遞迴公式，改寫成特徵方程式，根絕對值皆小於等於1。
特徵值／振幅／根的數學性質。總之都是單位圓裡面。
weak stability：弱穩定性。特徵值／振幅／根絕對值等於1至少兩個（重根）。
zero-stability：零穩定性。特徵值／振幅／根絕對值等於1至多一個。
strong stability：強穩定性。特徵值／振幅／根絕對值等於1恰好一個（簡單根）。
https://sites.math.rutgers.edu/~falk/math573/lecture21.pdf
fundamental theorem of numerical analysis特殊版本
數值分析基本定理的特殊版本，針對相對穩定性：
(1) Lax–Richtmyer equivalence theorem
IVP初始值問題、線性微分方程式、平滑解
線性方案、L₂遞減穩定性、L₂收斂性
(2) Dahlquist's equivalence theorem
IVP初始值問題、非線性微分方程式、平滑解
multistep method、零穩定性、L₂收斂性
(3) Lax–Wendroff consistency theorem & [DiPerna 1983]
IVP初始值問題、物理學之連續方程式／數學之守恆律、弱解
守恆型、L∞ 暨total variation受限穩定性、L₁收斂性
數值分析基本定理的特殊版本，針對絕對穩定性：
(1) Courant–Friedrichs–Lewy condition：
IVP初始值問題、物理學之連續方程式／數學之守恆律、平滑解
某些線性方案、CFL condition
(2) strong-stability preserving property：
IVP初始值問題、物理學之連續方程式／數學之守恆律、平滑解
Euler method家族、CFL condition
數值分析基本定理的符號版本。兩邊必須一致，包含收斂性。
(1) Cauchy—Lipschitz theorem
常微分方程ODE、初始值問題IVP，若是收縮函數，則存在唯一解。
(2) Cauchy—Kovalevskaya theorem
偏微分方程PDE、初始值問題IVP，若是收縮函數，則存在唯一解。
Lax–Richtmyer equivalence theorem參考資料：
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
線性遞迴公式，數值解，最大的奇異值不超過一
Lax–Wendroff consistency theorem參考資料：
描述 [Shu 1987] https://www.ams.org/journals/mcom/1987-49-179/S0025-5718-1987-0890256-5/S0025-5718-1987-0890256-5.pdf
證明 [LeVeque 1992] theorem 15.2
補充 https://math.stackexchange.com/questions/1529156/
numerical scheme
數值方案。面對各式各樣的微分方程式，想方設法調整遞迴公式，滿足一致性、滿足穩定性，導致收斂性。
具體來說有四個項目。目前沒有專有名詞。
(1) type of solution
(2) property of solution
(3) design of numerical method
(4) choice of numerical solver
一、觀察微分方程式／符號解的類型。
二、調整遞迴公式／數值解的數學性質。
三、設計遞迴公式。
四、解遞迴公式。
type of solution
微分方程式分為兩類。
(1) linear differential equation：一次。函數微分視作變數，形成一次函數。
(2) nonlinear differential equation：非一次。函數微分視作變數，形成非一次函數。
微分方程式的常見數學性質。
(a) quasilinear equation：擬一次。形如一次函數，惟係數不是常數函數。
(b) homogeneous equation：齊次。沒有常數項。導致符號解可以是零函數。
符號解分為四類。
(1) smooth solution：平滑解。無窮次可微函數。
(2) differentiable solution：可微解。可微函數。
(3) weak solution：弱解。也可以是不可微函數。
(4) multivalued solution：多值解。也可以不是函數。
線性微分方程式，初始條件是可微函數，符號解亦是可微函數。初始條件是不可微函數，符號解亦是不可微函數。
非線性微分方程式，初始條件是可微函數，符號解可能是可微函數、不可微函數、不是函數。
數值解沒有微分的概念。
https://math.stackexchange.com/questions/3649691/
smooth solution
「平滑解」。微分方程式的符號解，設定為平滑函數。
數學家從微積分基本定理出發，利用指數函數、三角函數，求得微分方程式的符號解。利用平滑函數，自然而然得到平滑解。
數學家從數值分析基本定理出發，利用微分運算的泰勒近似，求得微分方程式的數值解。數值解是近似解。全域誤差是近似誤差。想讓全域誤差趨近零，一種方式是採用平滑解。
differentiable solution
「可微解」。微分方程式的符號解，必須是可微函數。
詳細來說，內含N次微分的微分方程式的符號解，必須是N次可微函數。顯而易見的廢話。
理論上來說，可微解可以是平滑函數（處處無窮次可微），也可以是分段函數（銜接處N次可微）。然而目前無人討論分段函數。
weak solution
「弱解」。微分方程式的符號解，推廣為不可微函數。
原式乘上平滑函數、積分，得到新式。提升平滑程度。
針對各種平滑函數，新式一律出現的解，定義為原式的弱解。
balabala huluhulu
現實世界經常出現弱解。架空世界流行可微函數，現實世界流行不可微函數。藉由架空世界的可微的微分方程式，強行模擬現實世界的不可微的符號解。
知名範例是初始條件不是可微函數，導致解不是可微函數。
multivalued solution（set-valued solution）
「多值解」。非線性微分方程式的符號解，可能變得不是函數。
微分方程式的解，本來是函數；從某個時刻開始，變得不是函數，形成曲面，一個位置擁有多個函數值。
smooth approximation
不連續函數強行變成連續函數（弱解強行變成平滑解）。
數值解是不連續函數，導致Runge phenomenon。事先平滑化，成為可微函數。
初始值、邊界值是不連續函數。脈衝初始條件、對稱邊界條件。
conservative formulation
非函數強行變成函數（多值解強行變成單值解）。
垂直方向切上一刀，左右凹凸面積相等。
知名範例是交通車流問題。車輛不會互相穿越，整個過程都是函數。
實際範例是衝擊波shock wave。衝擊波不會互相交疊，整個過程都是函數。
property of solution
遞迴公式分為兩類。
(1) linear numerical scheme
遞迴公式是線性遞迴函數
（線性微分方程式＋線性差分）
(2) nonlinear numerical scheme
遞迴公式是非線性遞迴函數
遞迴公式暨數值解的常見數學性質。
(1) diminishing
遞減。數值解範數每回合變小或不變。
(2) local extremum diminishing
局部極值遞減。
遞迴公式的變數是自己和所有鄰居、係數是凸組合。
數值解極值數量不增加、極值大小不變廣。
(3) monotonicity preserving
保單調。數值解每回合維持單調。
(4) conservative
守恆。數值解每回合某種物理量總和維持定值。
微分方程式暨符號解的數學性質，遞迴公式暨數值解的數學性質，兩邊必須一致。簡單起見，此處只談後者。
diminishing（non-increasing）
「遞減」。數值解範數每回合變小或不變。
L∞ -diminishing max |f̂⁽ⁿ⁺¹⁾| ≤ max |f̂⁽ⁿ⁾|
L₁-diminishing ‖f̂⁽ⁿ⁺¹⁾‖₁ ≤ ‖f̂⁽ⁿ⁾‖₁
Lₚ-diminishing ‖f̂⁽ⁿ⁺¹⁾‖ₚ ≤ ‖f̂⁽ⁿ⁾‖ₚ
知名範例是微分不動點。【尚待檢查accuracy】
彈簧振動、琴弦振動、鼓膜振動的符號解，可以寫成複數的實部。儘管實部乍看反覆變大變小，但是複數其實保持變小或不變，滿足指數穩定性。【尚待確認】
d
—— f(t) = k f(t) k is a coefficient
dt
某些特殊情況下，收縮即是遞減。當微分方程式／遞迴公式沒有常數項，則符號解／數值解可以是零。此時收縮即是遞減（零作為減數）。
homogeneous & contractive => diminishing
diminishing的衍生性質。針對線性遞迴公式。
linear & contractive => diminishing => exponential decay
收縮即是遞減即是指數衰減。
local extremum diminishing（range diminishing）
專著《Computational Aerodynamics》
「局部極值遞減」。數值解同時滿足下述三點。
(1) nonincreasing local maximum 數值解局部極大值減少或不變
(2) nondecreasing local minimum 數值解局部極小值增加或不變
(3) no new extremum 數值解局部極值數量減少或不變
遞迴公式，變數是自己和所有鄰居，係數是凸組合。
一維空間是三點。二維空間是十字五點。
1D
{ fᵢ⁽ⁿ⁺¹⁾ = cᵢ₋₁ fᵢ₋₁⁽ⁿ⁾ + cᵢ fᵢ⁽ⁿ⁾ + cᵢ₊₁ fᵢ₊₁⁽ⁿ⁾
{ cᵢ₋₁ + cᵢ + cᵢ₊₁ = 1
{ cᵢ₋₁ , cᵢ , cᵢ₊₁ ≥ 0
2D
{ fᵢ,ⱼ⁽ⁿ⁺¹⁾ = cᵢ,ⱼ fᵢ,ⱼ⁽ⁿ⁾
{ + cᵢ₋₁,ⱼ fᵢ₋₁,ⱼ⁽ⁿ⁾ + cᵢ₊₁,ⱼ fᵢ₊₁,ⱼ⁽ⁿ⁾
{ + cᵢ,ⱼ₋₁ fᵢ,ⱼ₋₁⁽ⁿ⁾ + cᵢ,ⱼ₊₁ fᵢ,ⱼ₊₁⁽ⁿ⁾
{ cᵢ,ⱼ + cᵢ₋₁,ⱼ + cᵢ₊₁,ⱼ + cᵢ,ⱼ₋₁ + cᵢ,ⱼ₊₁ = 1
{ cᵢ,ⱼ , cᵢ₋₁,ⱼ , cᵢ₊₁,ⱼ , cᵢ,ⱼ₋₁ , cᵢ,ⱼ₊₁ ≥ 0
知名範例是Laplacian運算。【尚待檢查accuracy】
1D heat equation
∂ ∂²
—— f(t,x) = ——— f(t,x)
∂t ∂x²
forward Euler method
fᵢ⁽ⁿ⁺¹⁾ - fᵢ⁽ⁿ⁾ (fᵢ₋₁⁽ⁿ⁾ - fᵢ⁽ⁿ⁾) + (fᵢ₊₁⁽ⁿ⁾ - fᵢ⁽ⁿ⁾)
——————————————— = —————————————————————————————————————
Δt 2 Δx
forward Euler method
fᵢ⁽ⁿ⁺¹⁾ = (1 - C) fᵢ⁽ⁿ⁾ + ½ C fᵢ₋₁⁽ⁿ⁾ + ½ C fᵢ₊₁⁽ⁿ⁾
where C = Δt/Δx
LED的衍生性質。針對線性遞迴公式、擬線性遞迴公式。
local extremum diminishing (LED)
遞迴公式的變數是自己和所有鄰居、係數是凸組合
一、數值解局部極大值每回合減少或不變
二、數值解局部極小值每回合增加或不變
三、數值解局部極值數量每回合變少或不變
{ fᵢ⁽ⁿ⁺¹⁾ = sumⱼ cᵢⱼ fⱼ⁽ⁿ⁾
{ sumⱼ cᵢⱼ = 1
{ cᵢⱼ ≥ 0
{ cᵢⱼ = 0 , if j ∉ neighbor(i)
=> convex combination
遞迴公式的係數是凸組合
{ fᵢ⁽ⁿ⁺¹⁾ = sumⱼ cᵢⱼ fⱼ⁽ⁿ⁾
{ sumⱼ cᵢⱼ = 1
{ cᵢⱼ ≥ 0
=> positive combination (positivity)
遞迴公式的係數皆是正值
{ fᵢ⁽ⁿ⁺¹⁾ = sumⱼ cᵢⱼ fⱼ⁽ⁿ⁾
{ cᵢⱼ ≥ 0
=> order preserving (monotone)
遞迴公式是單調遞增函數（任一變數皆如此）
兩個數值解每回合維持高低關係
{ fᵢ⁽ⁿ⁺¹⁾ = Fᵢ(fⱼ₀⁽ⁿ⁾, fⱼ₁⁽ⁿ⁾, ..., fⱼₖ⁽ⁿ⁾)
{ fⱼₗ⁽ⁿ⁾ ≤ gⱼₗ⁽ⁿ⁾ => fᵢ⁽ⁿ⁺¹⁾ ≤ gᵢ⁽ⁿ⁺¹⁾ for all 0 ≤ l ≤ k
f̂⁽ⁿ⁾ ≤ ĝ⁽ⁿ⁾ => f̂⁽ⁿ⁺¹⁾ ≤ ĝ⁽ⁿ⁺¹⁾
[Harten 1983] [Osher 1985] [Jameson 1995]
LED能夠抑制振盪、漸趨平緩。此性質稱作maximum principle，源自convex combination。
振盪程度可以具體表示成數值大小L∞ -norm與total variation。此性質稱作L∞ -diminishing與total variation diminishing。
local extremum diminishing (LED)
=> convex combination
=> maximum principle
數值解下回合數值介於遞迴公式的變數極值之間
min {fⱼ⁽ⁿ⁾ | cᵢⱼ > 0} ≤ fᵢ⁽ⁿ⁺¹⁾ ≤ max {fⱼ⁽ⁿ⁾ | cᵢⱼ > 0}
=> L∞ -diminishing
數值解每回合振幅減少或不變
‖f̂⁽ⁿ⁾‖∞ ≤ ‖f̂⁽ⁿ⁺¹⁾‖∞
where ‖f‖∞ = maxᵢ |fᵢ|
local extremum diminishing (LED)
=> convex combination & neighborhood
=> maximum principle & neighborhood
數值解下回合數值介於本回合自己與鄰居的極值之間
（當邊界條件是指定數值，數值解最終回合介於邊界極值之間。）
min {fⱼ⁽ⁿ⁾ | j ∈ neighbor(i)} ≤ fᵢ⁽ⁿ⁺¹⁾ ≤ max {fⱼ⁽ⁿ⁾ | j ∈ neighbor(i)}
=> total variation diminishing (TVD)
數值解每回合總差值減少或不變
‖f̂⁽ⁿ⁾‖ᴛᴠ ≤ ‖f̂⁽ⁿ⁺¹⁾‖ᴛᴠ
where ‖f‖ᴛᴠ = sumᵢ |fᵢ₊₁ - fᵢ|
monotonicity preserving
「保單調」。數值解每回合維持單調。不會新增局部極值。
遞迴公式，自己和左右鄰居的中位數。然而數值解整體保持不變，僅極值被削平。
fᵢ⁽ⁿ⁺¹⁾ = median(fᵢ₋₁⁽ⁿ⁾, fᵢ⁽ⁿ⁾, fᵢ₊₁⁽ⁿ⁾)
我有找到一種做法是中位數。但是我沒有學會。
median(x,y,0) = minmod(x,y) = ½ (sgn(x) + sgn(y)) min(x,y)
median(x,y,z) = x + minmod(y-x, z-x)
vUL = vj + alpha * (vj - vj-1) vUL is upper bound of vj+1
where alpha ≥ 2 is slope 實務上≥4以避振
vUL - vj-1 ≤ (1 + alpha)(vj - vj-1)
vj+1/2 belongs [vj,vj+1]
w = vj - lambda * (Vj+l/2 - vj-l/2) Euler method
lambda ≤ 1/(1+alpha) 如此才能讓 w belongs [vj-1,vj]
v_mp = vj + minmod(vj+1 - vj, alpha*(vj - vj-1))
vj+1/2 = median(vj+1/2, vj, v_mp) monotoncity preserving scheme
[Suresh–Huynh 1997]
早期的數學家，將保單調遞迴公式硬是用於非單調數值解。數值解拆解成單調段落與局部極值段落，單調段落可以保持單調，局部極值段落則可能發生任何事情，例如振盪。調整係數以便緩解振盪。
保單調的衍生性質。針對線性遞迴公式。【尚待確認】
for linear numerical scheme,
monotonicity preserving => order preserving (monotone)
for linear numerical scheme with 3 nearest points,
monotonicity preserving => local extremum diminishing (LED)
One-dimensional Quasilinear Hyperbolic Systems of Conservation Laws and Their Applications
https://books.google.com.tw/books?id=AAqMDwAAQBAJ&pg=PA257
額外介紹一個我認為大家都搞錯的定理。
注意到，前者沒有設置邊界條件（數值解左右延伸至無限遠），後者卻設置邊界條件（數值解兩端固定），無法一路衍生。
L₁-contraction
=> total variation diminishing
=> monotonicity preserving ✘
[Harten 1983]
我認為正確的定理是這樣：
Dirichlet boundary condition & total variation diminishing (TVD)
數值解兩端固定＆數值解總差值縮小或不變（事實上只能保持不變）
{ f̂ᵢ⁽ⁿ⁾ = f̂ᵢ⁽ⁿ⁺¹⁾ for all i at boundary
{ ‖f̂⁽ⁿ⁾‖ᴛᴠ ≤ ‖f̂⁽ⁿ⁺¹⁾‖ᴛᴠ
=> monotonicity preserving
數值解每回合維持單調
(f̂⁽ⁿ⁾ is monotone => f̂⁽ⁿ⁺¹⁾ is monotone)
f is monotone:
min(fᵢ₋₁, fᵢ₊₁) ≤ fᵢ ≤ max(fᵢ₋₁, fᵢ₊₁) for all i
[Harten 1983]
conservative
「守恆」。數值解每回合某種物理量總和維持定值。
物理學已經建立許多守恆定律，例如質量守恆、動量守恆、能量守恆。各種微分方程式，擁有各種守恆數學式子。
最簡單的情況：數值解每一項總和是定值。
sumᵢ f̂ᵢ⁽ⁿ⁾ = sumᵢ f̂ᵢ⁽ⁿ⁺¹⁾
design of numerical method
遞迴公式的設計方式。
(1) finite difference method：有限差分法。求平滑解。使用原本的微分方程式。
(2) finite volume method：有限體積法。求弱解。改成通量（體積流速）。對空間積分、對時間平均。
(3) level set method：等高線法。求多值解。改成勢。
finite difference method
「有限差分法」。大家分開討論時間與空間。
temporal discretization：對時間變數偏微分（左式）的離散化
spatial discretization：對空間變數偏微分（右式）的離散化
temporal discretization：左式離散化。
(1) 取左方斜率、取右方斜率。
(2) multistep method：取多處斜率。
(3) multistage method：求多次斜率。
(4) multiderivative method：取多階微分。
多步法遞迴公式
f[i] = (a₁ f[i-1] + ... + aₘ f[i-m]) + k Δt (b₀ f[i] + ... + bₘ f[i-m])
Adams method: Maclaurin series
https://mathworld.wolfram.com/AdamsMethod.html
Richardson extrapolation
https://en.wikipedia.org/wiki/Richardson_extrapolation
圖表稱作Butcher tableau。數學家已經發明許多種計算公式，各有利弊。
spatial discretization：右式離散化。
(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/dt f(t) = f(t) + 1 ---> (f[i] - f[i-1]) / Δt = f[i] + 1
^^^^
explicit。取舊值。
d/dt f(t) = f(t) + 1 ---> (f[i] - f[i-1]) / Δt = f[i-1] + 1
^^^^^^
operator splitting。項次太多，逐步更新。
d/dt f(t) = f(t) + g(t) ---> 1. (f[i] - f[i-1]) / Δt = f[i]
2. (f[i] - f[i-1]) / Δt = g[i]
semi-implicit。式子太多，設定順序。
d/dt f(t) = f(t) + g(t) ---> 1. (f[i] - f[i-1]) / Δt = g[i]
2. (f[i] - f[i-1]) / Δt = f[i]
圖表稱作stencil。數學家已經發明許多種計算公式，各有利弊。為了應付非線性微分方程式，花招百出。
【尚待確認】
多步、多階段、多導數，主要有兩個功用。
一、暗中縮短Δt與Δx。一階差分重複n回合，其實就是n階差分、n倍Δt與Δx。反過來說，n階差分讓Δt與Δx變成1/n倍、一口氣做n次差分。
二、微調答案。取大量斜率的加權平均，讓數值解更符合符號解。直接來說，就是修改係數、通靈答案。缺乏科學根據。
取左方斜率，無論是一處或多處，無論是遞增或遞推，無論是顯式或隱式，只要強行展開數學式子，多階段就會變成單階段，缺乏討論意義。取右方斜率，才有討論意義。
重新歸納為左右兩種。宛如差分。
Euler method：取當地斜率。參考先前當前函數點。
multistep method：取左方斜率。參考先前幾個函數點。優點是快。
multistage method：取右方斜率。預測之後幾個函數點。優點是準。
Shu–Osher strong-stability preserving theorem
專著《Strong Stability Preserving Runge–Kutta and Multistep Time Discretizations》。
時間離散化可以推廣為多步多階段。
「保強穩定性定理」。如果一個數值方案，時間離散化採用forward Euler method，一旦滿足CFL condition，並且導致收縮穩定性，那麼該數值方案，時間離散化推廣為多步多階段顯式隱式，其係數是凸組合，一旦滿足CFL condition乘上適當倍率，也會導致收縮穩定性。
實務上會仔細調整時間離散化的係數，使得適當倍率是1，如同普通的CFL condition。
實務上將時間離散化推廣為Runge–Kutta method，形成strong-stability preserving Runge–Kutta method (SSPRK)。
「保強穩定性定理」原理簡單。foward Euler method只取一處斜率。多步多階段則是取多處斜率的加權平均數。視作多個foward Euler method各自加權，而多個CFL condition也跟著加權。多個CFL condition不等式聯立，簡化成最緊的那一個不等式。
此處僅證明一步多階段顯式線性遞迴公式，可以推廣到多步多階段隱式一般遞迴公式。
differential equation
d/dt f(t,x) = F(f(t,x))
forward euler method
fⁿ⁺¹ = fⁿ + Δt F(fⁿ) skip subscript
contractive stability of linear scheme
‖fⁿ⁺¹‖ = ‖fⁿ + Δt F(fⁿ)‖ ≤ ‖fⁿ‖ for simplicity
explicit multistage method
⎧ f⁽⁰⁾ = fⁿ skip subscript
⎨ f⁽ⁱ⁾ = sum { aᵢⱼ f⁽ʲ⁾ + Δt bᵢⱼ F(f⁽ʲ⁾) } 1 step
⎪ j=1⋯i-1 already expanded
⎩ fⁿ⁺¹ = f⁽ᵐ⁾ m stages
convex combination property of each stages
⎧ aᵢⱼ ≥ 0 for all i and j
⎨ bᵢⱼ ≥ 0 for all i and j
⎪ sum aᵢⱼ = 1 for all i
⎩ ʲ
strong-stability preserving property:
if forward euler method reaches contractive stability
under CFL condition Δt ≤ ΔtCFL ,
then explicit multistage method reaches contractive stability
under CFL condition Δt ≤ min (aᵢⱼ/bᵢⱼ) ΔtCFL
ⁱʲ
when a and b are non-negative
and a satisfies convex combination property.
proof:
‖f⁽ⁱ⁾‖
= ‖sum { aᵢⱼ f⁽ʲ⁾ + Δt bᵢⱼ F(f⁽ʲ⁾) }‖ ith stage
≤ ‖aᵢⱼ‖ ‖sum { f⁽ʲ⁾ + Δt (bᵢⱼ/aᵢⱼ) F(f⁽ʲ⁾) }‖ 三角不等式
= aᵢⱼ ‖sum { f⁽ʲ⁾ + Δt (bᵢⱼ/aᵢⱼ) F(f⁽ʲ⁾) }‖ aᵢⱼ ≥ 0
≤ aᵢⱼ sum ‖ f⁽ʲ⁾ + Δt (bᵢⱼ/aᵢⱼ) F(f⁽ʲ⁾) ‖ 正負抵銷
≤ aᵢⱼ sum ‖f⁽ʲ⁾‖ 原因如後
≤ max ‖f⁽ʲ⁾‖ 凸組合
⎧ ‖f⁽¹⁾‖ ≤ max(‖f⁽⁰⁾‖) = ‖f⁽⁰⁾‖
⎨ ‖f⁽²⁾‖ ≤ max(‖f⁽⁰⁾‖, ‖f⁽¹⁾‖) ≤ ‖f⁽⁰⁾‖
⎪ : :
⎩ ‖f⁽ᵐ⁾‖ ≤ max(‖f⁽⁰⁾‖, ‖f⁽¹⁾‖, ..., ‖f⁽ᵐ⁻¹⁾‖) ≤ ‖f⁽⁰⁾‖
‖fⁿ⁺¹‖ = ‖f⁽ᵐ⁾‖ ≤ ‖f⁽⁰⁾‖ = ‖fⁿ‖
for each j,
forward Euler method reaches contractive stability
‖f⁽ʲ⁾ + Δt (bᵢⱼ/aᵢⱼ) F(f⁽ʲ⁾)‖ ≤ ‖f⁽ʲ⁾‖
under CFL condition Δt (bᵢⱼ/aᵢⱼ) ≤ ΔtCFL
「保強穩定性定理」名稱很雲。此概念既非強穩定性、亦非保Ｘ變換。此概念是將時間離散化從一階泰勒展開推廣成多步多階段顯式隱式，並且考慮CFL condition。
「保強穩定性定理」內容很雲。CFL condition只是必要條件，此定理卻以充分條件作為前提。前提是一件未必成立的事情。
因此，實務上需要確認前提是否成立。穩定性分析，確認收縮穩定性等價於CFL condition。一旦前提成立，SSPRK所向無敵。
Godunov's order barrier theorem
空間離散化的截斷誤差受到侷限。
「階數屏障定理」。具備特殊性質的遞迴公式暨數值解，局部誤差的Δx的次方值有上限。誤差不可能更小了、精度不可能更高了。
階數屏障定理有各種版本，這裡記載其中三種。Godunov首開先例。
for linear numerical scheme
of one-dimensional scalar differential function,
(1) LED scheme is at most second-order accurate.
(since stencil has 3 nearest points)
[?]
(2) monotoncity preserving scheme is at most first-order accurate.
(since stencil has 2 nearest points)
[Godunov 1954]
(3) multistep method is at most second order accurate.
[Dahlquist 1978]
https://hal.science/hal-01620642
https://en.wikipedia.org/wiki/Godunov's_theorem
https://www.jstor.org/stable/2008046
Dahlquist's stability theorem
時間離散化的穩定性受到侷限。
「穩定性定理」。特殊的遞迴公式暨數值解，可能不穩定。
Dahlquist的書沒有給出證明。嗚嗚。
(1) implicit multistep method is A-stable
under balabala....
(2) explicit multistep method is not A-stable.
一、隱式左側遞迴公式，指數穩定性變成：特徵方程式分式是負數或零（實部為負或零）
二、顯式左側遞迴公式，無論如何不可能穩定。
https://www.math.unipd.it/~alvise/AN_2017/LETTURE/DAHLQUIST_STAB.pdf
one-dimensional scalar hyperbolic conservation law
專著《Computational Gasdynamics》。
各種連續方程式擁有各種守恆解。數學家正在探索當中。此處只談最簡單的情況：物理學之連續方程式／數學之雙曲線偏微分方程式。一維、純量。
shock fitting method：衝擊波斜切，凹凸面積相等，變成函數。
shock capturing method：衝擊波縱切，凹凸面積相等，滿足jump condition。
one-dimensional scalar hyperbolic conservation law:
有限差分法。
此處只談最簡單的情況：物理學之連續方程式／數學之雙曲線偏微分方程式。一維（單變數）、純量（單變量）。
Lax–Friedrichs method: 1st order
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (f(uᵢ₊₁)⁽ⁿ⁾ - f(uᵢ₋₁)⁽ⁿ⁾) / 2
uᵢ⁽ⁿ⁺¹⁾ = (uᵢ₊₁⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾) / 2
- C (f(uᵢ₊₁)⁽ⁿ⁾ - f(uᵢ₋₁)⁽ⁿ⁾) / 2
Lax–Wendroff method: 2nd order
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾
- C (f(uᵢ₊₁)⁽ⁿ⁾ - f(uᵢ₋₁)⁽ⁿ⁾) / 2
+ C² f′(uᵢ₊₁⸝₂ ; f(uᵢ₊₁⁽ⁿ⁾) - f(uᵢ⁽ⁿ⁾)) / 2
- C² f′(uᵢ₋₁⸝₂ ; f(uᵢ⁽ⁿ⁾) - f(uᵢ₋₁⁽ⁿ⁾)) / 2
flux averaging 凸加權平均 一階二階
flux limited fⁿ = f⁽¹⁾ + w (f⁽²⁾ - f⁽¹⁾)
flux corrected fⁿ = f⁽¹⁾ + diff(f⁽²⁾, f⁽¹⁾)
flux splitting: increasing + decreasing
flux splitting uⁿ = u⁽¹⁾ + C (f⁺ - f⁻)
speed splitting
upwind method: conditional function，確保形成凸組合。
uᵢ⁽ⁿ⁺¹⁾ = { uᵢ⁽ⁿ⁾ - C a (uᵢ⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾) , if a > 0
{ uᵢ⁽ⁿ⁾ - C a (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) , if a < 0
upwind method
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C a⁺ (uᵢ⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾)
- C a⁻ (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾)
where a⁺ = max(a,0)
a⁻ = min(a,0)
|a| = a⁺ - a⁻
one-dimensional scalar hyperbolic conservation law:
conservation law：此處只談最簡單的情況：物理學之連續方程式／數學之雙曲線偏微分方程式。一維、純量。
d/dt u = d/dx f(u)
「有限體積法」。
時間變化對空間積分＝空間變化對時間積分。
weak solution of conservation law：「守恆律的弱解」。方便起見，以下簡稱為「守恆解」。針對連續方程式，滿足物理量守恆的解。
守恆解是不連續函數。守恆解可以視作弱解。數學家一律使用弱解來解釋守恆解。數學家不採用守恆解這個名詞。
conservation form：守恆型。前綴和的相鄰差，必須等於原本數值。flux滿足一致性。
前綴和可以調整大小。衍生各式各樣的flux limiter。
順帶一提，例如最短路徑演算法Johnson's algorithm的調整路徑手法，把邊權重調整為無負環。也許這種手法也可以把流場調整為無旋場。然而我們對調整函數知之甚少。等你發表論文。
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (f(uᵢ₊₁⸝₂)⁽ⁿ⁾ - f(uᵢ₋₁⸝₂)⁽ⁿ⁾)
where C = Δt / Δx
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)
discretization method
(1) vertex-centered finite volume scheme
(2) cell-centered finite volume scheme
numerical method：前面小節的數值方案，通通可以改寫成守恆型，而且通通收斂至守恆解。然而，通通是單調函數，通通是一階精準。因此數學家又發明了高精度的數值方案，。
semi-Lagrangian method
(1) reconstruction：一個格子，重新設定水面形狀。分段內插（加權平均）。
(2) advection：根據各格速度，挪動各格水面。不可超過一格，即是CFL condition。
(3) remapping：重新統計每個格子獲得多少水（加權平均）。
Godunov method: upwind method推廣成守恆型與連續函數
(1) reconstruction: Godunov flux
(2) advection + remapping = recurrence
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)
fᵢ₊₁⸝₂⁽ⁿ⁾ = Godunov_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)
= { min { f(u*) | uᵢ⁽ⁿ⁾ ≤ u* ≤ uᵢ₊₁⁽ⁿ⁾ }, if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
{ max { f(u*) | uᵢ₊₁⁽ⁿ⁾ ≤ u* ≤ uᵢ⁽ⁿ⁾ }, if uᵢ⁽ⁿ⁾ > uᵢ₊₁⁽ⁿ⁾
for linear advection equation,
Godunov method = upwind method
monotonic upstream-centered scheme for conservation laws (MUSCL)
(1) reconstruction
(1-1) flux limiter curve between f(uᵢ₋₁⸝₂) and f(uᵢ₊₁⸝₂)
(1-2) Riemann solver fᵢ₊₁⸝₂⁽ⁿ⁾
(2) advection + remapping = recurrence
(2-2) temporal discretization
essentially non-oscillatory (ENO)
‖f̂⁽ⁿ⁾‖ᴛᴠ ≤ ‖f̂⁽⁰⁾‖ᴛᴠ + O((Δx)𐞥)
數值解總差值減少量受限
[Harten–Engquist–Osher–Chakravarthy 1987]
https://www.researchgate.net/publication/223605523
TVD-MUSCL flux splitting，得以套用TVD。
WENO 內插改成特殊權重的ENO。
flux reconstruction：一個格子，重新設定水面形狀。計算公式稱作通量限制器flux limiter，數量眾多，各有利弊。
constant flux limiter：水平面（常數函數）
van Leer flux limiter：斜平面（一次函數）
加權平均，數值是左鄰差、右鄰差，
權重是左鄰差絕對值、右鄰差絕對值。
https://amrvac.org/md_doc_limiter.html
https://en.wikipedia.org/wiki/Flux_limiter
local Riemann problem：局部黎曼問題。兩個相鄰格子，給定兩個高度、兩個速度，求出分界處通量。計算公式稱作黎曼求解器Riemann solver，數量眾多，各有利弊。
Roe：
HLLE：
HLLC：
https://princetonuniversity.github.io/Athena-Cversion/AthenaDocsUGRiemann.html
http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf
one-dimensional scalar hyperbolic conservation law:
專著《Computational Gasdynamics》。
1D linear advection equation
df df
-- + a -- = 0
dt dx
initial condition f(x,0) = f₀(x)
solution f(x,t) = f₀(x - at)
waveform f(x,t) at certain t
wave speed a
rightward a > 0
leftward a < 0
static a = 0
wavefront x - at = constant
characterics dx/dt = a
expansion a f(x,t) ≤ a f(y,t)
compression a f(x,t) ≥ a f(y,t)
for bʟ(t) ≤ x ≤ y ≤ bʀ(t)
one-dimensional scalar hyperbolic PDE
df dφ(f)
-- + ----- = 0
dt dx
chain rule
df dφ df df df
-- + -- -- = 0 => -- + φ′(f) -- = 0
dt df dx dt dx
initial condition f(x,0) = f₀(x)
solution f(x,t) = f₀(x - φ′(f(x,t)) t)
waveform f(x,t) at certain t
wave speed φ′(f(x,t))
rightward φ′(f(x,t)) > 0
leftward φ′(f(x,t)) < 0
static φ′(f) = 0
wavefront x - φ′(f)t = constant
characterics dx/dt = φ′(f)
expansion φ(f(x,t)) ≤ φ(f(y,t))
compression φ(f(x,t)) ≥ φ(f(y,t))
for bʟ(t) ≤ x ≤ y ≤ bʀ(t)
one-dimensional scalar hyperbolic PDE
d/dt f = - d/dx φ(f) where φ(x) is known
solution
f(t,x)
discontinuous position of solution
p(t)
move speed of discontinuous position
s(t) = p′(t)
function value of discontinuous position of solution
f⁻(t) = lim f(t,x)
x→p(t)⁻
f⁺(t) = lim f(t,x)
x→p(t)⁺
compression
φ(f⁻) ≥ φ(f⁺) at t
Rankine–Hugoniot jump condition
左側溢出等於右側溢出。物理量守恆。
φ(f⁻) - s f⁻ = φ(f⁺) - s f⁺ for all t
Rankine–Hugoniot jump condition
改寫數學式子，宛如平緩函數。斷點函數值變化，符合斷點位置變化。
φ(f⁻) - φ(f⁺)
------------- = s for all t
f⁻ - f⁺
Oleinik admissibility condition
左側斷點收縮率總是大於（小於）右側斷點收縮率。
φ(f⁻) - φ(f) φ(f) - φ(f⁺) for all t
------------ ≥ s ≥ ------------ for all f
f⁻ - f f - f⁺ that f⁻ > f > f⁺
Lax admissibility condition
f恰是凸函數的情況。
φ′(f⁻) > s > φ′(f⁺) for all t
Riemann initial condition
方便教學的簡易範例。不連續處只有一個，函數值是兩個常數。
f(0,x) = { f⁻(0) if x ≤ p(0)
{ f⁺(0) if x > p(0)
one-dimensional scalar hyperbolic PDE
weak solution existence <=> Rankine–Hugoniot jump condition
(conservative solution)
weak solution uniqueness <=> Rankine–Hugoniot jump condition
(shock wave) & Oleinik admissibility condition
weak solution uniqueness <=> Rankine–Hugoniot jump condition
(shock wave) & Lax admissibility condition
when φ is convex
Riemann problem: Rankine–Hugoniot jump condition
& Lax admissibility condition
& Riemann initial condition
(1) 稀疏波rarefaction wave：初始條件左速小於右速φ(f⁻) < φ(f⁺)
(2) 衝擊波shock wave ：初始條件左速大於右速φ(f⁻) > φ(f⁺)
https://web.stanford.edu/class/math220a/handouts/conservation.pdf
https://www.researchgate.net/publication/230676056
entropy pair (η,q)
η(u) is entropy. η″(u) > 0 (strictly convex)
q(u) is entropy flux. q′(u) = f′(u) η′(u)
entropy inequality
d/dt η(u) + d/dx q(u) ≤ 0 for all entropy pairs
所有凸函數皆使得不等式成立。
admissible solution (entropy solution)
d/dt η(u)ϕ + d/dx q(u)ϕ ≥ 0 for all test function ϕ ≥ 0
弱解的其中一個特例：很純的衝擊波
Kruzhkov entropy inequality = Oleinik admissibility condition
{ d/dt η(u) + d/dx q(u) ≤ 0
{ η(u, k) = |u - k|
{ q(u, k) = (f(u) - f(k)) sgn(u - k)
其中一種凸函數。導致L₁-contraction。
（用來證明可容解是唯一解。）
Lₚ entropy
η(u) = |u|ᵖ / p
其中一種凸函數。導致Lₚ-contraction。
admissible solution => Lₚ-diminishing & total variation diminishing (TVD)
可容解具備Lₚ-diminishing和TVD。
畢竟守恆律沒有常數項。
[DiPerna 1983]
一、當任何一個凸函數使得不等式成立，而且總變差受限，就保證符號解L∞ norm受限。
二、當符號解L∞ norm受限、總變差受限，就保證符號解L₁ norm收斂。（穩定即收斂）
single entropy inequality & bounded variation
{ d/dt η(u) + d/dx q(u) ≤ 0 for any (η,q)
{ ‖u⁽ⁿ⁾‖ᴛᴠ ≤ B when n≥0
=> bounded L∞ norm & bounded variation
符號解每回合振幅／總差值絕對值受限
{ ‖u⁽ⁿ⁾‖∞ ≤ M when n≥0
{ ‖u⁽ⁿ⁾‖ᴛᴠ ≤ B when n≥0
=> L₁-convergence
符號解絕對值總和收斂至零
‖u⁽ⁿ⁾‖₁ → 0 when n→∞
[DiPerna 1983]
https://ntrs.nasa.gov/api/citations/19840010891/downloads/19840010891.pdf
https://ocw.mit.edu/courses/16-920j-numerical-methods-for-partial-differential-equations-sma-5212-spring-2003/b83f5998bf019571a250acd9acb7aa46_lec11_notes.pdf
https://math.stackexchange.com/questions/4083618/
https://math.stackexchange.com/questions/18395/
https://math.stackexchange.com/questions/271658/
https://math.stackexchange.com/questions/782558/
https://math.stackexchange.com/questions/1529156/
bounded L1 => bounded L2
L2 convergence => L1 convergence
注意到上述定理無法直接推廣為多變數（高維）／多變量（向量）／多式聯立（方程組）、無法直接引入散度／旋度／梯散、無法推廣為不規則網格／特殊邊界條件。每種連續方程式都必須特地檢查符號解的存在性／唯一性／穩定性／收斂性，是數學界的大難題。
Convergence of approximate solutions of conservation laws
https://www.igpm.rwth-aachen.de/Download/reports/pdf/IGPM217.pdf
one-dimensional scalar hyperbolic conservation law:
專著《Numerical Methods for Conservation Laws
and Related Equation》。
consistency：三個條件。第三個條件只是純粹撞名。古人考慮不周，成為歷史共業。
flux averaging
1 ⌠ xᵢ₊₁⸝₂
f̂ᵢ = —— ⎮ f(tₙ,x) dx when Δx→0
Δx ⌡ xᵢ₋₁⸝₂
consistency φ̂(f̂ᵢ,f̂ᵢ,...,f̂ᵢ) = φ(f̂ᵢ)
conservation sumᵢ f̂ᵢ⁽ⁿ⁾ = sumᵢ f̂ᵢ⁽ⁿ⁺¹⁾
uniqueness：
conservation form
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (Fᵢ₊₁⸝₂⁽ⁿ⁾ - Fuᵢ₋₁⸝₂⁽ⁿ⁾)
where C = Δt / Δx
F(u,u) = φ(u)
Fᵢ₊₁⸝₂ = F(uᵢ,uᵢ₊₁)
discrete entropy inequality [Tadmor 1984]
η(uᵢ⁽ⁿ⁺¹⁾) ≤ η(uᵢ⁽ⁿ⁾) - C (Qᵢ₊₁⸝₂⁽ⁿ⁾ - Qᵢ₋₁⸝₂⁽ⁿ⁾)
where entropy pair (η,q)
Q(u,u) = q(u)
Qᵢ₊₁⸝₂ = Q(uᵢ,uᵢ₊₁)
數值方案一旦符合此式，就是可容解。
entropy satisfying scheme
目前有三種滿足discrete entropy inequality的方法
1. monotone scheme：僅滿足其中一種凸函數，達成L₁收斂性。至多一階精準。
2. Osher's E scheme：滿足所有凸函數，達成所有收斂性。至多一階精準。
3. 人工黏性項。
[Tadmor 1984]
upwind Lax–Friedrichs method符合此式，得到可容解。
[Osher 1985]
Godunov method符合此式，得到可容解。
incremental form(單調三點)符合此式，得到可容解。
MUSCL符合此式，得到可容解。
不過證明的是semi-discrete........
https://www.math.umd.edu/~tadmor/references/files/Tadmor%20entropy%20stable%20schemes-r%20handbook2016.pdf
https://ocw.mit.edu/courses/16-920j-numerical-methods-for-partial-differential-equations-sma-5212-spring-2003/bafff8bb84c0ee8ac362cea89c5fea38_lec12_notes.pdf
stability：單調遞增函數追加條件，得到穩定性。各種範數。
主角是單調遞增函數，原因是非線性。
conservation容易聯想到L₁範數。
consistent & conservative & monotone & 3pt
=> discrete entropy condition
=> Lₚ-diminishing (1 ≤ p < ∞)
consistent & monotone
=> maximum principle
=> L∞ -diminishing & total variation diminishing (TVD)
conservative & monotone
=> Kruzhkov entropy inequality
=> L₁-contraction
=> total variation diminishing (TVD)
[Harten–Hyman–Lax–Keyfitz 1976] [Crandall–Majda 1980]
一、f ≤ max(f,g) => F(f) ≤ F(max(f,g))
二、數值解相減，結果有正有負，拆開正負，分開累計。
[Harten 1983]
數值解相減，令被減數是減數往右移位，硬是湊成相鄰差。
https://www.researchgate.net/publication/285599438
https://www.researchgate.net/publication/230676056
E scheme => TVD
conservative TVD => weak solution
convergence：單調遞增函數配合守恆律，導致可容條件成立、L₁收縮穩定性成立，進一步導致數值分析基本定理成立。守恆型一致性、L₁收縮穩定性、L₁收斂性。
consistent & conservative & monotone => converge
[Lax 1971] 台大圖書館沒有，要通靈了
https://www.sciencedirect.com/science/article/abs/pii/B9780127758503500182
守恆型一致姓 & L₁收縮穩定性 => L₁收斂性
符號版本：兩邊必須一致。
專著《Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves》。
symbolic solution | numerical solution
----------------------------+----------------------------------
Lax admissibility condition | order preserving (monotone)
Lax entropy inequality | discrete entropy inequality
tame oscillation condition | total variation diminishing (TVD)
time continuity estimate | L₁-contraction
https://bpb-us-e1.wpmucdn.com/sites.psu.edu/dist/d/80666/files/2019/09/clawtut4-Oxford.pdf
level set method
我還沒有學會。關鍵字也許是capturing multivalued solution。Hamilton–Jacobi Equations。
Stanley Osher
https://link.springer.com/book/10.1007/b98879
James Sethian
https://math.berkeley.edu/~sethian/2006/Publications/Book/book.html
Hailiang Liu
https://faculty.sites.iastate.edu/hliu/
Ken Museth
https://ken.museth.org/Publications.html
Multi-Valued Solution and Level Set Methods in Computational High Frequency Wave Propagation
https://faculty.sites.iastate.edu/hliu/files/inline-files/Liu-osher-tsai-CICP06_1.pdf
Front Tracking for Hyperbolic Conservation Laws
https://www.uio.no/studier/emner/matnat/math/MAT4380/v15/pensumliste/chapter3.pdf
A Level Set Method for the Computation of Multivalued Solutions to Quasi-Linear Hyperbolic PDEs and Hamilton-Jacobi Equations
choice of numerical solver
遞迴公式的求解方式。
(1) finite element method：有限元素法。取樣擇鄰：隨機點／網格／基底內插。
（本站稱為函數資料結構mesh discretization。）
(2) meshfree method：無網格法。隨機點／半徑範圍／基底內插。
（本站稱為函數資料結構meshfree discretization。）
(3) spectral method：頻域法。空域轉頻域。
(4) multigrid method：多格法。矩陣求解採用scaling method。
https://en.wikipedia.org/wiki/Numerical_methods_for_partial_differential_equations
