functional equation - 演算法筆記

differential equation - dynamical system🚧

dynamical system（time-dependent differential equation）

本文介紹微分方程式的其中一種特例：動態系統。

一、未知函數擁有一個時間變數t、任意個空間變數xyz。

二、未知函數對時間微分，出現在方程式之中。

（簡單起見，下方數學式子省略係數、相乘、變換。）

         ∂           ∂
f(t,x) + —— f(t,x) + —— f(t,x)
         ∂t          ∂x

  ∂²           ∂²           ∂²
+ ——— f(t,x) + ——— f(t,x) + ———— f(t,x) + ...
  ∂t²          ∂x²          ∂t∂x

+ 1 + t + x + t² + tx + x² + ... = 0

t: time variable
x, y, z: space variables

動態系統的特色：解是動畫。解隨著時間改變。

注意到，動態系統的時間變數t和空間變數xyz，不一定非得是物理學的時間和空間，也可以是其他事物。取名時間變數和空間變數，是為了形成動畫。

重新整理方程式，以便形成動畫。

未知函數只對時間微分，一律挪到左式。

（簡單起見，左式恰好只有一項：未知函數對時間一次微分。）

∂                      ∂           ∂²           ∂²
—— f(t,x) = - f(t,x) - —— f(t,x) - ——— f(t,x) - ———— f(t,x) - ...
∂t                     ∂x          ∂x²          ∂t∂x

            - 1 - t - x - t² - tx - x² - ...

左式是前後兩幀的變化程度。右式重新包裝成泛函數。形成動態系統的標準式。

∂
—— f(t,x) = F(f(t,x), t, x)
∂t

symbolic solution

動態系統的符號解，已經形成一門複雜的學問。由於這不是我的專長，就不多提了。

numeric solution

動態系統的數值解，我分成四章。第一章實務，寫程式求得數值解。第二章理論，證明數值解符合符號解。第三章是進階實務，探討數值解的處理技巧。第四章進階理論，探討數值解的數學性質。

discretization

未知函數實施離散化，以便計算數值解。

先前文件已經介紹，不再贅述。簡單起見，本文討論：方格、單點數值、有限差分、浮點數運算。

(1) domain discretization           定義域離散化
  (a) grid (mesh-based method)      方格（網格法）
  (b) mesh (mesh-based method)      網格（網格法）
  (c) particles (meshfree method)   粒子（無網格法）
(2) codomain discretization         對應域離散化
  (a) function value at point (FDM) 單點數值（有限差分法）
  (b) average function value (FVM)  區間平均（有限體積法）
  (c) shape function (FEM)          形狀函數（有限元素法）
  (d) sinusoid (spectral method)    弦波（頻譜法）
(3) operator discretization         運算離散化
  (a) numerical arithmetic          四則運算離散化
  (b) numerical differentiation     微分運算離散化（有限差分）
  (c) numerical integration         積分運算離散化
(4) number discretization           數字離散化
  (a) floating-point arithmetic     四則運算離散化（浮點數運算）

time-stepping method

大家以動態系統的標準式求得數值解。給定當前時刻的數值解，求得下個時刻的數值解。這個手法稱作time-integration method或time-stepping method。宛如gradient descent。

∂
—— f(t,x) = F(f(t,x), t, x)                 dynamical system
∂t

f(t+Δt,x) - f(t,x)
—————————————————— = F(f(t,x), t, x)        discretization
        Δt

f(t+Δt,x) = f(t,x) + F(f(t,x), t, x) ⋅ Δt   time-stepping
                                            equation
^^^^^^^^^   ^^^^^^   ^^^^^^^^^^^^^^^   ^^
next        current  viewed as slope   time
frame       frame                      step

void numerical_simulation()
{
	// time duration
	t₀ := beginning_time;
	tₙ := ending_time;

	// initialize first frame
	f(t₀, ⋅) := initial_values;

	// draw first frame
	draw(f(t₀, ⋅));

	// time-stepping loop
	for (t = t₀; t <= tₙ; t += Δt)
	{
		// compute next frame
		f(t+Δt, ⋅) := f(t, ⋅) + F(i, t, ⋅) * Δt;

		// draw next frame
		draw(f(t+Δt, ⋅));
	}
}

numerical simulation

數值模擬。微分方程式，求得數值解。四個步驟：

(1) numerical scheme：數值方案。精心設計的遞迴公式，使得數值收斂。
(2) numerical iteration：數值遞推。套用遞迴公式，逐步求得數值解。
(3) numerical error：數值誤差。數值遞推過程，數值解與符號解的差距。
(4) numerical convergence：數值收斂。數值遞推過程，差距趨近零。

numerical scheme

數值方案。精心設計的遞迴公式，使得數值收斂。

大家分開討論時間與空間。考慮三個項目：

(1) temporal discretization：時間離散化。時間步進的步伐大小與斜率大小。
(2) spatial discretization：空間離散化。取樣與擇鄰。
(3) implicit method：隱方法。右式是下個時刻。左式右式都有未知項。
    explicit method：顯方法。右式是當前時刻。僅左式有未知項。

temporal discretization

時間離散化。就是先前章節提到的三種方式：

Euler method：歐拉法。使用當前時刻，求得下個時刻。
multistep method：多步法。使用當前時刻與多個過往時刻，求得下個時刻。
multistage method：多階段法。細分當前時刻至下個時刻，反覆求得這些時刻。

實務上採用歐拉法。想要減少數值誤差則換成多階段法。花招百出，各有利弊。

舉例來說，下述微分方程式：

d
—— f(t,x) = F(f(t), t, x)
dt

時間離散化。

Euler method:

f(t+Δt,x) - f(t,x)
—————————————————— = F(f(t), t, x)
        Δt

multistep method (AB2):

⎧ k₁ = F(f(t,x), t, x)
⎪ k₂ = F(f(t-Δt,x), t-Δt, x)
⎨ 
⎪ f(t+Δt,x) - f(t,x)
⎪ —————————————————— = (3/2) k₁ - (1/2) k₂
⎩         Δt

multistage method (RK4):

⎧ k₁ = F(f(t,x), t, x)
⎪ k₂ = F(f(t,x) + k₁ ½Δt, t + ½Δt, x)
⎪ k₃ = F(f(t,x) + k₂ ½Δt, t + ½Δt, x)
⎨ k₄ = F(f(t,x) + k₃ 1Δt, t + 1Δt, x)
⎪
⎪ f(t+Δt,x) - f(t,x)
⎪ —————————————————— = (1/6)k₁ + (2/6)k₂ + (2/6)k₃ + (1/6)k₄
⎩         Δt

改寫成陣列。

Euler method:

f[n+1][i] = f[n][i] + F(f[n][i], t, x) ⋅ Δt

multistep method (AB2):

⎧ k₁ = F(f[n][i], t, x)
⎪ k₂ = F(f[n-1][i], t-Δt, x)
⎨ 
⎪ f[n+1][i] - f[n][i]
⎪ ——————————————————— = (3/2)k₁ - (1/2)k₂
⎩          Δt

multistage method (RK4):

⎧ k₁ = F(f[n][i], t, x)
⎪ k₂ = F(f[n][i] + k₁⋅(1/2)Δt, t+(1/2)Δt, x)
⎪ k₃ = F(f[n][i] + k₂⋅(1/2)Δt, t+(1/2)Δt, x)
⎨ k₄ = F(f[n][i] + k₃⋅Δt, t+Δt, x)
⎪ 
⎪ f[n+1][i] - f[n][i]
⎪ ——————————————————— = (1/6)k₁ + (2/6)k₂ + (2/6)k₃ + (1/6)k₄
⎩          Δt

時間離散化圖表稱作「Butcher tableau」。

multistep method (AB2):

    n     n-1    n-2   n-3
┌───────────────────────────┐
│   1                       │ Euler method
╞═══════════════════════════╡
│  3/2   -1/2               │ AB2
╞═══════════════════════════╡
│ 23/12 -16/12   5/12       │ AB3
╞═══════════════════════════╡
│ 55/24 -59/24  37/24 -9/24 │ AB4
└───────────────────────────┘

multistage method (RK4):

      shift amount
  Δt     k₁  k₂  k₃  k₄  
┌───────────────────────┐
│  0  │  0   0   0   0  │ round 1
│ 1/2 │ 1/2  0   0   0  │ round 2
│ 1/2 │  0  1/2  0   0  │ round 3
│  1  │  0   0   1   0  │ round 4
└─────┼─────────────────│
      │ 1/6 2/6 2/6 1/6 │ round weight
      └─────────────────┘

spatial discretization

空間離散化。考慮微分運算離散化的三種方式：

backward difference：後向差分。左邊。(f[i] - f[i-1]) / Δx
forward difference：前向差分。右邊。(f[i+1] - f[i]) / Δx
central difference：中央差分。中間。(f[i+1] - f[i-1]) / 2Δx

時間的微分運算離散化，無論採用哪種方式，計算過程都是一模一樣。我們只需要考慮空間的微分運算離散化。

實務上採用後向差分。想要減少數值誤差則換成中央差分。甚至兩者混用。花招百出，各有利弊。

舉例來說，下述微分方程式：

d           d
—— f(t,x) = —— f(t,x)
dt          dx

空間離散化。（時間離散化採用歐拉法。）

backward difference:

f(t+Δt,x) - f(t,x)   f(t,x) - f(t,x-Δx)
—————————————————— = ——————————————————
        Δt                   Δx

forward difference:

f(t+Δt,x) - f(t,x)   f(t,x+Δx) - f(t,x)
—————————————————— = ——————————————————
        Δt                   Δx

central difference:

f(t+Δt,x) - f(t,x)   f(t,x+Δx) - f(t,x-Δx)
—————————————————— = —————————————————————
        Δt                    2 Δx

改寫成陣列。

backward difference:

f[n+1][i] - f[n][i]   f[n][i] - f[n][i-1]
——————————————————— = ———————————————————
         Δt                    Δx

forward difference:

f[n+1][i] - f[n][i]   f[n][i+1] - f[n][i]
——————————————————— = ———————————————————
         Δt                    Δx

central difference:

f[n+1][i] - f[n][i]   f[n][i+1] - f[n][i-1]
——————————————————— = —————————————————————
         Δt                    2 Δx

空間離散化圖表稱作「stencil」。

backward difference:
    ┌───────┐
n+1 │   ·   │
n   │ · ·   │
    └───────┘
    i-1 i i+1
forward difference:
    ┌───────┐
    │   ·   │
    │   · · │
    └───────┘
central difference:
    ┌───────┐
    │   ·   │
    │ ·   · │
    └───────┘

implicit method / explicit method

時間索引值有兩種設定方式：

implicit method：隱方法。右式是下個時刻。左式右式都有未知項。
explicit method：顯方法。右式是當前時刻。僅左式有未知項。

實務上採用隱方法。想要省時則換成顯方法，但是可能不收斂。甚至兩者混用。順帶一提，兩者名稱源自隱方程式和顯方程式。

舉例來說，下述微分方程式：

d           d
—— f(t,x) = —— f(t,x)
dt          dx

時間離散化（歐拉法）、空間離散化（後向差分）。

implicit method:

f(t+Δt,x) - f(t,x)   f(t+Δt,x) - f(t+Δt,x-Δx)
—————————————————— = ————————————————————————
        Δt                      Δx

explicit method:

f(t+Δt,x) - f(t,x)   f(t,x) - f(t,x-Δx)
—————————————————— = ——————————————————
        Δt                   Δx

改寫成陣列。

implicit method:

f[n+1][i] - f[n][i]   f[n+1][i] - f[n+1][i-1]
——————————————————— = ———————————————————————
         Δt                     Δx

explicit method:

f[n+1][i] - f[n][i]   f[n][i] - f[n][i-1]
——————————————————— = ———————————————————
         Δt                   Δx

移項整理成遞迴公式。

implicit method:

            f[n+1][i-1] ⋅ Δt - f[n][i] ⋅ Δx
f[n+1][i] = ———————————————————————————————
                        Δt - Δx

explicit method:

                       Δt
f[n+1][i] = f[n][i] + ———— (f[n][i] - f[n][i-1])
                       Δx

隱方法是標準方法。優點是正確，缺點是難以移項整理成遞迴公式。甚至無法移項整理，必須解方程式，導致時間複雜度升高。甚至無法求解，只好放棄隱方法、採用顯方法。

顯方法是投機取巧的方法。優點是一定可以移項整理成遞迴公式、不必解方程式，缺點是可能不收斂。

備註

temporal discretization和spatial discretization不是一致公認的名稱，不過由於容易記憶，近年來逐漸風行。以下補充其他名稱。

時間步進的步伐大小與斜率：temporal discretization、time-stepping discretization、time integration，以及上述單字重新組合。其方法稱作time-stepping scheme、numerical integrator。

時間暨空間的取樣與擇鄰：spatial discretization、domain discretization。

numerical iteration

數值遞推。套用遞迴公式，逐步求得數值解。考慮三個項目：

(1) boundary：邊界。遞迴公式計算範圍。
(2) initial condition：初始條件。遞迴公式已知數值。時間變數邊界數值。
    boundary condition：邊界條件。遞迴公式已知數值。空間變數邊界數值。
(3) boundary scheme：邊界方案。遞迴公式邊界處理。

以及兩個技巧：

iteration matrix：遞推矩陣。遞迴公式恰是一次遞迴函數，可以改寫成矩陣運算。
local linearization：局部線性化。遞迴公式強行近似成一次遞迴函數。

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 condition / boundary condition

動態系統當中，變數被分成時間變數、空間變數，邊界被分成初始值、邊界值。

initial value：初始值。時刻是t₀的時間變數邊界。
boundary value：邊界值。空間變數邊界。

動態系統當中，問題被分成初始值問題、邊界值問題。初始值問題較簡單，演算法如前所述。邊界值問題較困難，演算法更加複雜，不屬於本章範疇。

initial value problem：初始值問題。邊界值不是預定義數值，取決於鄰居數值。
boundary value problem：邊界值問題。邊界值是預定義數值。

初始值視實際問題而定。例如零函數、脈衝函數、隨機函數。

(1) zero initial condition         零初始條件（一律是零）
(2) delta-pulse initial condition  脈衝初始條件（一處非零）
(3) random initial condition       隨機初始條件（亂數）

最難纏的初始值是隨機函數。必須謹慎內插適當數值，以確保數值連續平滑。一種替代方案是smooth noise：大量弦波疊加，頻率呈倍數，振幅為隨機數，相位是零。

邊界值不是預定義數值，稱作初始值問題。邊界值視實際問題而定。例如無界函數、對稱函數、週期函數。

(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：亂槍打鳥、試誤法。邊界值問題簡化為初始值問題。嘗試各種初始值（發射地點），各自進行數值遞推（行進軌道），檢查邊界值是否恰好符合預定義數值（恰好射中邊界值）。

實際問題還有各式各樣的初始條件與邊界條件，但是目前尚未形成一套知識體系，無法為大家介紹。例如零散時刻零散地點已知數值、外因即時影響數值。等你發表論文。

boundary scheme

邊界附近，遞迴公式變數索引值超出邊界。

對稱邊界條件、週期邊界條件，只需微調索引值，取絕對值、取餘數。開放邊界條件，需要大動土木，兩種解法：

one-side difference：單邊差分。縮減遞迴公式，變數索引值只剩邊界內側。
ghost point：鬼點。擴大表格外圍，變數索引值可達邊界外側。

單邊差分。左邊界用前向差分、右邊界用後向差分。缺點是數值解跑偏。

鬼點。擴大表格外圍。數量：源自遞迴公式右式，又源自差分公式。數值：源自邊界內側數值。

表格外圍沒有一致公認的名稱。此處採用ghost point。其他名稱：extra/virtual/ghost接sample/point/cell/element。

ghost point需要設定數值。兩種解法：

一、內插／迴歸。例如多項式內插、弦波迴歸。需要指定內插函數／迴歸函數。缺點是不合邏輯。數值解理應取決於遞迴公式，不應取決於自行假定的內插函數／迴歸函數。你有經驗根據另當別論。

二、近似／外插。例如泰勒近似、歐拉外插。需要指定導數值，即是追加邊界條件Neumann boundary condition。缺點是不符現實。我們哪知道實際的導數值。你有辦法測量導數值另當別論。

撇開邊界不談，平行計算也需要ghost point。空間分割成多個區塊，各個區塊各自計算。想要計算一個區塊的邊界數值，需要擷取隔壁區塊的邊界數值。ghost point是隔壁區塊的邊界數值。

備註

內插interpolation、外插extrapolation，兩者都是在推測函數值。前者針對已知函數值之間，後者針對已知函數值之外。

time-stepping method就是外插，推測下個時刻的數值解。

ghost point也是外插，推測邊界外側的數值解。

最初的、最經典的外插演算法，就是先前章節提到的Euler method。又稱作Euler extrapolation。

由於這些概念全是外插，促使大家另創名稱。數值模擬time-stepping method等等、時間離散化temporal discretization等等、表格外圍ghost point等等。天花亂墜，其實都是外插。

iteration matrix
（matrix formulation of the numerical scheme）

當遞迴公式恰是線性遞迴函數，可以改寫成矩陣A。

當遞迴公式恰是一次遞迴函數，再追加向量b。

implicit method:
f̂⁽ⁿ⁾ = A f̂⁽ⁿ⁺¹⁾ + b     (resolution)

explicit method:
f̂⁽ⁿ⁺¹⁾ = A f̂⁽ⁿ⁾ + b     (evaluation)

隱方法，習慣改寫成矩陣求解。形成一次方程組。

顯方法，無須改寫成矩陣求值。計算時間反而變長。

（簡單起見，我用顯方法作為主要範例。）

對稱邊界條件、週期邊界條件，只需要調整矩陣數值。開放邊界條件，還需要調整矩陣尺寸，以囊括表格外圍。

（簡單起見，我寫成負數索引值。請自行改成非負數索引值。）

3-point explicit scheme in polynomial formulation:

f[n+1][i] = α f[n][i-1] + β f[n][i] + γ f[n][i+1]

3-point explicit scheme in matrix formulation:

f̂⁽ⁿ⁺¹⁾ = A f̂⁽ⁿ⁾

open boundary condition:

⎡ f[n+1][-1]   ⎤   ⎡  │                 │  ⎤⎡ f[n][-1]   ⎤
⎢──────────────⎥   ⎢──┼─────────────────┼──⎥⎢────────────⎥
⎢ f[n+1][0]    ⎥   ⎢ γ│β α              │  ⎥⎢ f[n][0]    ⎥
⎢ f[n+1][1]    ⎥   ⎢  │γ β α            │  ⎥⎢ f[n][1]    ⎥
⎢ f[n+1][2]    ⎥   ⎢  │  γ β α          │  ⎥⎢ f[n][2]    ⎥
⎢ f[n+1][3]    ⎥   ⎢  │    γ β α        │  ⎥⎢ f[n][3]    ⎥
⎢      :       ⎥ = ⎢  │        ⋱        │  ⎥⎢    :       ⎥
⎢ f[n+1][♯x-4] ⎥   ⎢  │        γ β α    │  ⎥⎢ f[n][♯x-4] ⎥
⎢ f[n+1][♯x-3] ⎥   ⎢  │          γ β α  │  ⎥⎢ f[n][♯x-3] ⎥
⎢ f[n+1][♯x-2] ⎥   ⎢  │            γ β α│  ⎥⎢ f[n][♯x-2] ⎥
⎢ f[n+1][♯x-1] ⎥   ⎢  │              γ β│α ⎥⎢ f[n][♯x-1] ⎥
⎢──────────────⎥   ⎢──┼─────────────────┼──⎥⎢────────────⎥
⎣ f[n+1][♯x]   ⎦   ⎣  │                 │  ⎦⎣ f[n][♯x]   ⎦

 (the first row and the last row can be anything, e.g. zeros)
 (since they are unrelated to calculation)

 where f[n][-1] is extrapolated from f[n][0...k]
       f[n][♯x] is extrapolated from f[n][♯x-1...♯x-k]

symmetric boundary condition:
                           
⎡ f[n+1][0]    ⎤   ⎡ β z               ⎤⎡ f[n][0]    ⎤
⎢ f[n+1][1]    ⎥   ⎢ γ β α             ⎥⎢ f[n][1]    ⎥
⎢ f[n+1][2]    ⎥   ⎢   γ β α           ⎥⎢ f[n][2]    ⎥
⎢ f[n+1][3]    ⎥   ⎢     γ β α         ⎥⎢ f[n][3]    ⎥
⎢      :       ⎥ = ⎢         ⋱         ⎥⎢    :       ⎥
⎢ f[n+1][♯x-4] ⎥   ⎢         γ β α     ⎥⎢ f[n][♯x-4] ⎥
⎢ f[n+1][♯x-3] ⎥   ⎢           γ β α   ⎥⎢ f[n][♯x-3] ⎥
⎢ f[n+1][♯x-2] ⎥   ⎢             γ β α ⎥⎢ f[n][♯x-2] ⎥
⎣ f[n+1][♯x-1] ⎦   ⎣               z β ⎦⎣ f[n][♯x-1] ⎦

 (let z = α+γ)

  where f[n][-1] = f[n][1]
        f[n][♯x] = f[n][♯x-2]

periodic boundary condition:

⎡ f[n+1][0]    ⎤   ⎡ β α             γ ⎤⎡ f[n][0]    ⎤
⎢ f[n+1][1]    ⎥   ⎢ γ β α             ⎥⎢ f[n][1]    ⎥
⎢ f[n+1][2]    ⎥   ⎢   γ β α           ⎥⎢ f[n][2]    ⎥
⎢ f[n+1][3]    ⎥   ⎢     γ β α         ⎥⎢ f[n][3]    ⎥
⎢      :       ⎥ = ⎢         ⋱         ⎥⎢    :       ⎥
⎢ f[n+1][♯x-4] ⎥   ⎢         γ β α     ⎥⎢ f[n][♯x-4] ⎥
⎢ f[n+1][♯x-3] ⎥   ⎢           γ β α   ⎥⎢ f[n][♯x-3] ⎥
⎢ f[n+1][♯x-2] ⎥   ⎢             γ β α ⎥⎢ f[n][♯x-2] ⎥
⎣ f[n+1][♯x-1] ⎦   ⎣ α             γ β ⎦⎣ f[n][♯x-1] ⎦

 (either the first row or the last row is redundant)

  where f[n][0] = f[n][♯x-1]

表格外圍是已知數值，不符合一次方程組的標準格式。表格外圍視作區塊矩陣，挪至向量b，改寫成標準格式。

⎡ y₀ ⎤   ⎡ A₀₀ │ A₀₁ │ A₀₂ ⎤ ⎡ x₀ ⎤
⎢────⎥   ⎢─────┼─────┼─────⎥ ⎢────⎥
⎢ y₁ ⎥ = ⎢ A₁₀ │ A₁₁ │ A₁₂ ⎥ ⎢ x₁ ⎥   block matrix
⎢────⎥   ⎢─────┼─────┼─────⎥ ⎢────⎥
⎣ y₂ ⎦   ⎣ A₂₀ │ A₂₁ │ A₂₂ ⎦ ⎣ x₂ ⎦

y₁ = A₁₀x₀ + A₁₁x₁ + A₁₂x₂            middle row

y₁ = A₁₁x₁ + (A₁₀x₀ + A₁₂x₂)          y = Ax + b
     ^^^^^   ^^^^^^^^^^^^^^^
       Ax           b

where x₀ and x₂ are known values that are extrapolated by x₁

local linearization

當遞迴公式不是一次遞迴函數，可以強行近似成一次遞迴函數。

遞迴公式的非一次函數們，實施泰勒近似，逐個改成一次函數。

微分運算離散化：未知函數的微分運算的泰勒近似。局部線性化：已知函數的泰勒近似。接著繼續做，做得更徹底。

優點是讓隱方法得以求解。缺點是每回合都要重新計算矩陣。

程式碼

現在馬上來一段程式碼，讓各位實際觀察numerical scheme和numerical iteration是怎麼回事。

零維空間微分方程式。implicit method。

d
—— f(t) = f(t) + 1
dt

(f(t+Δt) - f(t)) / Δt = f(t+Δt) + 1
(f[n+1] - f[n]) / Δt = f[n+1] + 1
f[n+1] - f[n] = f[n+1] ⋅ Δt + Δt
f[n+1] - f[n+1] ⋅ Δt = f[n] + Δt
f[n+1] ⋅ (1 - Δt) = f[n] + Δt
f[n+1] = (f[n] + Δt) / (1 - Δt)

零維空間微分方程式。explicit method。

d
—— f(t) = f(t) + 1
dt

(f(t+Δt) - f(t)) / Δt = f(t) + 1
(f[n+1] - f[n]) / Δt = f[n] + 1
(f[n+1] - f[n]) = (f[n] + 1) ⋅ Δt
f[n+1] = f[n] + (f[n] + 1) ⋅ Δt

零維空間微分方程組。implicit method。

Lotka–Volterra equation:
⎰ d/dt x(t) = 1.1 x(t) - 0.4 x(t) y(t)
⎱ d/dt y(t) = 0.1 x(t) y(t) - 0.4 y(t)

⎰ (x[n+1] - x[n]) / Δt = 1.1 x[n+1] - 0.4 x[n+1] y[n+1]
⎱ (y[n+1] - y[n]) / Δt = 0.1 x[n+1] y[n+1] - 0.4 y[n+1]

離散化之後，形成方程組，需要求解。

一次方程組：手工推導反矩陣、高斯喬登消去法、鬆弛法。
多項式方程組：Gröbner basis。
一般的方程組：自求多福。

此例是特殊的方程組，可以簡化為一元二次多項式方程式。雖然得到兩個候選解，但是無法確認哪一個是正解。死局了。

由於難以求解，只好放棄隱方法、採用顯方法。

零維空間微分方程組。explicit method。

Lotka–Volterra equation:
⎰ d/dt x(t) = 1.1 x(t) - 0.4 x(t) y(t)
⎱ d/dt y(t) = 0.1 x(t) y(t) - 0.4 y(t)

⎰ (x[n+1] - x[n]) / Δt = 1.1 x[n] - 0.4 x[n] y[n]
⎱ (y[n+1] - y[n]) / Δt = 0.1 x[n] y[n] - 0.4 y[n]

一維空間微分方程式。explicit method。

constant-coefficient advection equation:
∂                ∂        
—— f(t,x) = -0.4 —— f(t,x)
∂t               ∂x

(f(t+Δt,x) - f(t,x)) / Δt = -0.4 (f(t,x) - f(t,x-Δx)) / Δx
(f[n+1][i] - f[n][i]) / Δt = -0.4 (f[n][i] - f[n][i-1]) / Δx
f[n+1][i] = f[n][i] - 0.4 (f[n][i] - f[n][i-1]) ⋅ Δt / Δx
f[n+1][i] = f[n][i] - C (f[n][i] - f[n][i-1])
where C = 0.4 (Δt / Δx)

二維空間微分方程式。explicit method。

heat equation:
∂                   ∂²             ∂²
—— f(t,x,y) = 0.1 ( ——— f(t,x,y) + ——— f(t,x,y) )
∂t                  ∂x²            ∂y²

空間離散化：中央版本的laplace。
(f[n+1][x][y] - f[n][x][y]) / Δt
  = 0.1 (f[n][x][y+1] + f[n][x][y-1] + f[n][x+1][y]
  + f[n][x-1][y] - 4 f[n][x][y]) / Δx / Δy

移項整理
f[n+1][i][j] = f[n][i][j] + 0.1 (f[n][i][j+1] + f[n][i][j-1]
  + f[n][i+1][j] + f[n][i-1][j] - 4 f[n][i][j]) ⋅ Δt / Δx / Δy

兩個陣列fₙₑₓₜ[i][j]和f[i][j]輪流使用。
fₙₑₓₜ[i][j] = f[i][j] + 0.1 (f[i][j+1] + f[i][j-1]
  + f[i+1][j] + f[i-1][j] - 4 f[i][j]) ⋅ Δt / Δx / Δy

簡寫成C和laplace，比較清爽。
fₙₑₓₜ[i][j] = f[i][j] + C laplace(f[i][j])
where C = 0.1 (Δt / Δx / Δy)

const int X = 5, Y = 5, T = 10;
float array1[X+2][Y+2];			// 兩個陣列輪流使用
float array2[X+2][Y+2];

float (*f)[Y+2] = array1;		// 解的一小部分（當前時刻）
float (*fnext)[Y+2] = array2;	// 解的一小部分（下個時刻）

float Δx = 1.0, Δy = 1.0, Δt = 1.0;	// 取樣間距
float v = 0.1;						// 擴散速度

void heat_equation()
{
	// initial condition
	memset(f, 0, sizeof(array1));	// 最初一片虛空
	f[1][1] = 1;					// 左上角熱量1

// boundary condition
	for (int i=0; i<=X+1; ++i) f[i][0] = f[i][Y+1] = 0;
	for (int j=0; j<=Y+1; ++j) f[0][j] = f[X+1][j] = 0;
	for (int i=0; i<=X+1; ++i) fnext[i][0] = fnext[i][Y+1] = 0;
	for (int j=0; j<=Y+1; ++j) fnext[0][j] = fnext[X+1][j] = 0;

draw_array(f);	// 一幕一幕，畫成動畫。

// 沿著主軸方向填充陣列
	for (int n=1; n<=T; ++n)
	{
		// 熱擴散，侷限於鄰居。
		for (int i=1; i<=X; ++i)
			for (int j=1; j<=Y; ++j)
			{
				float C = v * Δt / Δx / Δy;
				float laplace = f[i][j-1] + f[i][j+1]
							  + f[i-1][j] + f[i+1][j]
							  - 4 * f[i][j];
				fnext[i][j] = f[i][j] + C * laplace;
			}

swap(f, fnext);	// 兩個陣列輪流使用
		draw_array(f);	// 一幕一幕，畫成動畫。
	}
}

void draw_array(float f[X+2][Y+2])
{
	for (int i=1; i<=X; i++)
	{
		for (int j=1; j<=Y; j++)
			cout << setiosflags(ios::fixed)
				 << setprecision(3)
				 << f[i][j] << ' ';
		cout << '\n';
	}
	cout << '\n';
}

二維空間微分方程式。implicit method。

heat equation:
∂                   ∂²             ∂²
—— f(t,x,y) = 0.1 ( ——— f(t,x,y) + ——— f(t,x,y) )
∂t                  ∂x²            ∂y²

隱方法。
fₙₑₓₜ[i][j] = f[i][j] + C laplace(fₙₑₓₜ[i][j])

簡寫成sum，比較清爽。
fₙₑₓₜ[i][j] = f[i][j] + C (sumₙₑₓₜ - 4 fₙₑₓₜ[i][j])

移項整理，新值通通挪至左式。
(1 + 4C) fₙₑₓₜ[i][j] - C sumₙₑₓₜ = f[i][j]

一次方程組 A fₙₑₓₜ = b，已知 A 和 b，求解 fₙₑₓₜ。
採用鬆弛法，得到計算公式。
fₙₑₓₜ[i][j] = (f[i][j] + C sum) / (1 + 4C)

const int X = 5, Y = 5, T = 10;
float array1[X+2][Y+2];			// 兩個陣列輪流使用
float array2[X+2][Y+2];

float (*f)[Y+2] = array1;		// 解的一小部分（當前時刻）
float (*fnext)[Y+2] = array2;	// 解的一小部分（下個時刻）

float Δx = 1.0, Δy = 1.0, Δt = 1.0;	// 取樣間距
float v = 0.1;						// 擴散速度

void heat_equation()
{
	// initial condition
	memset(f, 0, sizeof(array1));	// 最初一片虛空
	f[1][1] = 1;					// 左上角熱量1

draw_array(f);	// 一幕一幕，畫成動畫。

for (int n=1; n<=T; ++n)
	{
		solver(f, fnext);	// 解一次方程組
		swap(f, fnext);		// 兩個陣列輪流使用
		draw_array(f);		// 一幕一幕，畫成動畫。
	}
}

// Gauss–Seidel iteration
void solver(float f[X+2][Y+2], float fnext[X+2][Y+2])
{
	// 初始化，理論上可以是任意值。
//	for (int i=1; i<=X; ++i)
//		for (int j=1; j<=Y; ++j)
//			fnext[i][j] = f[i][j];

// 意思意思遞推100回合。熱擴散，不再侷限於鄰居。
	for (int k=0; k<100; ++k)
		for (int i=1; i<=X; ++i)
			for (int j=1; j<=Y; ++j)
			{
				float C = v * Δt / Δx / Δy;
				float sum = fnext[i][j-1] + fnext[i][j+1]
						  + fnext[i-1][j] + fnext[i+1][j];
				fnext[i][j] = (f[i][j] + sum * C)
							/ (1 + 4 * C);
			}
}

differential equation - dynamical system🚧

numerical error

數值誤差。數值遞推過程，符號解與數值解的差距。

symbolic solution / numeric solution

符號解、數值解，數學符號的寫法。

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ᵢ)

local error / global 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。global error當中，一步的誤差持續影響往後每一步。

error metric

局部誤差、全域誤差的距離函數，通常設定成平方距離。

兩個函數的距離：黎曼積分。

‖f⁽ⁿ⁾ - g⁽ⁿ⁾‖₂ = | f(tₙ) - g(tₙ) |                    0D space

‖f⁽ⁿ⁾ - g⁽ⁿ⁾‖₂ = √ ∫ (f(tₙ,x) - g(tₙ,x))² dx          1D space

‖f⁽ⁿ⁾ - g⁽ⁿ⁾‖₂ = √ ∫∫ (f(tₙ,x,y) - g(tₙ,x,y))² dxdy   2D space

兩個離散化函數的距離：黎曼和。

我們習慣在邊界取樣，導致黎曼和在左右邊界各多出½Δx寬度。無傷大雅。如果你很在意，你可以扣掉它。

‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖₂ = √ sum (f̂⁽ⁿ⁾(xᵢ) - ĝ⁽ⁿ⁾(xᵢ))² Δx
               = √ sum (f[i] - g[i])² Δx

一個函數與一個離散化函數的距離：黎曼和，事情比較簡單。

‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖₂ = √ sum (f⁽ⁿ⁾(xᵢ) - f̂⁽ⁿ⁾(xᵢ))² Δx

numerical convergence

數值收斂。數值遞推過程，全域誤差趨近零，答案精準。

中文譯作「收斂」，原意應是「合流」。

數值收斂的數學定理的進化過程：

(1) problem: series convergence
    convergence: monotone convergence theorem
                 (if monotone, then bounded = convergent)

(2) problem: equation solving
    algorithm: fixed-point iteration
    convergence: Banach fixed-point theorem
                 (if Lipschitz, 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)

discretization of differential equation
L̂(f(t,x)) = 0

numeric 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̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ 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̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ → 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

一致性：函數方程式、其離散化（遞迴公式），兩者趨近相等。
穩定性：兩組數值解（初始條件不同），各回合距離受限。
一致性且穩定性：符號解、數值解，各回合距離受限。
收斂性：兩組數值解（初始條件不同），最終趨近相等。
一致性且收斂性：符號解、數值解，各回合趨近相等。

一致性：截斷誤差趨近零。
一致性且穩定性：全域誤差受限。
一致性且收斂性：全域誤差趨近零。

有些文獻又進一步修改定義，只討論動態系統d/dt f(t,x) = F(f(t,x))，只針對右式F(f(t,x))。

dynamical system
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))，事情更加棘手，我沒有研究。

如此一來，穩定性宛如平緩函數Lipschitz function。數值分析基本定理宛如不動點定理Banach fixed-point theorem：平緩函數套用不動點遞推法收斂至不動點。請見本站文件「fixed point」。

數值分析基本定理的證明，我還沒有學會。等我學完泛函分析不知何年何月。關鍵字也許是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 (with first-order accuracy)
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 (with first-order accuracy)
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)

first-order central difference (with second-order accuracy)
⎰ 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′[i] = (f[i+1] - f[i-1]) / (2 Δx) + O((Δx)²)

二次微分的差分公式。此處提供中央版本。

second-order central difference (with second-order accuracy)
⎰ 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)²)

三次微分的差分公式。若有興趣請自行深究。

https://en.wikipedia.org/wiki/Finite_difference_coefficient

大家習慣使用大O符號呈現截斷誤差。大家習慣用Δx次方值作為截斷誤差的數量級。

nth-order error      n階誤差：截斷誤差是O((Δx)ⁿ)
nth-order accuracy   n階精度：截斷誤差是O((Δx)ⁿ)

負面詞彙「誤差error」、正面詞彙「精度accuracy」，兩者概念相同，都有人用。

注意到，n階差分，可能是指二次微分，可能是指二階精度。有些文獻將微分次數放在前綴，有些文獻將精度數量級放在前綴。

泰勒級數各項，分母是階乘、分子是指數。階乘成長速度遠大於指數成長速度。Δx次方值越高，該項數值越小。泰勒級數保留越多低次方項（消滅越少高次方項），截斷誤差越小。

想讓截斷誤差趨近零，一種方式是採用解析函數。如此一來，泰勒級數每一項通通存在、而且通通趨近零。無論消滅多少高次方項，最後保留下來的低次方項，仍然通通趨近零。

想讓截斷誤差趨近零，也可以採用其他函數。找到截斷誤差的數學式子，證明其趨近零。比較麻煩就是了。

consistency of numerical simulation
（numerical consistency）

數值模擬當中的一致性：當取樣間距Δt與Δx趨近零，微分方程式、其離散化（遞迴公式），兩者相符。

注意到，數學家只有討論「Δt與Δx趨近零則趨近相符」，並未討論「Δt與Δx越微小則越相符」。不太務實就是了。

注意到，數學家只有考慮截斷誤差，並未考慮四捨五入誤差。馬馬虎虎就是了。

微分方程式重新整理成泛函數的模樣。微分方程式、其離散化（遞迴公式），兩個泛函數相減得到差距O((Δt)ᵖ + (Δx)𐞥)。

differential equation:
d/dt f(t,x) = d/dx f(t,x)
d/dt f(t,x) - d/dx f(t,x) = 0

discretization:
  (f(t+Δt,x) - f(t,x)) / Δt + O((Δt)ᵖ)
= (f(t,x+Δx) - f(t,x)) / Δx + O((Δx)𐞥)

  (f(t+Δt,x) - f(t,x)) / Δt - (f(t,x+Δx) - f(t,x)) / Δx
= O((Δt)ᵖ + (Δx)𐞥)

functional:
let L(f(t)) = d/dt f(t,x) - d/dx f(t,x)
let L̂(f(t)) = (f(t+Δt,x) - f(t,x)) / Δt
            - (f(t,x+Δx) - f(t,x)) / Δx

consistency:
L(f(t,x)) =    lim    L̂(f(t,x))
            ^{Δt→0,Δx→0}
   lim    { L(f(t,x)) - L̂(f(t,x)) } = 0
^{Δt→0,Δx→0}
L(f(t,x)) - L̂(f(t,x)) → 0 when Δt→0, Δx→0

truncation error:
L(f(t,x)) - L̂(f(t,x)) = O((Δt)ᵖ + (Δx)𐞥)

有些文獻又進一步修改定義，只討論動態系統d/dt f(t,x) = F(f(t,x))，只針對右式F(f(t,x))。

dynamical system:
d/dt f(t,x) = F(f(t,x))

consistency:
F(f(t,x)) =    lim    F̂(f(t,x))
            ^{Δt→0,Δx→0}
   lim    { F(f(t,x)) - F̂(f(t,x)) } = 0
^{Δt→0,Δx→0}
F(f(t,x)) - F̂(f(t,x)) → 0 when Δt→0, Δx→0

truncation error, that equals local error:
F(f(t,x)) - F̂(f(t,x)) = O((Δx)𐞥)

時間的微分運算離散化，無論採用哪種方式，計算過程都是一模一樣。左式採用一階前向差分，即是p = 1。因此只針對右式。

「數值代入右式」即是「遞推一步」。截斷誤差即是局部誤差！O((Δx)𐞥)得以描述‖f⁽¹⁾ - f̂⁽¹⁾‖。

詳細來說，動態系統右式的各個微分運算的截斷誤差，組成了右式的截斷誤差，即是遞推一步的局部誤差。我們可以使用大O符號描述遞推一步的誤差／精度。各種遞迴公式擁有各種誤差／精度。

截斷誤差趨近零，那麼局部誤差也趨近零。想讓局部誤差趨近零，一種方式是符號解f(t,x)與右式F(f(t,x))都採用解析函數。

stability

「穩定性」。兩個數列距離受限。

穩定性細分為許多類型。

Lyapunov stability：內部穩定性。距離小於定值，定值取決於最初距離。
asymptotic stability：漸進穩定性。距離收斂至零。
exponential stability：指數穩定性。距離指數衰減。

https://en.wikipedia.org/wiki/Lyapunov_stability

stability of numerical simulation
（numerical stability）

數值模擬當中的穩定性：（兩種初始條件所得到的）兩種數值解每回合距離受限。

穩定性細分為許多類型。

bounded stability：受限穩定性。每回合距離小於等於定值。
Lipschitz stability：平緩穩定性。每回合距離小於等於最初距離的特定倍率。
contractive stability：收縮穩定性。每回合距離小於等於最初距離。

https://en.wikipedia.org/wiki/Contraction_mapping

數值分析基本定理，採用受限穩定性。針對動態系統，換成平緩穩定性。工程數學教科書習慣挑選簡易的微分方程式作為教學範例（指數方程式、守恆型），並且再換成收縮穩定性。

數值分析基本定理，穩定性還能追加等於。先前章節已經提及。

fundamental theorem of numerical analysis特殊版本

數值分析基本定理，符號解版本。

(1) Peano existence theorem
    常微分方程ODE、初始值問題IVP，一定存在解。

(2) Cauchy—Lipschitz theorem
    常微分方程ODE、初始值問題IVP，若是平緩函數，則存在唯一解。

(3) Cauchy—Kovalevskaya theorem
    偏微分方程PDE、初始值問題IVP，若是平緩函數，則存在唯一解。

數值分析基本定理，數值解版本：

(1) Lax–Richtmyer equivalence theorem
    初始值問題IVP、線性微分方程式、解析解
    線性方案、Lₚ-diminishing穩定性、Lₚ收斂性

(2) Dahlquist's equivalence theorem
    初始值問題IVP、齊次微分方程式、解析解
    multistep method、根條件、Lₚ收斂性

(3) Lax–Wendroff consistency theorem & [DiPerna 1983]
    初始值問題IVP、守恆律、弱解
    一致守恆方案、bounded L_∞暨bounded variation穩定性、L₁收斂性

stability analysis

接下來介紹實務上經常使用的數學理論。

穩定性有三種類型。

zero-stability：只考慮遞迴公式暨數值解的穩定性。
              　當Δt與Δx趨近零，滿足穩定性。
absolute stability：只考慮遞迴公式暨數值解的穩定性。
                  　特定大小的Δt與Δx，滿足穩定性。
A-stability：也考慮微分方程式暨符號解的穩定性。

穩定性分析。使用數學工具推導穩定性。

穩定性分析有三種手法。總之都是頻譜。

(1) stability analysis with matrix
    (Lax–Richtmyer stability analysis)
    線性遞迴公式，改寫成矩陣，特徵值絕對值皆小於等於1。

(2) stability analysis with Fourier transform
    (von Neumann stability analysis)
    遞迴公式，改寫成弦波，振幅絕對值皆小於等於1。

(3) stability analysis with characteristic equation
    (Dahlquist stability analysis)
    遞迴公式，改寫成特徵方程式，根絕對值皆小於等於1。

absolute stability

數值分析基本定理，要求Δt與Δx趨近零。然而，實務上浮點數位數有限、陣列大小有限，因此Δt與Δx不可能趨近零。實務上必須檢查特定大小的Δt與Δx是否滿足穩定性。嚴謹起見，數學家另造一詞，稱作絕對穩定性。簡單起見，大家還是稱作穩定性。

zero-stability：當Δt與Δx趨近零，滿足穩定性。
absolute stability：特定大小的Δt與Δx，滿足穩定性。

數值分析基本定理，採用受限穩定性或平緩穩定性。然而，實務上難以判斷受限穩定性或平緩穩定性。方便起見，大家換成更強更嚴格的收縮穩定性。

順帶一提，收縮穩定性得以消滅或抑制數值誤差。數值解變小或不變，數值誤差也隨之變小或不變。

contractive stability

「收縮穩定性」。數值解每回合距離小於等於最初距離。

‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖         contractive stability

當遞迴公式不因時刻而變，每回合的數值解都可以重新作為初始條件。收縮穩定性變成：數值解每回合距離變小或不變。

‖f̂⁽ⁿ⁺¹⁾ - ĝ⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖     time-invariant

當遞迴公式恰是線性遞迴函數，兩種數值解相加減仍是數值解。收縮穩定性變成：數值解每回合範數變小或不變。

‖f̂⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾‖                     linearity

exponential stability

「指數穩定性」。數值解每回合距離指數衰減或不變。

‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ exp(-c n) ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ exponential stability
c is a constant and 0 ≤ exp(-c n) ≤ 1

當遞迴公式不因時刻而變，指數穩定性變成：數值解每回合距離乘上固定倍率，小於等於1。

‖f̂⁽ⁿ⁺¹⁾ - ĝ⁽ⁿ⁺¹⁾‖ ≤ k ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖   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
（Lax–Richtmyer stability analysis）

當遞迴公式恰是線性遞迴函數，可以改寫成矩陣運算。

順帶一提，還有一種改寫方法，將多個回合的數值解併成一個向量，請見本站文件「linear recurrence」。該方法不適用於此處。

專論《Some insights into the stability of difference approximations for hyperbolic initial-boundary-value problems》。

differential equation:

d           d
—— f(t,x) = —— f(t,x)
dt          dx

explicit backward Euler method:
                      
(f[n+1][i] - f[n][i]) / Δt = (f[n][i] - f[n][i-1]) / Δx
f[n+1][i] = f[n][i] + (Δt/Δx) (f[n][i] - f[n][i-1])
f[n+1][i] = (1 + Δt/Δx) f[n][i] - (Δt/Δx) f[n][i-1]

linear recurrence in polynomial formulation:

f[n+1][i] = α f[n][i] + β f[n][i-1]     where α = 1 + Δt/Δx
                                              β = - Δt/Δx

linear recurrence in matrix formulation:

⎡ f[n+1][0]    ⎤   ⎡ α β 0 0 ... 0 0 0 0 ⎤⎡ f[n][0]    ⎤
⎢ f[n+1][1]    ⎥   ⎢ β α 0 0 ... 0 0 0 0 ⎥⎢ f[n][1]    ⎥
⎢ f[n+1][2]    ⎥   ⎢ 0 β α 0 ... 0 0 0 0 ⎥⎢ f[n][2]    ⎥
⎢ f[n+1][3]    ⎥   ⎢ 0 0 β α ... 0 0 0 0 ⎥⎢ f[n][3]    ⎥
⎢      :       ⎥ = ⎢ : : : :     : : : : ⎥⎢    :       ⎥
⎢      :       ⎥   ⎢ : : : :     : : : : ⎥⎢    :       ⎥
⎢ f[n+1][♯x-3] ⎥   ⎢ 0 0 0 0 ... β α 0 0 ⎥⎢ f[n][♯x-3] ⎥
⎢ f[n+1][♯x-2] ⎥   ⎢ 0 0 0 0 ... 0 β α 0 ⎥⎢ f[n][♯x-2] ⎥
⎣ f[n+1][♯x-1] ⎦   ⎣ 0 0 0 0 ... 0 0 β α ⎦⎣ f[n][♯x-1] ⎦

     f̂⁽ⁿ⁺¹⁾                   A                f̂⁽ⁿ⁾

(symmetric boundary condition)

當遞迴公式恰是線性遞迴函數，收縮穩定性改寫成矩陣運算。

‖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                contractive stability

順帶一提，當遞迴公式恰是線性遞迴函數，收縮穩定性顯然滿足收斂性，畢竟指數衰減至零。大可不必煩惱數值分析基本定理。

amplification factor

「假的放大因子」。下個時刻的數值解、當前時刻的數值解，兩者對應項相除。這種定義方式導致嚴重問題。

    f̂⁽ⁿ⁺¹⁾
a = ——————                            amplification factor ✘
     f̂⁽ⁿ⁾

假的放大因子與矩陣無關。假的放大因子與收縮穩定性無關。

       ║f̂⁽ⁿ⁺¹⁾║    ‖f̂⁽ⁿ⁺¹⁾‖₂
‖a‖₂ = ║——————║  ≠ ————————— ≤ ‖A‖₂ ≤ 1
       ║ f̂⁽ⁿ⁾ ║₂    ‖f̂⁽ⁿ⁾‖₂

當遞迴公式恰是線性遞迴函數，而且，微分方程式恰是零維空間（常微分方程式ODE），假的放大因子與矩陣皆退化成相同數值a。此時假的放大因子與收縮穩定性才勉強有關聯。

      │f̂⁽ⁿ⁺¹⁾│   |f̂⁽ⁿ⁺¹⁾|
|a| = │——————│ = ———————— = |a| ≤ 1   contractive stability
      │ f̂⁽ⁿ⁾ │    |f̂⁽ⁿ⁾|

「放大因子」。數值解實施傅立葉轉換得到頻譜。下個時刻的數值解頻譜、當前時刻的數值解頻譜，兩者對應項相除得到放大因子。這是比較合理的定義方式。

    DFT(f̂⁽ⁿ⁺¹⁾)
a = ———————————                       amplification factor
     DFT(f̂⁽ⁿ⁾)

放大因子的地位宛如先前小節的矩陣特徵值。

當遞迴公式恰是線性遞迴函數，放大因子的每一個數值的絕對值皆小於等於1，才會滿足收縮穩定性。【待補證明】

      │DFT(f̂⁽ⁿ⁺¹⁾)│   |DFT(f̂⁽ⁿ⁺¹⁾)|
|a| = │———————————│ = ————————————— ≤ 1   contractive stability
      │ DFT(f̂⁽ⁿ⁾) │    |DFT(f̂⁽ⁿ⁾)|

stability analysis with Fourier transform
（von Neumann stability analysis）

當遞迴公式恰是線性遞迴函數，可以逐項轉換到頻域。

舉例來說，下述微分方程式：

d             d
—— f(t,x) + k —— f(t,x) = 0     k is a parameter
dt            dx

遞迴公式：

f[n+1][i] = f[n][i] - ½ C (f[n][i+1] - f[n][i-1])
where C = k Δt / Δx is CFL number without absolute value

forward difference of t
central difference of x

傅立葉轉換：

DFT(f[n+1]) = DFT(f[n]) - ½ C (exp(+𝑖ωi) DFT(f[n])
                             - exp(-𝑖ωi) DFT(f[n]))
DFT(f[n+1]) = (1 - ½ C (exp(+𝑖ωi) + exp(-𝑖ωi)) DFT(f[n])

此遞迴公式無論如何不可能穩定。

a[i] = 1 - ½ C (exp(+𝑖ωi) - exp(-𝑖ωi))   amplification factor
a[i] = 1 - C 𝑖 sin(ωi)                   amplification factor
numerical scheme is stable iff |a[i]| ≤ 1
this numerical scheme is unstable since |a[i]| > 1

修改遞迴公式，採用Lax–Friedrichs method，就有機會穩定。

f[n+1][i] = ½ (f[n][i+1] + f[n][i-1])
          - ½ C (f[n][i+1] - f[n][i-1])

a[i] = ½ (exp(+𝑖ωi) + exp(-𝑖ωi))         amplification factor
     - ½ C (exp(+𝑖ωi) - exp(-𝑖ωi))
a[i] = cos(ωi) - C 𝑖 sin(ωi)             amplification factor

let |C| ≤ 1 is CFL condition
thus |a[i]|² = cos²(ωi) + C² sin²(ωi)
             ≤ cos²(ωi) + sin²(ωi)
             = 1
thus |a[i]| ≤ 1

a numerical scheme is stable iff |a[i]| ≤ 1
this numerical scheme is stable under CFL condition |C| ≤ 1

順帶一提，當空間邊界無窮大，只好分析一小段區間，以此代表一整片範圍。此時大家習慣採用週期邊界條件、並且採用連續傅立葉轉換。如此一來，一小段區間以外的地方，通通都是同理可證。

dispersion analysis

數值解拆解成多個弦波，其頻率呈倍數。下個時刻各個弦波的縮放比例，即是放大因子。

放大因子是多個複數，分別代表各個弦波。一個複數可以拆解成長度和角度，分別稱作強度magnitude和相位phase，進一步討論強度變化比例、相位變化比例。若有興趣請自行深究。

專著《Numerical Methods for Fluid Dynamics: with Applications to Geophysics》。

magnitude = sqrt(Re[a]² + Im[a]²) = |a|
phase = tan⁻¹(Im[a] / Re[a])

disturbance amplification

數值誤差源自四捨五入誤差、截斷誤差。當滿足穩定性，數值誤差被抑制在一定程度，數值解符合符號解。當不滿足穩定性，數值解突然出現一個特定頻率的弦波，振幅呈指數成長。這個現象稱作擾動放大。俗稱爆掉。

von Neumann stability analysis可以解釋這個現象。觀察放大因子的強度頻譜，強度頻譜一般來說呈尖峰狀。當不滿足穩定性，表示尖峰處的強度大於一，導致該頻率的弦波的振幅呈指數成長，成長速度領先其他頻率的弦波。

region of absolute stability

滿足絕對穩定性的情況下，推導Δt與Δx的所有範圍，叫做region of absolute stability，簡稱stable region。

舉例來說，下述微分方程式：

d
—— f(t) = f(t)
dt

套用implicit Euler method，選擇適當的Δt，就會穩定。

(f(t + Δt) - f(t)) / Δt = f(t + Δt)   implicit Euler method
(f[n+1]    - f[n]) / Δt = f[n+1]      array indexing
f[n+1] = f[n] / (1 - Δt)              recurrence

f[n+1] = a f[n]   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

套用explicit Euler method，無論如何不可能穩定。

(f(t + Δt) - f(t)) / Δt = f(t)        explicit Euler method
(f[n+1]    - f[n]) / Δt = f[n]        array indexing
f[n+1] = f[n] (1 + Δt)                recurrence

f[n+1] = a f[n]   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

套用implicit Euler method，將數值a換成函數ϕ(Δt)。

f[n+1] = a f[n]       where a = 1 / (1 - Δt)
f[n+1] = ϕ(Δt) f[n]   where ϕ(Δt) = 1 / (1 - Δt)

數學家甚至追加參數。微分方程式／遞迴公式附帶參數，判斷各種參數數值的穩定性。

舉例來說，下述微分方程式：

d
—— f(t) = k f(t)     k is a parameter
dt

套用implicit Euler method，將數值a換成函數ϕ(k Δt)。

f[n+1] = a f[n]         where a = 1 / (1 - k Δt)
f[n+1] = ϕ(k Δt) f[n]   where ϕ(k Δt) = 1 / (1 - k Δt)

k可以是複數，形成複數函數ϕ(z)。函數輸入是z = k Δt。

ϕ(z) = 1 / (1 - z)

數學家描述等高線，習慣使用符號ϕ。我也採用符號ϕ。

diagram of absolute stability

有些教科書喜歡把stable region畫出來。此圖毫無實際用處，只是想讓讀者轉換心情。

舉例來說，下述微分方程式：

d
—— f(t) = k f(t)     k is a parameter
dt

套用implicit Euler method，得到左圖。

套用explicit Euler method，得到右圖。

穩定區域|a| ≤ 1即是|ϕ(k Δt)| ≤ 1即是|ϕ(z)| ≤ 1畫在複數平面。座標軸是z = k Δt的實部虛部。

【圖片尚待修正】

https://personal.math.vt.edu/embree/math5466/nanotes.pdf

k可能是複數，於是畫在複數平面、視作二維向量。Δt只能是實數，於是視作向量縮放倍率。穩定區域|a| ≤ 1與射線k的交點，觀察縮放倍率，得到適當的Δt。

數學家利用ϕ(z)定義穩定性。此例當中：

zero-stability：穩定區域包含原點。
absolute stability：原點放射線可以碰到穩定區域、且不是碰到原點。
stiff stability：承上，當ϕ(z)趨近零，z趨近負無限大。

stiffness與stiff stability

「剛硬性」。微分方程式的造型太過特殊，選擇很小很小的Δt與Δx才能滿足穩定性，需要很多很多的遞推回合才能到達指定時刻，導致很久很久的執行時間才能跑完數值模擬。硬篤、歹舞。

「剛硬穩定性」。滿足絕對穩定性，而且，當Δt與Δx趨近無限大，則收斂速度趨近無限大（指數趨近負無限大）。

剛硬穩定性得以避免上述困擾。剛硬穩定性一點都不剛硬。

A-stability

先前只有考慮遞迴公式暨數值解的穩定性，現在還要考慮微分方程式暨符號解的穩定性。

微分方程式暨符號解的數學性質、遞迴公式暨數值解的數學性質，兩邊必須一致，包含穩定性。

「A穩定性」。一旦微分方程式暨符號解滿足穩定性，那麼遞迴公式暨數值解必須滿足絕對穩定性。

zero-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

甚至深入觀察特殊的穩定性。若有興趣請自行深究。

專著《Numerical Methods for Ordinary Differential Equations》。

A-stability：A穩定性。A是角度angle，詞不達意。
           　一旦微分方程式暨符號解滿足穩定性，
           　那麼遞迴公式暨數值解必須滿足絕對穩定性。
L-stability：L穩定性。L是左邊left，詞不達意。
           　一旦微分方程式暨符號解滿足穩定性，
           　那麼遞迴公式暨數值解必須滿足剛硬穩定性。
B-stability：B穩定性。B是作者Butcher，我猜的。
           　一旦微分方程式暨符號解滿足單調遞減性（等價於收縮穩定性），
           　那麼遞迴公式暨數值解必須滿足收縮穩定性。

甚至深入觀察經典的演算法的穩定性。若有興趣請自行深究。

G-stability：滿足A穩定性的情況下，
           　multistep method的數學性質。
algebraic stability：滿足A穩定性的情況下，
                   　multistage method的數學性質。
AN-stability：滿足B穩定性的情況下，
            　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 parameter
dt

符號解：exp(k t)收斂至零或受限，(k t)是負數或零，才會滿足穩定性。

f(t) = f(0) exp(k t)

數值解：特定的數值方案，得到特定的ϕ(z)。

f[n+1] = ϕ(z) f[n]   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
（Dahlquist stability analysis）

關鍵字linear difference equation。

遞迴公式：

f[n] = c₁ f[n-1] + ... + cₘ f[n-m]

特徵方程式：

λᵐ = c₁ λᵐ⁻¹ + ... + cₘ λᵐ⁻ᵐ

移項整理，左式是特徵多項式：

λᵐ - c₁ λᵐ⁻¹ - ... - cₘ λᵐ⁻ᵐ = 0

特徵多項式分解：

(x - λ₁)...(x - λₘ) = 0

遞迴公式改寫成一般公式：

f[n] = c₁ λ₁ⁿ + ... + cₘ λₘⁿ

想要檢查收縮穩定性，只需將遞迴公式改寫成特徵方程式，然後檢查所有根絕對值皆小於等於1。

根絕對值又細分為許多情況。總之都在單位圓裡面。

root condition：根條件。根絕對值皆小於等於1，而且絕對值等於1者皆不相同。
weak root condition：弱根條件。承上，而且絕對值等於1者至少兩個。
strong root condition：強根條件。承上，而且絕對值等於1者恰好一個。

https://sites.math.rutgers.edu/~falk/math573/lecture21.pdf

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[n] + ... + aₘ f[n-m])
　　　　　　　　= k Δt (b₀ f[n] + ... + bₘ f[n-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 theorem

此穩定性分析衍生許多定理。若有興趣請自行深究。

專論《Dahlquist's classical papers on stability theory》。

專著《Numerical Mathematics》。

(1) Dahlquist's equivalence theorem
    multistep method is convergent iff consistent and stable.
    multistep method is stable iff it satisfies root condition.
    [Dahlquist 1956]

(2) Dahlquist's stability theorem
    implicit multistep method is A-stable under certain condition.
    explicit multistep method is never A-stable.
    [Dahlquist 1963]

(3) Dahlquist's order barrier theorem
    implicit multistep method with m steps is at most
    m+1 (m is odd) / m+2 (m is even) order accurate.
    explicit multistep method with m steps is at most
    m order accurate. 
    [Dahlquist 1978]

(4) Dahlquist's order barrier theorem
    A-stable implicit multistep method is at most
    second order accurate.
    [Dahlquist 1963]

differential equation - dynamical system🚧

numerical simulation

數值模擬。微分方程式，求得數值解。四個步驟：

(1) numerical scheme
(2) numerical iteration
(3) numerical error
(4) numerical convergence

先前已經提過，只是複製貼上而已。以下針對每個項目，補充更多詞彙。二周目追加新要素，希望各位喜歡。

numerical scheme

數值方案。面對各式各樣的微分方程式，想方設法調整遞迴公式，滿足一致性、滿足穩定性，導致收斂性。

(1) temporal discretization
(2) spatial discretization
(3) implicit method / explicit method

temporal discretization

時間離散化也可以管得更細、管得更寬。牛屄。

(1) multistep method：多步法。取多處斜率。都是左方斜率。
(2) multistage method：多階段法。求多次斜率。都是右方斜率。
(3) operator splitting method：算子分裂法。項數太多，分成多步。
(4) multiderivative method：多導數法。步數太多，併成一步。

時間離散化也可以採用其他級數。雞肋。

(1) Euler extrapolation：時間離散化採用Taylor series。
(2) Adams extrapolation：時間離散化採用Maclaurin series。
(3) Richardson extrapolation：時間離散化採用geometric series。

https://mathworld.wolfram.com/AdamsMethod.html
https://en.wikipedia.org/wiki/Richardson_extrapolation

經典的基礎的方法。物理模擬經常提及。理都懂然并卵。

forward time method：單步法。亦是單階段法。
 1. 取當前時刻斜率，走一刻到下個時刻。（時間前向差分）
leapfrog method：多步法。
 1. 取當前時刻斜率，從上個時刻走兩刻到下個時刻。（時間中央差分）
predictor–corrector method：兩階段法。
 1. 取當前時刻斜率，走一刻到下個時刻。
 2. 取當前時刻與下個時刻的斜率平均數，從頭重新走一刻。
midpoint method：兩階段法。
 1. 取當前時刻斜率，走半刻到中點時刻。
 2. 取中點時刻斜率，從頭重新走一刻。

隱性多階段法的重要流派。超潮正夯。

implicit Runge–Kutta method
 └ linearly implicit Runge–Kutta method
    └ diagonally implicit Runge–Kutta method

operator splitting method

算子分裂法。分割步驟。項數太多，逐步更新。

多步法、多階段法針對右式整體。算子分裂法針對右式各項。

(1) operator splitting：第一項從頭走一刻，第二項從頭走一刻。
(2) Strang splitting：第一項從頭走半刻，第二項從中點走一刻，第一項從中點走半刻。

operator splitting。第一項從頭走一刻，第二項從頭走一刻。
     d/dt f(t) = a(t) + b(t)
---> 1. (f[n+1] - f[n]) / Δt = a(t)
     2. (f[n+1] - f[n]) / Δt = b(t)

Strang splitting。第一項從頭走半刻、第二項從中點走一刻、第一項從中點走半刻。
     d/dt f(t) = a(t) + b(t)
---> 1. (f[n+1] - f[n]) / ½Δt = a(t)
     2. (f[n+1] - f[n]) / Δt  = b(t+½Δt)
     3. (f[n+1] - f[n]) / ½Δt = a(t+½Δt)

https://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book018.html

multiderivative method

多導數法。時間離散化也可以採用高階差分公式。

多步法、多階段法的數學式子強行展開，宛如高階差分。高階差分主要有兩個功用。

一、偷食步。一階差分重複n回合，其實就是n階差分、Δt與Δx變成n倍。反過來說，n階差分，其實就是Δt與Δx變成1/n倍、一口氣做n次差分。

二、微調數值解。取大量斜率的加權平均，讓數值解更符合符號解。強行修改遞迴公式的係數，嘗試通靈答案。

補個範例

spatial discretization

專著《Fundamentals of Engineering Numerical Analysis》。

(1) higher-order accuracy finite difference：追加項數，升高精度。
(2) compact finite difference：減少項數，維持精度。
(3) modified wave number：修改係數，達成絕對一致性。

higher-order accuracy finite difference：追加項數，逐項使用泰勒近似，通通加總湊出係數，形成線性方程組。推導過程在此。

也可以視作多項式內插，然後微分。大家習慣使用Lagrange interpolation。推導過程在此。

first-order central difference (with 4th-order accuracy)
f′[i] = (- 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：使用泰勒近似求得差分公式，進一步使用Padé approximation簡化差分公式並且維持精度。

https://en.wikipedia.org/wiki/Compact_finite_difference
https://people.bath.ac.uk/ensdasr/COMPACT/dasr.compact.pdf

modified wave number：針對振盪方程式，針對一種數值方案，修改其係數，使得零一致性變成絕對一致性。

differential equation (oscillation equation):

d           
—— f(x) = 𝑖 k f(x)     k is a parameter (wave number)
dx

symbolic solution:

f(x) = exp(𝑖kx)

numerical scheme (leapfrog method):

f[i+1] - f[i-1]
——————————————— = 𝑖 k f[i]
      2 Δx

substitute LHS with symbolic solution:

f[i+1] - f[i-1]   exp(𝑖k(i+1)Δx) - exp(𝑖k(i-1)Δx)
——————————————— = ——————————————————————————————
      2 Δx                     2 Δx

                  exp(𝑖kΔx) - exp(-𝑖kΔx)
                = ————————————————————— exp(𝑖kiΔx)
                           2 Δx

                  exp(𝑖kΔx) - exp(-𝑖kΔx)
                = ————————————————————— f[i]
                           2 Δx

                    sin(kΔx)
                = 𝑖 ———————— f[i]
                       Δx

modified wave number:

     sin(kΔx)
k' = ————————
        Δx

modified numerical scheme (leapfrog method):

f[i+1] - f[i-1]
——————————————— = 𝑖 k' f[i]
      2 Δx

https://www.nas.nasa.gov/assets/nas/pdf/ams/2018/introtocfd/Intro2CFD_Lecture2_Lecture3_Pulliam_Chap3_Modk.pdf

implicit method / explicit method

決定時刻。

(1) implicit method：右式取新值。
(2) explicit method：右式取舊值。
(3) implicit–explicit method：有些取新值、有些取舊值。
(4) Crank–Nicolson method：新值加舊值然後除以二。

implicit method。取新值。
     d/dt f(t) = f(t) + 1
---> (f[n+1] - f[n]) / Δt = f[n+1] + 1
                            ^^^^^^
explicit method。取舊值。
     d/dt f(t) = f(t) + 1
---> (f[n+1] - f[n]) / Δt = f[n] + 1
                            ^^^^
implicit–explicit method：有些取新值、有些取舊值。
     d/dt f(t) = f(t) + g(t) + 1
---> (f[n+1] - f[n]) / Δt = f[n+1] + g[n] + 1
                            ^^^^^^   ^^^^
Crank–Nicolson method。新值加舊值然後除以二。
     d/dt f(t) = f(t) + 1
---> (f[n+1] - f[n]) / Δt = (f[n+1] + f[n]) / 2 + 1
                            ^^^^^^^^^^^^^^^^^^^

微分方程式，只對同一個變數微分（常微分方程式ODE）的情況下，無法區分差分方法、隱顯方法。古人採用差分方法來命名。

for ordinary differential equation,
forward method     <=> implicit method
backward method    <=> explicit method
trapezoidal method <=> Crank–Nicolson method

numerical iteration

先前介紹了三個項目、兩個技巧。現在追加兩個要素：

(1) matrix solver
(2) domain decomposition method

matrix solver

矩陣求解器。調整演算法，減少計算時間。

0D space：不需要矩陣
1D space：
nD space：

矩陣求解演算法，稱作矩陣求解器。根據矩陣特性，採用適合的矩陣求解器。

一、零維空間動態系統，就一純量，哪來矩陣。要嘛手工移項推導，要嘛函數求解。經典演算法是二分法、牛頓法、割線法。

二、一維空間動態系統，形成雙對角矩陣、三對角矩陣、常對角矩陣。矩陣求解擁有高速演算法。然而本站尚未整理，沒辦法為大家介紹。自己看著辦吧。

三、多維空間動態系統，無法形成上述矩陣。大家直接視作大型稀疏矩陣。經典演算法是鬆弛法、共軛梯度法。經典加速技巧是多格法、預條件法。不適合使用高斯消去法、LU分解。

補充說明一下。循環矩陣，擁有更快的演算法，但是實務上鮮少遇到循環邊界條件。對稱矩陣，擁有更快的演算法，但是遞迴公式係數未必對稱。

domain decomposition method（Schwarz method）

矩陣求解技巧。調整矩陣尺寸，減少疊代次數，加快收斂速度，減少計算時間。

domain decomposition method：領域分解法。矩陣求解切成多個區塊輪流計算。
multigrid method：多格法。矩陣求解採用scaling method。

https://scicomp.stackexchange.com/questions/540/

numerical error

coupling：微分方程組，多個未知函數，取樣地點重合／錯開。

collocated mesh：重合網格。每個未知函數，取樣地點重合。
staggered mesh：錯開網格。有些在網格中央，其餘在網格邊緣。

discretization：固定取樣、動態取樣、固定擇鄰、動態擇鄰。

fixed mesh refinement：固定網格細緻化。事先指定空間各處網格精細程度。
adaptive mesh refinement：自適應網格細緻化。各處根據誤差調整網格精細程度。

error suppression：調整數值解，強行抑制數值誤差。

smooth transition：不連續處強行平滑化，避免Runge phenomenon。
artificial smoothing：整體強行平滑化，避免不穩定。

property preservation：調整數值解，強行維持數學性質。

lake at rest：零／負數強行改成machine epsilon，避免歸零，強行保正。

coupling

(1) forward-backward method
(2) alternating-direction implicit method

smooth transition

不可微函數，實施高次多項式內插、高次多項式近似，導致Runge phenomenon：不可微處，內插函數抖動、近似函數抖動。

遞迴公式是泰勒近似，導致Runge phenomenon，形成disturbance。

當初始條件、邊界條件是不可微函數（例如對稱邊界條件），導致符號解是不可微函數，導致數值解變化劇烈（數值解無法定義連續、可微），導致Runge phenomenon。

波浪撞牆反彈demo

當初始條件、邊界條件是不可微函數，事先強行實施「smooth approximation」，避免Runge phenomenon。

artificial smoothing

概念等同於訊號處理的平滑化、離散微分幾何的平滑化。

平滑化可以強行降低數值解範數，抑制disturbance。平滑化也可以強行降低放大因子，滿足穩定性。

大家習慣使用Laplacian smoothing：追加二次微分項。

temporal smoothing:

f̄[n][i] - f[n][i]     f[n+1][i] - 2 f[n][i] + f[n-1][i]
————————————————— = α —————————————————————————————————
        Δx                          (Δt)²

spatial smoothing:

f̄[n][i] - f[n][i]     f[n][i+1] - 2 f[n][i] + f[n][i-1]
————————————————— = α —————————————————————————————————
        Δt                          (Δx)²

where α is smoothing factor (0 ≤ α ≤ 1)

有些人把分母Δx和(Δt)²併入平滑係數α。有些人直接拋棄。

temporal smoothing:

f̄[n][i] = f[n][i] + α (f[n+1][i] - 2 f[n][i] + f[n-1][i])

spatial smoothing:

f̄[n][i] = f[n][i] + α (f[n][i+1] - 2 f[n][i] + f[n][i-1])

where α is smoothing factor (0 ≤ α ≤ 1)

時間平滑化，經典數值方案是Robert–Asselin filter。

Euler method:

f[n+1][i] = f[n][i] + F(f[n][i], t, x) ⋅ Δt

forward time method:

f[n+1][i] = f[n][i] + F(f[n][i], t, x) ⋅ Δt

leapfrog method (central time method):

f[n+1][i] = f[n-1][i] + F(f[n][i], t, x) ⋅ 2Δt

Robert–Asselin filter:

1. f[n+1][i] = f̄[n-1][i] + F(f[n][i], t, x) ⋅ 2Δt
2. f̄[n][i] = f[n][i] + α (f[n+1][i] - 2 f[n][i] + f̄[n-1][i])

空間平滑化，經典數值方案是Lax–Friedrichs method。

advection equation:

d             d
—— f(t,x) + k —— f(t,x) = 0     k is a parameter
dt            dx

forward-time central-space method (FTCS):

f[n+1][i] = f[n][i] - ½ C (f[n][i+1] - f[n][i-1])
where C = k Δt / Δx is CFL number without absolute value

Lax–Friedrichs method:

1. f[n+1][i] = f̄[n][i] - ½ C (f̄[n][i+1] - f̄[n][i-1])
2. f̄[n+1][i] = f[n+1][i] + ½ (f̄[n][i+1] - 2 f̄[n][i] + f̄[n][i-1])

Lax–Friedrichs method:

f[n+1][i] = f[n][i] - ½ C (f[n][i+1] - f[n][i-1])
          + ½ (f[n][i+1] - 2 f[n][i] + f[n][i-1])

Lax–Friedrichs method:

f[n+1][i] = ½ (f[n][i+1] + f[n][i-1])
          - ½ C (f[n][i+1] - f[n][i-1])

補充說明一下。空間追加二次微分項，沒有一致公認的名稱：artificial/numerical接damping/dissipation/diffusivity/viscosity/resistivity。牛頓流體當中，這些物理現象都是對空間二次微分，即是Laplacian。

時間平滑化已經進化成多步法和多階段法。空間平滑化已經進化為有限元素法。上述兩個數值方案，曾經是豐功偉業，現在卻是陳腔濫調。就當作是學個想法吧。

最後放一篇Lax的文章。

lake at rest

a scheme has C-property if it preserves the steady state
https://hal.science/hal-03893632/file/ms.pdf

numerical convergence

CFL condition：鳩佔鵲巢的假貨偽物。

characteristics / characteristic curve
Courant–Friedrichs–Lewy condition
Courant number 
strong stability preserving temporal discretization

characteristics / characteristic curve

偏微分方程的b² - 4ac ⪋ 0分類法。

Courant–Friedrichs–Lewy condition

專著《Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations》。

一旦微分方程式出現一種特徵曲線走向，那麼遞迴公式也必須出現相同的特徵曲線走向。如此一來，符號解和數值解才有機會相符。

【待補圖片】

另一種比較粗糙的說法：符號解的domain of dependence小於等於數值解的domain of dependence。

採用domain of dependence，不採用region of influence。特徵曲線收斂於一點，得以界定走向範圍。

region of influence：去到哪裡。
domain of dependence：來自哪裡。

【待補圖片】

當遞迴公式恰好形成斜截式，特徵曲線走向恰是斜截式的斜率。

以常係數平流方程式為例。CFL條件是取樣間距Δt與Δx、倍率a，三者構成的不等式。

characteristic curves全是直線，斜率1 / a。數值解必須遵循直線，斜率Δt / Δx。因此a Δt / Δx = 1。

constant-coefficient advection equation:
df     df
—— + c —— = 0
dt     dx

explicit backward Euler method:
                      c Δt
f[n+1][i] = f[n][i] + ———— (f[n][i] - f[n][i-1])
                       Δx

CFL condition for constant-coefficient advection equation:
c Δt
———— = 1
 Δx

以平流方程式為例。CFL條件修改成數值模擬過程當中的倍率的極值。硬撐過去。

characteristic curves走向不一。找到斜率上限與斜率下限，訂立斜率範圍。因此|vₘₐₓ| Δt / Δx ≤ 1。

advection equation:
df     df
—— + v —— = 0
dt     dx

explicit backward Euler method:
                      v[n][i] Δt
f[n+1][i] = f[n][i] + —————————— (f[n][i] - f[n][i-1])
                          Δx

CFL condition for advection equation:
|vₘₐₓ| Δt
————————— ≤ 1     where vₘₐₓ is the extremum value of v
    Δx            in numerical simulation

CFL條件可以視作：移動距離v Δt不得超過取樣間距Δx，避免微分運算離散化抓錯函數值。

微分運算離散化必須抓到相鄰函數值。當移動距離太大，則微分運算離散化誤抓到從遠處跨越鄰居而來的函數值。

【待補圖片】

CFL condition ≠ numerical stability

CFL條件跟穩定性是兩回事！既不充分也不必要！

當運氣很好，遞迴公式萬中選一，CFL條件才會碰巧是穩定性。歷史上首次撞大運的人，就是CFL這三人。

綜觀歷史，先有CFL條件，才有數值分析基本定理。古人先入為主認為CFL條件就是穩定性。即便後來發現數值分析基本定理，大家仍未察覺自身認知有誤。當今所有文獻仍將CFL條件當作敘事主角，成為歷史共業。

雖然CFL條件跟穩定性是兩回事，但是兩者可能有某種關聯。等你發表論文。

順帶一提，針對特定的微分方程式，我們可以利用收縮穩定性，強行定義全新的CFL條件。然而沒有意義。

一種微分方程式有各種遞迴公式，各種遞迴公式有各種收縮穩定性條件，從中找到數學式子最漂亮那一個作為CFL條件。

heat equation       2 |vₘₐₓ| Δt / (Δx)² ≤ 1
wave equation       |vₘₐₓ| Δt / Δx ≤ 1
advection equation  |vₘₐₓ| Δt / Δx ≤ 1

https://people.math.ethz.ch/~grsam/SS21/NAII/resources/slides/ODE-Lecture6.pdf

順帶一提，針對特定的遞迴公式，我們可以將收縮穩定性條件，強行改寫成原始的CFL條件的模樣。然而沒有意義。

一種微分方程式有各種遞迴公式，各種遞迴公式有各種收縮穩定性條件，得到不同的CFL條件右式數值。

CFL condition for numerical scheme of
constant-coefficient advection equation:
時間Euler、空間後向差分（一階精度）　|c| Δt / Δx ≤ 1
時間Euler、空間中央差分（四階精度）　|c| Δt / Δx ≤ 0.728...
時間RK4  、空間中央差分（二階精度）　|c| Δt / Δx ≤ 2.83...

https://scicomp.stackexchange.com/questions/25398/
https://scicomp.stackexchange.com/questions/2927/
https://books.google.com.tw/books?id=npxiU81KY-UC&pg=PA114
https://books.google.com.tw/books?id=npxiU81KY-UC&pg=PA68
https://acfdlab.miami.edu/publications/AIAA2003_0248.pdf
https://www.researchgate.net/publication/267467969

Courant number

CFL條件的左式（遞迴公式的斜截式的斜率），代入v、Δt、Δx的實際數值，求得左式數值，此數值稱作「Courant number」或「CFL number」。

heat equation       C = 2 |vₘₐₓ| Δt / (Δx)²
wave equation       C = |vₘₐₓ| Δt / Δx
advection equation  C = |vₘₐₓ| Δt / Δx

就我所知，Courant number有兩個用途。

(1) Courant number as stability condition
(2) Courant number as amplification factor

針對一個數值方案，我們想要找出確切的穩定性條件。

很不幸地，大多數的數值方案，難以手工推導不等式右側數值。大家只好使用試誤法，枚舉不等式左側數值。

所幸的是，試誤法可以採用binary search，減少試誤時間。

一、設計數值方案。
二、設定初始條件、邊界條件。
三、逐步修改Courant number，重新實施數值模擬，直到數值收斂。
四、公布數值方案、初始條件、邊界條件、穩定性條件（Courant number）。

針對一個數值方案，我們想要調整放大因子：每種頻率的弦波的變化比例。

放大因子內含Courant number。調整Courant number藉以調整放大因子。例如讓高頻弦波變化小、低頻弦波變化大，進而抑制數值誤差。Courant number的地位宛如平滑係數。

strong stability preserving temporal discretization

專著《Strong Stability Preserving Runge–Kutta and Multistep Time Discretizations》。

時間離散化可以推廣為多步多階段。

「保強穩定性時間離散化」。如果一個數值方案，時間離散化採用explicit Euler method，一旦滿足CFL condition，並且導致收縮穩定性，那麼該數值方案，時間離散化推廣為隱式顯式多步多階段，其係數是凸組合，一旦滿足CFL condition乘上適當倍率，也會導致收縮穩定性。

實務上會仔細調整時間離散化的係數，使得適當倍率是1，如同普通的CFL condition。

實務上將時間離散化推廣為多階段，形成strong stability preserving Runge–Kutta method (SSPRK)。

multistage method (SSPRK43):

      shift amount
  Δt     k₁  k₂  k₃  k₄  
┌───────────────────────┐
│  0  │  0   0   0   0  │ round 1
│ 1/2 │ 1/2  0   0   0  │ round 2
│  1  │ 1/2 1/2  0   0  │ round 3
│ 1/2 │ 1/6 1/6 1/6  0  │ round 4
└─────┼─────────────────│
      │ 1/6 1/6 1/6 1/2 │ round weight
      └─────────────────┘

「保強穩定性時間離散化」原理簡單。explicit Euler method只取一處斜率。多步多階段則是取多處斜率的加權平均數。視作多個explicit Euler method各自加權，而多個CFL condition也跟著加權。多個CFL condition不等式聯立，簡化成最緊的那一個不等式。

此處僅證明顯式多階段。可以推廣到隱式顯式多步多階段。

方便起見，取消顯式多階段的最後一個計算步驟：每個回合的計算結果的加權總和。取而代之，將權重併入每個回合的計算公式，稱作Shu–Osher formulation。

differential equation
d/dt f(t,x) = F(f(t,x))

explicit Euler method
f⁽ⁿ⁺¹⁾ = f⁽ⁿ⁾ + Δt F(f⁽ⁿ⁾)                  skip subscripts

contractive stability of linear scheme
‖f⁽ⁿ⁺¹⁾‖ = ‖f⁽ⁿ⁾ + Δt F(f⁽ⁿ⁾)‖ ≤ ‖f⁽ⁿ⁾‖

explicit multistage method (Shu–Osher formulation)
⎧ k₀ = f⁽ⁿ⁾
⎪ :
⎨ kᵢ =  sum  { aᵢⱼ kⱼ + Δt bᵢⱼ F(kⱼ) }
⎪ :   j=1⋯i-1
⎪ :
⎩ f⁽ⁿ⁺¹⁾ = kₘ                               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 explicit Euler method reaches contractive stability
  under CFL condition Δt ≤ Δt_CFL,
  then explicit multistage method reaches contractive stability
  under CFL condition Δt ≤ min (aᵢⱼ/bᵢⱼ) Δt_CFL
                           ⁱʲ
  when a and b are non-negative
  and a satisfies convex combination property.

proof:

  ‖kᵢ‖
= ‖sum { aᵢⱼ kⱼ + Δt bᵢⱼ F(kⱼ) }‖             ith stage
≤ ‖aᵢⱼ‖ ‖sum { kⱼ + Δt (bᵢⱼ/aᵢⱼ) F(kⱼ) }‖     三角不等式
=  aᵢⱼ  ‖sum { kⱼ + Δt (bᵢⱼ/aᵢⱼ) F(kⱼ) }‖     aᵢⱼ ≥ 0
≤  aᵢⱼ   sum ‖ kⱼ + Δt (bᵢⱼ/aᵢⱼ) F(kⱼ) ‖      正負抵銷
≤  aᵢⱼ   sum ‖kⱼ‖                             CFL條件，原因如後。
≤        max ‖kⱼ‖                             凸組合

⎧ ‖k₁‖ ≤ max(‖k₀‖)                    = ‖k₀‖
⎨ ‖k₂‖ ≤ max(‖k₀‖, ‖k₁‖)              ≤ ‖k₀‖
⎪  :          :                          :
⎩ ‖kₘ‖ ≤ max(‖k₀‖, ‖k₁‖, ..., ‖kₘ₋₁‖) ≤ ‖k₀‖

‖f⁽ⁿ⁺¹⁾‖ = ‖kₘ‖ ≤ ‖k₀‖ = ‖f⁽ⁿ⁾‖

for each j,
explicit Euler method reaches contractive stability
‖kⱼ + Δt (bᵢⱼ/aᵢⱼ) F(kⱼ)‖ ≤ ‖kⱼ‖
under CFL condition Δt (bᵢⱼ/aᵢⱼ) ≤ Δt_CFL

「保強穩定性」名稱很雲。此概念既非強穩定性、亦非保Ｘ變換。此概念是將時間離散化從一階泰勒級數推廣成隱式顯式多步多階段，並且考慮CFL condition。

「保強穩定性」內容很雲。前提是一件未必成立的事情。CFL condition不一定導致收縮穩定性。

因此，實務上需要確認前提是否成立。穩定性分析，確認收縮穩定性等價於CFL condition。一旦前提成立，SSPRK所向無敵。

differential equation - dynamical system🚧

differential equation

本章討論微分方程式、遞迴公式的數學性質。

(1) property of differential equation
(2) property of dynamical system

一、微分方程式的數學性質。
二、動態系統的數學性質。

property of differential equation

微分方程式分為兩類。

(1) linear differential equation
    一次。函數微分視作變數，形成一次函數。
(2) nonlinear differential equation
    非一次。函數微分視作變數，形成非一次函數。

微分方程式的常見細類。

(a) homogeneous equation
    齊次。函數微分視作變數，沒有常數函數項。
    導致符號解可以是零函數。
(b) quasilinear equation
    擬一次。函數微分視作變數，係數不是常數函數而是泛函數，
    泛函數的輸入變數不含最高次微分。

符號解分為兩類。

(1) single-valued solution
    單值解。符號解是函數。一個位置擁有一個函數值。
(2) multivalued solution (set-valued solution)
    多值解。也可以不是函數。一個位置擁有多個函數值。

符號解的常見細類。

(a) analytic solution
    解析解。符號解是解析函數。
(b) weak solution
    弱解。也可以是不可微函數。

微分方程式暨符號解的數學性質，遞迴公式暨數值解的數學性質，兩邊必須一致。下文只講前者，講一個就等於兩個都講了。

順帶一提，數值解無法定義連續、可微。數值解不存在所謂的解析解、弱解。

【名稱尚待查詢】

線性微分方程式，初始條件是單值，符號解亦是單值。

非線性微分方程式，初始條件是單值，符號解可能是多值。從某個時刻開始，單值變多值，函數變流形。

differentiable solution

「可微解」。微分方程式的符號解，必須是可微函數。

詳細來說，內含N次微分的微分方程式的符號解，必須是N次可微函數。顯而易見的廢話。

宏觀來說，符號解可以是下述函數。然而大家只討論解析函數。

N次可微函數：處處N次可微。
N次分段函數：銜接處N次可微，其餘處更高次可微。
平滑函數：處處無窮次可微。
解析函數：處處無窮次可微，且處處導數是泰勒級數。
　　　　　（處處導數是無窮多項式函數，其輸入移位。）

analytic solution

「解析解」。微分方程式的符號解，設定為解析函數。

數學家從微積分基本定理出發，利用指數函數、三角函數，求得微分方程式的符號解。利用解析函數，自然而然得到解析解。

數學家從數值分析基本定理出發，利用微分運算的泰勒近似，求得微分方程式的數值解。數值解是近似解。全域誤差是近似誤差。而截斷誤差（一致性）影響了全域誤差（收斂性）。想讓截斷誤差趨近零，一種方式是採用解析解。先前章節已經提及。

weak solution

「弱解」。微分方程式的符號解，推廣為不可微函數。

現實世界流行不可微函數。初始條件是不可微函數，符號解亦是不可微函數。此時便需要弱解。

首先介紹一個轉換手法。目前沒有正式學術名稱。

一個函數L(t,x)，替它找一個平滑函數ϕ(t,x)，相乘、積分，重新包裝成一個可微函數∫∫ L(t,x) ϕ(t,x) dt dx。

接著，窮舉各種平滑函數，得到各種可微函數。

不同的函數L(t,x)，實施轉換，得到不同的結果。結果必定獨一無二，畢竟窮舉了各種平滑函數。

           ⎧ ⌠+∞ ⌠+∞                     │                  ⎫
L(t,x) <-> ⎨ ⎮   ⎮   L(t,x) ϕ(t,x) dt dx │ ϕ(t,x) is smooth ⎬
           ⎩ ⌡-∞ ⌡-∞                     │                  ⎭

舉例來說，下述微分方程式：

df(t,x)   df(t,x)
——————— + ——————— = 0
  dt        dx

整個式子經過包裝、展開、Fubini定理、提項，得到新式。

⌠+∞ ⌠+∞ ⎛ df(t,x)   df(t,x) ⎞
⎮   ⎮   ⎜ ——————— + ——————— ⎟ ϕ(t,x) dt dx = 0
⌡-∞ ⌡-∞ ⎝   dt        dx    ⎠

⌠+∞ ⌠+∞ df(t,x)                ⌠+∞ ⌠+∞ df(t,x)
⎮   ⎮   ——————— ϕ(t,x) dt dx + ⎮   ⎮   ——————— ϕ(t,x) dt dx = 0
⌡-∞ ⌡-∞   dt                   ⌡-∞ ⌡-∞   dx

⌠+∞ ⌠+∞        dϕ(t,x)         ⌠+∞ ⌠+∞        dϕ(t,x)
⎮   ⎮   f(t,x) ——————— dt dx + ⎮   ⎮   f(t,x) ——————— dt dx = 0
⌡-∞ ⌡-∞          dt            ⌡-∞ ⌡-∞          dx

⌠+∞ ⌠+∞ ⎛        dϕ(t,x)          dϕ(t,x) ⎞
⎮   ⎮   ⎜ f(t,x) ——————— + f(t,x) ——————— ⎟ dt dx = 0
⌡-∞ ⌡-∞ ⎝          dt               dx    ⎠

再舉例來說，下述初始值問題：

⎧ df(t,x)   df(t,x)
⎪ ——————— + ——————— = 0   where t ≥ 0
⎨   dt        dx
⎪ 
⎩ f(0,x) = f₀(x)

叭啦叭啦叭啦，得到新式：

⌠+∞ ⌠+∞ ⎛        dϕ(t,x)          dϕ(t,x) ⎞
⎮   ⎮   ⎜ f(t,x) ——————— + f(t,x) ——————— ⎟ dt dx
⌡0  ⌡-∞ ⎝          dt               dx    ⎠

  ⌠+∞
+ ⎮   f₀(x) ϕ(0,x) dx = 0
  ⌡-∞

原式稱作strong formulation。新式稱作weak formulation。

窮舉每一個平滑函數，新式一律出現的解，定義為原式的弱解。新式的解f(t,x)可以是不可微函數，原式的弱解可以是不可微函數。

property of dynamical system

引入時間變數與空間變數之後，符號解的常見細類。

(1) diminishing
    遞減。符號解範數每回合變小或不變。
(2) monotonicity preserving
    保單調。符號解每回合保持單調。
(X) total variation diminishing
    總變差遞減。符號解相鄰峰谷高低差總和每回合變小或不變。
(3) local extremum diminishing
    局部極值遞減。符號解每回合極值數量不增加、極值大小不變廣。
(4) positivity preserving
    保正。符號解每回合保持非負數（保持正負號）。

微分方程式暨符號解的數學性質，遞迴公式暨數值解的數學性質，兩邊必須一致。下文只講後者，講一個就等於兩個都講了。

這些詞彙可以是動名詞、也可以是形容詞。形容詞可以放在be動詞之後、也可以放在名詞之前。放在名詞之前的情況下，常見用法是接上scheme，例如monotonicity-preserving scheme，表示該數值方案具備此數學性質。

diminishing

「遞減」。數值解範數每回合變小或不變。

L_∞-diminishing max |f̂⁽ⁿ⁺¹⁾| ≤ max |f̂⁽ⁿ⁾|
L₁-diminishing ‖f̂⁽ⁿ⁺¹⁾‖₁ ≤ ‖f̂⁽ⁿ⁾‖₁
Lₚ-diminishing ‖f̂⁽ⁿ⁺¹⁾‖ₚ ≤ ‖f̂⁽ⁿ⁾‖ₚ

知名範例是指數衰減微分方程式。微分不動點。

exponential decay differential equation
d
—— f(t) = k f(t)     k is a parameter and k < 0
dt

explicit Euler method
fᵢ⁽ⁿ⁺¹⁾ - fᵢ⁽ⁿ⁾
——————————————— = k fᵢ⁽ⁿ⁾
       Δt

fᵢ⁽ⁿ⁺¹⁾ = (1 + Δt k) fᵢ⁽ⁿ⁾

diminishing if (1 + Δt k) ≤ 1

知名範例是振盪方程式。彈簧振動、琴弦振動、鼓膜振動的符號解，可以寫成複數的實部。儘管實部乍看反覆變大變小，但是複數範數其實保持變小或不變，滿足遞減。【尚待確認】

diminishing數學定理

一、齊次遞迴公式，收縮則遞減。

證明很簡單。當微分方程式／遞迴公式沒有常數項，則符號解／數值解可以是零。零作為減數，那麼收縮導致遞減。

for homogeneous scheme,
contractive => diminishing

二、線性遞迴公式，收縮即是遞減即是指數衰減。

證明很簡單。兩個符號解／數值解相加減，依然形成符號解／數值解，數學性質依然成立。

for linear scheme,
contractive <=> diminishing <=> exponential decay

monotonicity preserving

「保單調」。數值解每回合保持單調。

當遞迴公式不因地點而變，數值解的任意區段皆具備相同數學性質。保單調變成：數值解拆解成單調區間們與局部極值們，單調區間保持單調區間，局部極值保持局部極值。方便起見，單調區間們與局部極值們可以移動地點。

或者簡單一句話：數值解每回合都不會新增局部極值。

知名範例是淺水方程組。【方程組如何定義單調】

知名範例是shallow water equations，數值解由高度和速度組成，深水導致波浪降低加速，淺灘導致波浪升高減速。波浪隨著水深起伏，形狀大致不變，高處恆高、低處恆低。

monotonicity preserving數學定理

Godunov's monotonicity preserving theorem：針對線性遞迴公式，保單調就是正組合。

Godunov是第一位著眼於保單調的人。他也發明了上風方案，確保正組合，進而確保保單調。留待後面章節介紹。

    monotonicity preserving
<=> positive combination (positivity)

    monotonicity preserving
    數值解每回合保持單調

    (fᵢ⁽ⁿ⁾ is locally monotone) => (fᵢ⁽ⁿ⁺¹⁾ is locally monotone)
    fᵢ₋₁⁽ⁿ⁾ ≤ fᵢ⁽ⁿ⁾ ≤ fᵢ₊₁⁽ⁿ⁾ => fᵢ₋₁⁽ⁿ⁺¹⁾ ≤ fᵢ⁽ⁿ⁺¹⁾ ≤ fᵢ₊₁⁽ⁿ⁺¹⁾

<=> positive combination (positivity)
    遞迴公式的係數皆是正值

    ⎰ fᵢ⁽ⁿ⁺¹⁾ = sumⱼ cᵢⱼ fⱼ⁽ⁿ⁾
    ⎱ cᵢⱼ ≥ 0

f is monotone:
min(fᵢ₋₁, fᵢ₊₁) ≤ fᵢ ≤ max(fᵢ₋₁, fᵢ₊₁) for all i
fᵢ₋₁ ≤ fᵢ ≤ fᵢ₊₁ for all i

f is locally monotone:
fᵢ₋₁ ≤ fᵢ ≤ fᵢ₊₁ at i

https://en.wikipedia.org/wiki/Godunov's_theorem
https://books.google.com.tw/books?id=AAqMDwAAQBAJ&pg=PA257

Godunov's order barrier theorem：針對線性遞迴公式，保單調的精度至多一階。進一步推理，正組合、單調、保正，精度都是至多一階。

階數屏障定理有各種版本。Godunov首開先例。

monotonicity preserving scheme
is at most first-order accurate.
[Godunov 1954]

https://en.wikipedia.org/wiki/Godunov's_theorem
https://hal.science/hal-01620642
https://www.jstor.org/stable/2008046

local extremum diminishing（range diminishing）

專著《Computational Aerodynamics》。

「局部極值遞減」。數值解同時滿足下述三點。

(1) nonincreasing local maximum 數值解局部極大值減少或不變
(2) nondecreasing local minimum 數值解局部極小值增加或不變
(3) no new extremum             數值解局部極值數量減少或不變

當遞迴公式不因地點而變，局部極值遞減導致保單調且遞減。

反方向不一定正確。例如全域收縮、局部放大。

   local extremum diminishing (LED)
=> monotonicity preserving & diminishing

針對線性遞迴公式，變數是自己和所有鄰居，係數是凸組合。

for linear scheme,
a LED scheme is a neighborhood scheme.

一維空間是三點。二維空間是十字五點。

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運算。

1D heat equation
∂           ∂²        
—— f(t,x) = ——— f(t,x)
∂t          ∂x²       

explicit Euler method
fᵢ⁽ⁿ⁺¹⁾ - fᵢ⁽ⁿ⁾   (fᵢ₋₁⁽ⁿ⁾ - fᵢ⁽ⁿ⁾) + (fᵢ₊₁⁽ⁿ⁾ - fᵢ⁽ⁿ⁾)
——————————————— = —————————————————————————————————————
       Δt                          2 Δx

explicit Euler method
fᵢ⁽ⁿ⁺¹⁾ = (1 - λ) fᵢ⁽ⁿ⁾ + ½ λ fᵢ₋₁⁽ⁿ⁾ + ½ λ fᵢ₊₁⁽ⁿ⁾
where λ = Δt/Δx

local extremum diminishing數學定理

針對線性遞迴公式，LED的衍生性質。

專著《Computational Gasdynamics》。

   local extremum diminishing (LED)
   遞迴公式的變數是自己和所有鄰居、係數是凸組合
   一、數值解局部極大值每回合減少或不變
   二、數值解局部極小值每回合增加或不變
   三、數值解局部極值數量每回合變少或不變

   ⎧ fᵢ⁽ⁿ⁺¹⁾ = sumⱼ cᵢⱼ fⱼ⁽ⁿ⁾
   ⎨ sumⱼ cᵢⱼ = 1
   ⎪ cᵢⱼ ≥ 0
   ⎩ cᵢⱼ = 0 , if j ∉ neighbor(i)
   
=> convex combination (convexity)
   遞迴公式的係數是凸組合

   ⎧ 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ⱼₖ⁽ⁿ⁾ ≤ gⱼₖ⁽ⁿ⁾ => fᵢ⁽ⁿ⁺¹⁾ ≤ gᵢ⁽ⁿ⁺¹⁾  for all 1 ≤ k ≤ m

   f̂⁽ⁿ⁾ ≤ ĝ⁽ⁿ⁾ => f̂⁽ⁿ⁺¹⁾ ≤ ĝ⁽ⁿ⁺¹⁾

[Harten 1983] [Osher 1985] [Jameson 1995]

LED抑制振盪、漸趨平緩。此性質稱作maximum principle。

振盪程度可以具體表示成數值大小L_∞-norm與total variation。此性質稱作L_∞-diminishing與total variation diminishing。

   order preserving (monotone) & diminishing

=> maximum principle
   數值解下回合數值介於遞迴公式的變數極值之間
   證明手法：構造另一個數值解，恰是變數極值。

   min {fⱼ⁽ⁿ⁾ | cᵢⱼ > 0} ≤ fᵢ⁽ⁿ⁺¹⁾ ≤ max {fⱼ⁽ⁿ⁾ | cᵢⱼ > 0}

=> L_∞-diminishing
   數值解每回合振幅減少或不變

   ‖f̂⁽ⁿ⁾‖_∞ ≤ ‖f̂⁽ⁿ⁺¹⁾‖_∞

   where ‖f‖_∞ = maxᵢ |fᵢ|

   order preserving (monotone) & diminishing & neighborhood (3-point)

=> maximum principle & neighborhood (3-point)
   數值解下回合數值介於本回合自己與鄰居的極值之間
   （當邊界條件是預定義數值，數值解最終回合介於邊界極值之間。）

   min {fⱼ⁽ⁿ⁾ | j ∈ neighbor(i)} ≤ fᵢ⁽ⁿ⁺¹⁾ ≤ max {fⱼ⁽ⁿ⁾ | j ∈ neighbor(i)}

=> total variation diminishing (TVD)
   數值解每回合總差值減少或不變

   ‖f̂⁽ⁿ⁾‖ᴛᴠ ≤ ‖f̂⁽ⁿ⁺¹⁾‖ᴛᴠ

   where ‖f‖ᴛᴠ = sumᵢ |fᵢ₊₁ - fᵢ|

Harten's TVD theorem：當遞迴公式不因地點而變，數值解邊界維持不變，數值解內部滿足total variation diminishing，那麼保單調。應該說是LED。【尚待確認】

   Dirichlet boundary condition
 & total variation diminishing (TVD)
=> monotonicity preserving
[Harten 1983]

我認為大家都搞錯這個定理。傳世版本如下：

   L₁-contraction
=> total variation diminishing (TVD)
=> monotonicity preserving   ✘

　 證明手法：數值解相減，令被減數是減數往右移位，湊出相鄰差。

第一個衍生是開放邊界條件（數值解左右延伸至無限遠）。第二個衍生是預定義邊界條件（數值解兩端固定）。邊界條件必須統一，才能一路衍生。

   L₁-diminishing
=> L₁-diminishing & total variation diminishing (TVD)
=> monotonicity preserving

, under Dirichlet boundary condition at infinity

positivity preserving

「保正」。數值解每回合保持非負數（保持正負號）。

微分方程組，全部變數需要保正，可以採用LED。針對線性遞迴公式，LED則maximum principle則保正。

maximum principle => positivity preserving

微分方程組，部分變數需要保正，例如質量保正、動量不保正，事情相當複雜。我還沒有學會。

Chi-Wang Shu 舒其望
https://apps.dtic.mil/sti/tr/pdf/ADA273991.pdf

Xiangxiong Zhang 张翔雄
https://www.math.purdue.edu/~zhan1966/research/paper/review.pdf
https://www.math.purdue.edu/~zhan1966/research/paper/euler.pdf
https://www.math.purdue.edu/~zhan1966/research/paper/fd-weno.pdf
https://www.math.purdue.edu/~zhan1966/research/index.html

Kailiang Wu 吴开亮
https://arxiv.org/pdf/2410.05173
https://sites.google.com/site/klwuhomepage/

differential equation - conservation law🚧

conservation law（conservative differential equation）

本文介紹動態系統的其中一種特例：守恆律。

物理學之守恆定律：物理量總和固定。

例如質量守恆、動量守恆、能量守恆。不生不滅，相加為常數。

m₁ + m₂ = mₜₒₜₐₗ             conservation of mass
m₁v₁ + m₂v₂ = pₜₒₜₐₗ         conservation of momentum
½m₁v₁² + ½m₂v₂² = Eₜₒₜₐₗ     conservation of energy

進一步考慮物理量的時間變化。

例如力平衡。不生不滅，相加為零。

d/dt mv = f₁ + f₂ = 0        static equilibrium of force
d/dt mv + f₁ + f₂ = 0        dynamic equilibrium of force

進一步考慮物理量的空間變化。

例如質量傳輸、動量傳輸、能量傳輸。不生不滅，相加為零。

d       d       d       d
—— ρV + —— ρV + —— ρV + —— ρV = 0           transport of mass
dt      dx      dy      dz

數學之守恆律：物理量的時間變化等於空間變化。

教科書習慣將未知函數寫成u(t,x,y,z)，簡寫為u。

d      d      d      d
—— u + —— u + —— u + —— u = 0               conservation law
dt     dx     dy     dz

進一步讓空間變化乘上已知倍率a。時間變化不必乘上倍率。倍率取決於未知函數的輸出數值大小，因此寫成函數a(u)。

d            d            d            d
—— u + a₁(u) —— u + a₂(u) —— u + a₃(u) —— u = 0
dt           dx           dy           dz
                                            conservation law

進一步利用微分連鎖律，將倍率a改寫成變換f。形成守恆律的標準式。

df(u)   df(u) du         du         du
————— = ————— —— = f′(u) —— := a(u) ——      chain rule
 dx      du   dx         dx         dx

d      d          d          d
—— u + —— f₁(u) + —— f₂(u) + —— f₃(u) = 0   conservation law
dt     dx         dy         dz

順帶一提，未知函數可以推廣成多變量。倍率a(u)隨之推廣成矩陣A(u)。數學家習慣討論其中一種特例：A(u)可對角化、而且特徵值是實數。

主要原因是可對角化矩陣容易推導數學公式。透過特徵分解和座標轉換，將一個多變量函數化為多個單變量函數（一道方程組化為多道方程式），其變量不相依，得以各自討論。工程數學領域，用特徵分解對付多自由度系統，是名聞遐邇的手法。

古人將這種特例稱作hyperbolic conservation law，沿襲至今。注意到，此處的hyperbolic，偏微分方程的b² - 4ac ⪋ 0分類法的hyperbolic，兩者意義不同。此處的守恆律，在b² - 4ac ⪋ 0分類法當中是屬於parabolic，畢竟沒有二階偏微分。

取名hyperbolic是因為古人沒有考慮清楚。成為歷史共業。

scalar conservation law in one-dimensional space

專著《Finite Volume Methods for Hyperbolic Problems》。

專著《Riemann Problems and Jupyter Solutions》。

本文只談最簡單的情況：一維純量守恆律。

d/dt u + d/dx f(u) = 0

(1) one-dimensional space (1D):
    the input of u are two values t and x.
(2) scalar-valued (scalar):
    the output of u is a value.

一維：雙變數。函數輸入是時間變數t、空間變數x。
純量：單變量。函數輸出是數值、而不是向量。

以下列出一維純量守恆律的其中幾個經典範例。

constant-coefficient                 常係數平流方程式
advection equation                   q: physical quantity
u ≔ q , f(u) ≔ cq                    c: velocity (constant)

inviscid Burgers equation            無黏性漢堡堡方程式
u ≔ q , f(u) ≔ ½q²                   q: physical quantity

traffic-flow equation                交通流量方程式
u ≔ ρ , f(u) ≔ ρ(1−ρ)                ρ: density of cars

continuity equation                  連續方程式
u ≔ q , f(u) ≔ qv                    q: physical quantity
                                     v: velocity (known)

shallow water equations              淺水方程組
⎰ u₁ ≔ h  , f₁(u₁) ≔ hv              h: water depth
⎱ u₂ ≔ hv , f₂(u₂) ≔ hv²+½gh²        v: velocity
                                     g: gravity (constant)

Euler equations of gas dynamics      氣體動力學之歐拉方程組
⎧ u₁ ≔ ρ  , f₁(u₁) ≔ ρv              ρ: density
⎨ u₂ ≔ ρv , f₂(u₂) ≔ ρv²+p           v: velocity
⎩ u₃ ≔ E  , f₃(u₃) ≔ (E+p)v          p: pressure
                                     E: energy (E = ρe + ½ρv²)
                                     e: internal energy
                                        per unit mass

順帶一提，最後兩個例子是方程組。方程組可以進一步改寫成多變量方程式，恰好形成hyperbolic conservation law。

shallow water equations              淺水方程組
u ≔ ⎡ h  ⎤ , f(u) ≔ ⎡ hv       ⎤     h: water depth
    ⎣ hv ⎦          ⎣ hv²+½gh² ⎦     v: velocity
                                     g: gravity (known)

Euler equations of gas dynamics      氣體動力學之歐拉方程組
    ⎡ ρ  ⎤          ⎡ ρv     ⎤       ρ: density
u ≔ ⎢ ρv ⎥ , f(u) ≔ ⎢ ρv²+p  ⎥       v: velocity
    ⎣ E  ⎦          ⎣ (E+p)v ⎦       p: pressure
                                     E: energy (E = ρe + ½ρv²)
                                     e: internal energy
                                        per unit mass

solution of scalar conservation law in one-dimensional space

線性守恆律，初始條件是函數，那麼時時刻刻都是函數，形成單值解。例如constant-coefficient advection equation：

非線性守恆律，途中可能變得不是函數，形成多值解。例如inviscid Burgers equation：

shock wave

「衝擊波」或「震波」是一種物理現象。

例如空氣衝擊波，空氣高速移動、向前推擠、密度陡升。下圖是戰鬥機引起的空氣衝擊波，密度陡升之處形成黑線，黑線隨著戰鬥機前進。

「shock tube」是一種實驗設備，用於觀察封閉長管內部的空氣衝擊波，可以視作一維空間的衝擊波。「Sod shock tube」是知名的實驗模型暨實驗數據，用來檢驗數學模型暨計算結果是否屬實。

例如正在破碎的波浪，海水流動、向前推擠、高度和速度上升、前緣破碎。下圖是風力引起的波浪，高度和速度上升之處正在破碎。

「wave tank」是一種實驗設備，用於觀察封閉水道內部的波浪，可以視作一維空間的波浪。網路上可以找到許多實驗模型暨實驗數據，用來檢驗數學模型暨計算結果是否屬實。

按理來說，衝擊波滿足守恆定律。理所當然，大家運用守恆律描述衝擊波。衝擊波是守恆律的其中一種解，而且是不可微函數。空氣衝擊波是Euler equations of gas dynamics其中一解。正在破碎的波浪是shallow water equations其中一解。

以往發展的微分方程式求解技巧，只能求得解析解，不能求得衝擊波。數學家為此發展新的求解技巧。

shock wave in symbolic solution

專著《Computational Gasdynamics》。

常係數守恆律。

constant-coefficient conservation law
d             d
—— u(t,x) + a —— u(t,x) = 0
dt            dx

initial condition u(0,x) = u₀(x)
solution          u(t,x) = u₀(x - at)
waveform          u(t,x) at certain t
wave speed        a = constant
 rightward        a > 0
 leftward         a < 0
 static           a = 0
wavefront         x - at = constant
characteristics   dx/dt = a
expansion         a u(t,x) ≤ a u(t,y)   for x₀ ≤ x ≤ y ≤ y₀
compression       a u(t,x) ≥ a u(t,y)   for x₀ ≤ x ≤ y ≤ y₀

一般的守恆律。

conservation law
d                     d
—— u(t,x) + a(u(t,x)) —— u(t,x) = 0
dt                    dx

d                      d
—— u(t,x) + f′(u(t,x)) —— u(t,x) = 0     define f′ = a
dt                     dx

d           d
—— u(t,x) + —— f(u(t,x)) = 0             chain rule
dt          dx

initial condition u(0,x) = u₀(x)
solution          u(t,x) = u₀(x - a(u(x,t))⋅t)
waveform          u(t,x) at certain t
wave speed        a(u(t,x))
 rightward        a(u(t,x)) > 0
 leftward         a(u(t,x)) < 0
 static           a(u(t,x)) = 0
wavefront         x - a(u)t = constant
characteristics   dx/dt = a(u)
expansion         f(u(t,x)) ≤ f(u(t,y))   for x₀ ≤ x ≤ y ≤ y₀
compression       f(u(t,x)) ≥ f(u(t,y))   for x₀ ≤ x ≤ y ≤ y₀

符號解可以畫成函數曲線圖（橫軸位置x、縱軸函數值u）、特徵曲線圖（橫軸位置x、縱軸時間t）。

(1) waveforms (function curves)   u(t,x) for each t
(2) characteristic curves         dx/dt = a(u) for each x

不同初始條件，得到不同符號解。大家習慣討論最簡易的情況Riemann problem：初始條件是步進函數（兩段常數函數）。

scalar conservation law in one-dimensional space
d           d
—— u(t,x) + —— f(u(t,x)) = 0
dt          dx

solution
u(t,x)

Riemann initial condition
方便教學的簡易範例。不連續處只有一個，函數值是兩個常數。
u(0,x) = ⎰ uʟ(0)   if x ≤ p(0)
         ⎱ uʀ(0)   if x > p(0)

discontinuous position of solution
p(t)

move speed of discontinuous position
s(t) = p′(t)

function value of discontinuous position of solution
uʟ(t) =  lim  u(t,x)
       x→p(t)⁻
uʀ(t) =  lim  u(t,x)
       x→p(t)⁺

compression
f(uʟ) ≥ f(uʀ)     at certain t

Rankine–Hugoniot jump condition
左側溢出等於右側溢出。物理量守恆。
f(uʟ) - s uʟ = f(uʀ) - s uʀ     for all t

Rankine–Hugoniot jump condition
改寫數學式子，宛如平緩函數。
f(uʟ) - f(uʀ)
————————————— = s     for all t
   uʟ - uʀ

Oleinik admissibility condition
一種特殊的條件，得到唯一一種特殊的符號解：純粹的衝擊波。
左側一律大於（小於）右側斷點收縮率。總是往右推進。
f(uʟ) - f(u*)       f(u*) - f(uʀ)     for all t
————————————— ≥ s ≥ —————————————     for all u*
   uʟ - u*             u* - uʀ        that uʟ > u* > uʀ

Lax admissibility condition
f恰是凸函數，那麼可以簡化數學式子。
f′(uʟ) > s > f′(uʀ)     for all t

符號解的存在性、唯一性。

衝擊波有無限多種造型。大家限定造型：characteristic curves只合併、不分叉。白話解釋是衝擊波力道不會渙散。形成唯一解。

一種初始條件對應一種純粹的衝擊波。

(1) weak solution existence <=>
    Rankine–Hugoniot jump condition

(2) weak solution uniqueness <=>
    Rankine–Hugoniot jump condition &
    Oleinik admissibility condition

(3) weak solution uniqueness when f is convex <=>
    Rankine–Hugoniot jump condition &
    Lax admissibility condition

符號解細分為兩種：衝擊波、稀疏波。

Riemann problem:
(1) Riemann initial condition
(2) Rankine–Hugoniot jump condition
(3) Lax admissibility condition

two types of solutions of Riemann problem
(1) rarefaction wave      f(uʟ) ≤ f(uʀ)
(2) shock wave            f(uʟ) ≥ f(uʀ)

initial condition of Riemann problem
u(0,x) = ⎰ uʟ , if x < 0
         ⎱ uʀ , if x > 0

solution: shock wave
u(t,x) = ⎰ uʟ , if x < st
         ⎱ uʀ , if x > st
where s is speed s = (f(uʀ) - f(uʟ)) / (uʀ - uʟ)

solution: rarefaction wave
u(t,x) = ⎧ uʟ        , if x ≤ f′(uʟ)t
         ⎨ f′⁻¹(x/t) , if f′(uʟ)t ≤ x ≤ f′(uʀ)t
         ⎩ uʀ        , if x ≤ f′(uʀ)t

[Osher 1983]
since waveform is monotonicity preserving
u(x/t) = ⎰ argmin {f(u*) - (x/t)u* : uʟ ≤ u* ≤ uʀ} , if uʟ ≤ uʀ
         ⎱ argmax {f(u*) - (x/t)u* : uʀ ≤ u* ≤ uʟ} , if uʟ > uʀ

shock wave
(1) function curve           u(t,x) for each t
u uʟ           u              u
^ --           ^ ----         ^ ------
|  | uʀ    --> |    |     --> |      |
|  ------      |    ----      |      --
--------->x    --------->x    --------->x
(2) characteristic curves    dx/dt = f′(u) for each x
t                         t
^      ／/|               ^||||/／／／／／
|   ／／/||               ||||/／／／／
|／／／/|||               |||/／／／
--------->x               --------->x
entropy-consistent        entropy-violating
run into shock            emerge from shock
(中線記得對準discontinuity)

rarefaction wave
(1) function curve
u uʟ           u              u
^ --           ^ --           ^ --
|  | uʀ    --> |   \      --> |    ＼
|  ------      |    ----      |       --
--------->x    --------->x    --------->x
(2) characteristic curves
t                         t
^||  /  ／／／            ^||||/／／／／／
||| / ／／／              ||||/／／／／
|||/／／／                |||/／／／
--------->x               --------->x
entropy-consistent        entropy-violating
expansion fan             create from center
(是咧畫啥物潲)

https://www3.nd.edu/~dbalsara/numerical-pde-course/Appendix_LesHouches/LesHouches_Lecture_4_NonLinScalr.pdf
https://web.stanford.edu/class/math220a/handouts/conservation.pdf
https://www.researchgate.net/publication/266329534
https://math.stackexchange.com/questions/2539265/

符號解的波形變化過程細分為三種：消滅、延展、創造。

three types of behavior of waveforms
(1) destroy
(2) expand
(3) create

符號解的波形是monotonicity preserving、L_∞-diminishing、total variation diminishing (TVD)。

waveforms are monotonicity preserving
https://books.google.com.tw/books?id=pcwLAQAAQBAJ&pg=PA67

waveforms are L_∞-diminishing and TVD [Glimm–Lax 1970]
https://books.google.com.tw/books?id=-gHkBwAAQBAJ&pg=PA291  

waveforms satisfy tame oscillation condition
https://bpb-us-e1.wpmucdn.com/sites.psu.edu/dist/d/80666/files/2019/09/clawtut4-Oxford.pdf

一維純量守恆律的方程組，符號解出現新種類。此處省略。

Riemann problem of system of conservation laws in 1D space
has new type discontinuity.
i.e. contact discontinuity, tangential discontinuity.

高維純量守恆律，存在性、唯一性，必須砍掉重練。此處省略。

Riemann problem of conservation law in 2D space is not easy.
see books of Lax and 劉太平.
《Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves》
《Shock Waves》

shock wave in numeric solution

古代人將多值解對切對拼得到衝擊波，令兩半面積相等以便守恆。現代人將衝擊波視作容許解，並且發明容許解的數值方案。

容許解的數值方案，目前通通都是LED。線性遞迴公式，又符合下述三點，等價於LED。

(1) monotonicity preserving
    衝擊波符號解保單調。因此衝擊波數值解也必須保單調。
(2) diminishing
    大家習慣採用收縮穩定性。因此符號解／數值解必須遞減。
(3) neighborhood
    微分離散化總是取相鄰數值。因此遞迴公式總是取相鄰數值。

analytic solution：「解析解」。解是連續平滑函數。解析解不包含衝擊波、稀疏波。

weak solution：「弱解」。解可以是不連續函數。弱解包含了衝擊波、稀疏波、各種奇形怪狀的不連續函數。

admissible solution (entropy solution)：近期稱作「容許解」。早期稱作「熵解」。滿足Oleinik admissibility condition的符號解。

admissible measure-valued solution (entropy process solution)：抱歉我學不會。我選擇放棄。

differential equation - numerical method🚧

numerical method

數值模擬，目前已有六大類方法。

(1) finite difference method：有限差分法。局部近似。取樣。
(2) finite volume method：有限體積法。局部近似。通量。
(3) finite element method：有限元素法。全域近似。引力。
(4) spectral method：頻譜法。全域近似。諧波。
(5) method of characteristics：特徵線法。速度向量場。軌跡。
(6) level set method：等高線法。位勢純量場。前沿。

https://en.wikipedia.org/wiki/Numerical_methods_for_partial_differential_equations

finite difference method

有限差分法。可以求得解析解。有限差分法是一類數值方案的泛稱。先前的動態系統章節，整篇文章都是採用有限差分法。

domain discretization: grid
codomain discretization: function value at point
operator discretization: finite difference

定義域離散化：方格
對應域離散化：單點數值
微分運算離散化：有限差分（後向差分／前向差分／中央差分）

數值方案Lax–Friedrichs method與Lax–Wendroff method。

Lax–Friedrichs method: 1st order accuracy
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (f(uᵢ₊₁)⁽ⁿ⁾ - f(uᵢ₋₁)⁽ⁿ⁾) / 2
uᵢ⁽ⁿ⁺¹⁾ = (uᵢ₊₁⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾) / 2
        - C (f(uᵢ₊₁)⁽ⁿ⁾ - f(uᵢ₋₁)⁽ⁿ⁾) / 2

Lax–Wendroff method: 2nd order accuracy
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 splitting。

flux averaging (flux correction) (flux limiting)
f⁽ⁿ⁾ = wˡᵒʷ fˡᵒʷ + wʰⁱᵍʰ fʰⁱᵍʰ = lerp(fˡᵒʷ, fʰⁱᵍʰ)
     = fˡᵒʷ + w (fʰⁱᵍʰ - fˡᵒʷ) = fˡᵒʷ + diff(fʰⁱᵍʰ, fˡᵒʷ)
where wˡᵒʷ and wʰⁱᵍʰ forms convex combination
      0 ≤ w ≤ 1 forms convex combination

flux splitting
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ + C (f⁺ - f⁻)
where f⁺ is increasing and f⁻ is decreasing

upwind sheme根據流向決定差分地點，差分地點總是位於上風處。upwind scheme用來確保正組合，進而確保保單調。upwind scheme可以視作一種特殊的flux splitting。

upwind scheme: conditional function
uᵢ⁽ⁿ⁺¹⁾ = ⎰ uᵢ⁽ⁿ⁾ - C a (uᵢ⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾) , if a > 0
          ⎱ uᵢ⁽ⁿ⁾ - C a (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) , if a < 0

upwind scheme: flux splitting
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C a⁺ (uᵢ⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾)
                - C a⁻ (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾)
where a⁺ = max(a,0)
      a⁻ = min(a,0)
      a⁺ - a⁻ = |a|

finite volume method

有限體積法。可以求得守恆律的弱解。而且保證物理量守恆。

domain discretization: grid/mesh
codomain discretization: average function value
operator discretization: finite difference

定義域離散化：方格／網格
對應域離散化：區間平均
微分運算離散化：有限差分（後向差分／前向差分／中央差分）

有限差分法，無法求得守恆律的弱解。函數跳躍處，泰勒級數某些項不存在，截斷誤差不受控制，收斂性不受控制。

有限體積法，可以求得守恆律的弱解。未知函數切成小段區間，每一段區間的定積分的平均數，作為遞迴公式的變數。定積分的平均數，不受跳躍影響。

另一方面，讓這些區間的定積分的平均數互通有無、彼增我減，導致未知函數的全域的定積分維持定值，導致物理量守恆。

volume / quantity / flux：切割空間，形成體積。體積內部函數值總和，形成數量。體積之間函數值流動，形成通量。

volume：體積。物理量對空間定積分，定積分的範圍。
quantity：數量。物理量對空間定積分，定積分的結果。
flux：通量。數量對時間微分。

volume discretization：體積有兩種設定方式，第一種是鄰邊圍住的範圍，第二種是Voronoi diagram。此處採用第一種。

(1) cell-centered finite volume method
(2) vertex-centered finite volume method

conservation law：守恆律。物理量守恆。函數值總和是定值。函數值四處流動、不生不滅。

conservation law
d           d
—— u(t,x) + —— f(u(t,x)) = 0
dt          dx

conservation law
d/dt u = d/dx f(u)

integral form：積分型。守恆律做定積分，形成常微分方程式。定積分範圍，以xᵢ為中心，左至x_i-1/2，右至x_i+1/2，形成體積。定積分結果，形成數量和通量。

守恆律被改寫成：先後數量差異＝相鄰通量差異。

積分型沒有實際用途。積分型只是用來理解體積數量通量。

conservation law
d           d
—— u(t,x) + —— f(u(t,x)) = 0
dt          dx

definite integral with respect to x from xᵢ₋₁⸝₂ to xᵢ₊₁⸝₂
⌠xᵢ₊₁⸝₂ ⎛ d           d            ⎞
⎮       ⎜ —— u(t,x) + —— f(u(t,x)) ⎟ dx = 0
⌡xᵢ₋₁⸝₂ ⎝ dt          dx           ⎠

derivation
d  ⎛ ⌠xᵢ₊₁⸝₂           ⎞
—— ⎜ ⎮       u(t,x) dx ⎟ + f(u(t,xᵢ₊₁⸝₂)) - f(u(t,xᵢ₋₁⸝₂)) = 0
dt ⎝ ⌡xᵢ₋₁⸝₂           ⎠

integral form (semi-discretization)
d 
—— uᵢ = fᵢ₋₁⸝₂ - fᵢ₊₁⸝₂
dt

                     ⌠xᵢ₊₁⸝₂
define quantity uᵢ = ⎮       u(t,x) dx
                     ⌡xᵢ₋₁⸝₂

       flux fᵢ₊₁⸝₂ = f(u(t,xᵢ₊₁⸝₂))

conservation form：守恆型。守恆律做離散化，形成遞迴公式。空間取樣，以xᵢ為中心，左取x_i-1/2，右取x_i+1/2。

從積分型的觀點來看，uᵢ⁽ⁿ⁾沒有做定積分，形成數量平均數。

（CFL number的倍率a被挪至通量。對f微分以求得a。）

conservation law
d           d
—— u(t,x) + —— f(u(t,x)) = 0
dt          dx

discretization
u(t+Δt,x) - u(t,x)   f(u(t,x+½Δx)) - f(u(t,x-½Δx))
—————————————————— + ————————————————————————————— = 0
        Δt                         Δx

alternative notation style
uᵢ⁽ⁿ⁺¹⁾ - uᵢ⁽ⁿ⁾   f(uᵢ₊₁⸝₂⁽ⁿ⁾) - f(uᵢ₋₁⸝₂⁽ⁿ⁾)
——————————————— + ——————————————————————————— = 0
       Δt                     Δx

conservation form (full-discretization)
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - λ (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)

define volumetric average of quantity uᵢ⁽ⁿ⁾ = u(t,x)
       flux fᵢ₊₁⸝₂⁽ⁿ⁾ = f(uᵢ₊₁⸝₂⁽ⁿ⁾) = f(u(t,x+½Δx))
       CFL number without speed λ = Δt / Δx

conservation：守恆性。數值解每回合數值總和保持相同。

sum uᵢ⁽ⁿ⁾ = sum uᵢ⁽ⁿ⁺¹⁾

守恆型當中，通量的相鄰差的區間和等於頭尾相減，其餘兩兩抵銷。頭尾是零導致守恆性。

assume f₋₁⸝₂⁽ⁿ⁾ = fₜₐᵢₗ₊₁⸝₂⁽ⁿ⁾ = 0

  sum uᵢ⁽ⁿ⁺¹⁾
= sum { uᵢ⁽ⁿ⁾ - λ (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾) }
= sum uᵢ⁽ⁿ⁾ - λ sum { fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾ }
= sum uᵢ⁽ⁿ⁾ - λ (fₜₐᵢₗ₊₁⸝₂⁽ⁿ⁾ - f₋₁⸝₂⁽ⁿ⁾)
= sum uᵢ⁽ⁿ⁾

也就是說，我們可以自訂通量，反正都會抵消。方才沒有提到怎麼計算通量f_i+1/2⁽ⁿ⁾。我們可以自訂計算公式。例如先算左右通量、再算平均數。例如先算左右函數值平均數、再算通量。

nearest 2 points: average of fluxes
⎰ fᵢ₊₁⸝₂⁽ⁿ⁾ = (fᵢ⁽ⁿ⁾ + fᵢ₊₁⁽ⁿ⁾) / 2
⎱ fᵢ⁽ⁿ⁾ = f(uᵢ⁽ⁿ⁾)

nearest 2 points: flux of average of u
⎰ fᵢ₊₁⸝₂⁽ⁿ⁾ = f(uᵢ₊₁⸝₂⁽ⁿ⁾)
⎱ uᵢ₊₁⸝₂⁽ⁿ⁾ = (uᵢ⁽ⁿ⁾ + uᵢ₊₁⁽ⁿ⁾) / 2

nearest 4 points: flux of weighted average of u
⎰ fᵢ₊₁⸝₂⁽ⁿ⁾ = f(uᵢ₊₁⸝₂⁽ⁿ⁾)
⎱ uᵢ₊₁⸝₂⁽ⁿ⁾ = (uᵢ₋₁⁽ⁿ⁾ + 3uᵢ⁽ⁿ⁾ + 3uᵢ₊₁⁽ⁿ⁾ + uᵢ₊₂⁽ⁿ⁾) / 8

計算公式改寫成通量函數。輸入是兩個相鄰變數。輸入可以推廣成任意個相鄰變數。

conservation form
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - λ (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)

conservation form: flux function
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - λ (f̃(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾) - f̃(uᵢ₋₁⁽ⁿ⁾, uᵢ⁽ⁿ⁾))

define flux fᵢ₊₁⸝₂ = f̃(uᵢ, uᵢ₊₁)
       flux function f̃(uᵢ, uᵢ₊₁) = ......

nearest 2 points: average of fluxes
f̃(uᵢ, uᵢ₊₁) = (f(uᵢ) + f(uᵢ₊₁)) / 2

nearest 2 points: flux of average of u
f̃(uᵢ, uᵢ₊₁) = f((uᵢ + uᵢ₊₁) / 2)

nearest 4 points: flux of weighted average of u
  f̃(uᵢ₋₁, uᵢ, uᵢ₊₁, uᵢ₊₂)
= f((uᵢ₋₁ + 3uᵢ + 3uᵢ₊₁ + uᵢ₊₂) / 8)

consistency：自訂通量必須滿足一致性。

consistency
f̃(uᵢ, uᵢ₊₁) = f(uᵢ₊₁⸝₂)   when Δx→0

consistency
f̃(uᵢ₊₁⸝₂, uᵢ₊₁⸝₂) = f(uᵢ₊₁⸝₂)

consistency
f̃(u*, u*) = f(u*)   where u* is any constant

consistent conservative scheme：一致守恆方案。遞迴公式，同時滿足守恆性與一致性，保證得到弱解。詳情請見後面小節。

（詞性不正確。有人添加and。也有人添加,。）

admissible scheme (entropy satisfying scheme)：容許方案（熵滿足方案）。承上，遞迴公式，同時又滿足離散熵不等式，保證得到容許解（熵解）。詳情請見後面小節。

（詞義不到位。原意是成分單純無添加物的衝擊波。）

numerical flux：如何自訂通量。有限體積法的精髓在於自訂通量。除了方才提到的方式，大家也針對特定的微分方程式，發明特定的通量。例如大家最近開始考慮Euler equations of gas dynamics的kinetic energy守恆、Gibbs entropy守恆。若有興趣請自行深究。

stable and non-dissipative
kinetic-energy and entropy preserving scheme (KEEP)

https://www.iccfd.org/iccfd11/assets/pdf/papers/ICCFD11_Paper-0401.pdf
https://arxiv.org/pdf/2205.05217

numerical sheme：如何操作通量。有限體積法的精髓在於操作通量。嘗試各種操作手法，得到各種數值方案。本文介紹其中四類。

(1) finite difference method & flux limiter
(2) Godunov's method
(3) Godunov's method & slope limiter (MUSCL)
(4) essentially non-oscillatory method (ENO)

finite difference method：先前小節的數值方案可以改寫成守恆型。順帶一提，upwind scheme不守恆，無法改寫成守恆型。

Lax–Friedrichs method
可以收斂至容許解。
fᵢ₊₁⸝₂⁽ⁿ⁾ = (fᵢ⁽ⁿ⁾ + fᵢ₊₁⁽ⁿ⁾) / 2
          - (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) / λ / 2

Lax–Wendroff method
無法收斂至容許解。
fᵢ₊₁⸝₂⁽ⁿ⁾ = (fᵢ⁽ⁿ⁾ + fᵢ₊₁⁽ⁿ⁾) / 2
          - (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) / λ / 2
          + f′ᵢ₊₁⸝₂⁽ⁿ⁾ (fᵢ₊₁⁽ⁿ⁾ - fᵢ⁽ⁿ⁾) * λ / 2

flux limiter：有限差分法，引入flux averaging。微調權重以微調通量。特殊權重可以確保保單調。不幸的是，大家沒有仔細檢查這些數值方案是否收斂至容許解。

專著《Numerical Methods for Fluid Dynamics: with Applications to Geophysics》。

Lax–Wendroff
Warming–Beam
Sweby

Kurganov–Tadmor

https://en.wikipedia.org/wiki/Flux-corrected_transport

Harten's incremental form：

upwind scheme：有限差分法利用upwind scheme，確保正組合，進而確保保單調。有限體積法可以如法炮製。差分地點改成通量地點，以便守恆。

upwind scheme of finite volume method
fᵢ₊₁⸝₂⁽ⁿ⁾ = upwind_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)
          = ⎰ f(uᵢ)   , if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
            ⎱ f(uᵢ₊₁)  , if uᵢ⁽ⁿ⁾ > uᵢ₊₁⁽ⁿ⁾

CFL condition
λ max {|f′(u*)| : uᵢ⁽ⁿ⁾ ≥ u* ≥ uᵢ₊₁⁽ⁿ⁾} ≤ 1

https://dblp.org/pid/151/3708.html

Godunov's method：兩兩相鄰格子皆視作Riemann problem。

可以視作upwind scheme的改良版本。

一、差分地點改成通量地點。以便守恆。

二、單點數值改成區間極值。形成Riemann problem的正解。

順帶一提，當守恆律是advection equation，區間極值退化為單點數值，證明在此。

講義《Numerical Methods for Conservation Laws and Related Equation》。

Godunov's method [Osher 1983]
fᵢ₊₁⸝₂⁽ⁿ⁾ = Godunov_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)

          = ⎧ min {f(u*) : uᵢ⁽ⁿ⁾ ≤ u* ≤ uᵢ₊₁⁽ⁿ⁾}
            ⎨                         , if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
            ⎪ max {f(u*) : uᵢ⁽ⁿ⁾ ≥ u* ≥ uᵢ₊₁⁽ⁿ⁾}
            ⎩                         , if uᵢ⁽ⁿ⁾ > uᵢ₊₁⁽ⁿ⁾

Engquist–Osher method
fᵢ₊₁⸝₂⁽ⁿ⁾ = Engquist_Osher_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)

            f(uᵢ⁽ⁿ⁾) + f(uᵢ₊₁⁽ⁿ⁾)    1  ⌠ uᵢ₊₁⁽ⁿ⁾
          = ————————————————————— − ——— ⎮        |f′(θ)| dθ
                     2               2  ⌡ uᵢ⁽ⁿ⁾

          = f(max(uᵢ⁽ⁿ⁾, ω)) + f(min(uᵢ₊₁⁽ⁿ⁾, ω)) - f(ω)

for convex fluxes with minimum at ω,
Engquist_Osher_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾) = f⁺(uᵢ⁽ⁿ⁾) + f⁻(uᵢ₊₁⁽ⁿ⁾)
where f⁺(u) = f(max(u, ω))
      f⁻(u) = f(min(u, ω))

Rusanov's method: derived from Lax–Friedrichs method
fᵢ₊₁⸝₂⁽ⁿ⁾ = Rusanov_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)

            f(uᵢ⁽ⁿ⁾) + f(uᵢ₊₁⁽ⁿ⁾)
          = —————————————————————
                     2
            uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾
          - ——————————————— max {f′(u*) : uᵢ⁽ⁿ⁾ ≤ u* ≤ uᵢ₊₁⁽ⁿ⁾}
                   2

for convex fluxes,
fᵢ₊₁⸝₂⁽ⁿ⁾ = Rusanov_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)

            f(uᵢ⁽ⁿ⁾) + f(uᵢ₊₁⁽ⁿ⁾)
          = —————————————————————
                     2
            uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾ 
          - ——————————————— max(f′(uᵢ⁽ⁿ⁾), f′(uᵢ₊₁⁽ⁿ⁾))
                   2

Roe's method: approximated and simplified Godunov's method
fᵢ₊₁⸝₂⁽ⁿ⁾ = Roe_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)

          = ⎰ min(f(uᵢ⁽ⁿ⁾), f(uᵢ₊₁⁽ⁿ⁾)) , if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
            ⎱ max(f(uᵢ⁽ⁿ⁾), f(uᵢ₊₁⁽ⁿ⁾)) , if uᵢ⁽ⁿ⁾ > uᵢ₊₁⁽ⁿ⁾

https://ocw.mit.edu/courses/16-920j-numerical-methods-for-partial-differential-equations-sma-5212-spring-2003/bafff8bb84c0ee8ac362cea89c5fea38_lec12_notes.pdf
https://arxiv.org/pdf/2004.02258

Godunov's method：推廣版本。仔細調整水面。三個步驟。

講義《Computational Magnetohydrodynamics》。

Godunov's method [van Leer 1979]
(1) cell reconstruction：每個格子，求出水面形狀。
    use uᵢ⁽ⁿ⁾ get curve from uᵢ₋₁⸝₂⁺ to uᵢ₊₁⸝₂⁻
(2) cell interface flux evaluation：每個相鄰格子介面，求出通量。
    fᵢ₊₁⸝₂⁽ⁿ⁾
(3) update：套用遞迴公式，求出新的水面平均值。
    uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - λ (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)

initialize:
 - define the number of cells (N) and the cell width (Δx).
 - set the initial conditions in the cells.

time-stepping loop:
  for (t = beginning_time; t <= ending_time; t += Δt)
    for each cell
      - compute quantity curve in cell face
        using interpolation.
    for each cell interface
      - compute fluxes at cell interfaces
        using Riemann solvers.
    for each cell
      - update the solution using the fluxes
        at the left and right interfaces of the cell.
        Uᵢ := Uᵢ - (Δt / Δx) * (Fᵢ₊₁⸝₂ - Fᵢ₋₁⸝₂)
        where Uᵢ is the quantity in cell i.
        where Fᵢ₋₁⸝₂ and Fᵢ₊₁⸝₂ are the fluxes
        at the left and right interfaces of cell i.

https://en.wikipedia.org/wiki/MUSCL_scheme
https://www.ita.uni-heidelberg.de/~dullemond/lectures/num_fluid_2011/Chapter_4.pdf
https://core.ac.uk/download/pdf/82677472.pdf
https://www.math.sciences.univ-nantes.fr/~berthon/articles/hancock.pdf

cell reconstruction：一個格子，重新設定水面形狀。各種內插方式，達成各種精度。

(1) piecewise constant interpolation
    形成線性遞迴公式。
    由於保單調，根據Godunov's order barrier theorem，一階精度。

(2) piecewise linear interpolation and slope limiter
    形成非線性遞迴公式。二階精度。

for linear scheme,
(monotonicicy preserving <=> positivity) => monotone

for conservative linear scheme,
(monotonicity preserving & neighborhood <=> LED) => monotone

slope limiter：斜率限制器。每個格子，微調水面形狀，進而確保monotonicity preserving。局部極值們，斜率是零。單調區間們，斜率是左右差分取平緩者。

此方式滿足Osher's E scheme，收斂至容許解，並且滿足L_∞-diminishing與TVD。（引入slope limiter之後成為非線性遞迴公式，需要另外證明。證明預計在此。）【尚待確認】

數學式子有兩種寫法。早期寫成左右差分比值，近期寫成左中差分比值，分別是下面兩篇文章。我個人推薦左中差分比值。

專論《High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws》。

專論《Understand Slope Limiter - Graphically》。

minmod
superbee
van Leer

補零中位數
median(x,y,0) = minmod(x,y)
              = ½ (sgn(x) + sgn(y)) min(x,y)
中位數
median(x,y,z) = x + minmod(y-x, z-x)

https://arxiv.org/abs/2102.04435
https://en.wikipedia.org/wiki/Flux_limiter
https://amrvac.org/md_doc_limiter.html

cell interface flux evaluation：兩個相鄰格子，求出介面通量。視作Riemann problem，給定介面左右兩個高度，求得介面通量。計算公式稱作黎曼求解器Riemann solver，數量眾多，各有利弊。視情況收斂至容許解。

專著《Riemann Solvers and Numerical Methods for Fluid Dynamics》。

fᵢ₊₁⸝₂⁽ⁿ⁾ = f(uᵢ₊₁⸝₂⁽ⁿ⁾(0))
where uᵢ₊₁⸝₂⁽ⁿ⁾(x) is the solution of Riemann problem:
  ũᵢ₊₁⸝₂(x) = ⎰ uᵢ⁺   , if x < xᵢ₊₁⸝₂
              ⎱ uᵢ₊₁⁻ , if x > xᵢ₊₁⸝₂

Godunov solver：正解。前面提到的Godunov flux。判別不連續的類型，然後決定答案。
Harten–Lax–van Leer–Einfeldt solver (HLLE)：近似解。保正。
Harten–Lax–van Leer contact solver (HLLC)：近似解。改良版本。

大家也針對特定的微分方程式，發明特定的Riemann solver。例如Euler equations of gas dynamics的Riemann solver：

Godunov solver：正解。
Roe solver：近似解。需要計算微分矩陣的反矩陣。
Osher solver：近似解。
advection upstream splitting method (AUSM)：近似解。

http://www.ae.metu.edu.tr/tuncer/ae546/docs/NumericalMethods/section3_7.html
https://princetonuniversity.github.io/Athena-Cversion/AthenaDocsUGRiemann.html
http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf

2D Godunov's method：推廣到二維空間。此處不介紹。

講義《Computational Astrophysics》。

semi-Lagrangian method：設定水面形狀、移動整塊水體、分配各格水量。Godunov's method可以視作此演算法的特例。

semi-Lagrangian method
(1) reconstruction：一個格子，重新設定水面形狀。分段內插（加權平均）。
(2) advection：根據各格速度，挪動各格水體。不可超過一格，即是CFL condition。
(3) remapping：重新統計每個格子獲得多少水（加權平均）。

essentially non-oscillatory method (ENO method)：不再內插水面，直接內插通量，Neville interpolation暨upwind scheme。不再考慮L₂精度，直接考慮TVD收縮速度，突破天際。

essentially non-oscillatory method (ENO method)
‖f̂⁽ⁿ⁾‖ᴛᴠ ≤ ‖f̂⁽⁰⁾‖ᴛᴠ + O((Δx)𐞥)
數值解總差值減少量受限
[Harten–Engquist–Osher–Chakravarthy 1987]
https://www.researchgate.net/publication/223605523

weighted essentially non-oscillatory method (WENO method)
內插時，使用特殊權重。
https://en.wikipedia.org/wiki/WENO_methods

multi-dimensional optimal order detection (MOOD)：我沒有研究。

https://hal.science/hal-00566023/file/FVCA6_MOOD_Clain-Diot-Loubere.pdf
https://hal.science/hal-00637123/file/moodCaF.pdf
https://www.lamfa.u-picardie.fr/desveaux/eMOOD.pdf

temporal discretization：習慣採用SSPRK。diagonally implicit strong-stability-preserving Runge-Kutta method現正流行中。

spatial discretization：習慣採用一致守恆方案、容許方案。進一步細分為如何自訂通量、如何操作通量。

consistent conservative sheme

本節討論得到弱解的條件。

根據下面兩個定理，一致守恆方案可以得到弱解。【尚待確認】

一、Lax–Wendroff consistency theorem：一致守恆方案，當取樣間距Δt與Δx趨近零，則數值解趨近弱解。

證明手法是先前章節提到的包裝。守恆律（微分方程式）經過包裝得到新式，一致守恆方案（遞迴公式）經過包裝得到新式，兩道新式趨近相等，導致兩道新式的解也趨近相等。

注意到，此定理有個前提：必須有解。換句話說：符號解／數值解必須滿足存在性。換句話說：符號解／數值解必須收斂至定值。

證明過程請見這本專著：

conservation law with initial condition

   ⎧ du(t,x)   df(u(t,x))
   ⎪ ——————— + —————————— = 0   where t ≥ 0
   ⎨   dt          dx
   ⎪ 
   ⎩ u(0,x) = u₀(x)

after wrapping (u(t,x) is weak solution)

   ⌠+∞ ⌠+∞ ⎛        dϕ(t,x)             dϕ(t,x) ⎞
   ⎮   ⎮   ⎜ u(t,x) ——————— + f(u(t,x)) ——————— ⎟ dt dx
   ⌡0  ⌡-∞ ⎝          dt                  dx    ⎠
 
   ⌠+∞
 + ⎮   u₀(x) ϕ(0,x) dx
   ⌡-∞

 = 0

   (denoted by L(u(t,x)) = 0)

consistent conservative scheme

                     Δt
   uᵢ⁽ⁿ⁺¹⁾ - uᵢ⁽ⁿ⁾ + —— (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾) = 0
                     Δx

   where fᵢ₊₁⸝₂ = f̃(uᵢ,uᵢ₊₁) and f̃(u,u) = f(u)

after wrapping

  +∞ ⌠+∞        ϕ(t+Δt,x) - ϕ(t,x)
   ∑ ⎮   u(t,x) —————————————————— dx Δt
   0 ⌡-∞                 Δt
 
   ⌠+∞
 + ⎮   u₀(x) ϕ(0,x) dx
   ⌡-∞
 
  +∞ ⌠+∞                           ϕ(t,x+½Δx) - ϕ(t,x-½Δx)
 + ∑ ⎮   f̃(u(t,x-½Δx), u(t,x+½Δx)) ——————————————————————— dx Δt
   0 ⌡-∞                                      Δx          

 = 0

   (assume ϕ(t,x)→0 when x→±∞)
   (denoted by L̂(u(t,x)) = 0)

Lax–Wendroff consistency theorem
L(u(t,x)) - L̂(u(t,x)) → 0 when Δt→0, Δx→0

https://hal.science/hal-03257774/document
https://indico.math.cnrs.fr/event/8103/attachments/3538/4802/tlse.pdf

二、[DiPerna 1983]：一致守恆方案，若數值解滿足bounded L_∞與bounded variation穩定性，則數值解滿足L₁收斂性。

一維空間的情況下，bounded variation ⇒ bounded L_∞。光是bounded variation足以滿足L₁收斂性。

實務上大家難以判斷bounded。方便起見，大家換成更強更嚴格的diminishing：一致守恆方案，若數值解滿足bounded L_∞與total variation diminishing，則數值解滿足L₁收斂性。

一維空間的情況下，光是TVD足以滿足L₁收斂性。注意到TVD ⇏ L_∞-diminishing。例如全域收縮、局部放大。

做個總結。一致守恆方案收斂至弱解，當數值解滿足TVD。

Finite volume methods: foundation and analysis
https://ivv5hpp.uni-muenster.de/u/mohlb_01/postscript/finvol_script.pdf

Convergence of approximate solutions of conservation laws
https://www.igpm.rwth-aachen.de/Download/reports/pdf/IGPM217.pdf

[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/
L1-bounded => L2-bounded
L2-convergence => L1-convergence

fundamental theorem of numerical analysis：數值分析基本定理。一致時，穩定即收斂。

根據方才兩個定理，我們可以建立數值分析基本定理的新版本，針對守恆律的弱解。

theorem:     weak solution of conservation law
consistency: consistent conservative scheme
stability:   bounded L_∞ and bounded variation stability
convergence: L₁-convergence

守恆律暨弱解之數值分析基本定理
一致性：一致守恆方案
穩定性：bounded L_∞暨bounded variation
收斂性：L₁收斂性

consistent conservative scheme：守恆型的flux必須等於原本數值。

conservative:

consistent:
fᵢ₊₁⸝₂⁽ⁿ⁾ = flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)
where flux(uᵢ⁽ⁿ⁾, uᵢ⁽ⁿ⁾) = f(uᵢ⁽ⁿ⁾)

bounded L_∞ and bounded variation stability：數值解每回合L_∞-norm和total variation大小受限。

⎰ ‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖_∞  ≤ B_∞   when n≥0
⎱ ‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ᴛᴠ ≤ Bᴛᴠ

where ‖f‖_∞  = maxᵢ |fᵢ|
where ‖f‖ᴛᴠ = sumᵢ |fᵢ₊₁ - fᵢ|

L₁-convergence：數值解與符號解每回合L₁距離趨近零。

‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖₁ → 0 when n≥0

where ‖f‖₁ = sumᵢ |fᵢ|

L₁-norm是絕對值總和。若數值非負，L₁-norm是總和。

當符號解／數值解非負，L₁-norm趨近相等可以視作物理量守恆。差強人意。

u                        u       |
|     | |                |     | | |
|   | | | | |            |   | | | | |
| | | | | | | | |        | | | | | | |     
|-----------------> x    |-----------------> x
 symbolic solution        numeric solution

注意到，此定理僅適用於一維純量守恆律。此定理無法直接用於多維空間、向量、方程組，無法直接引入散度、旋度、梯散，無法直接指定取樣間距、網格造型、邊界條件。必須另行證明或反駁。

守恆律的數值分析基本定理是數學界的大難題。每種守恆律都必須特地檢查存在性、唯一性、一致性、穩定性、收斂性。

admissible scheme（entropy satisfying scheme）

本節討論得到容許解的條件。

entropy inequality：符號解滿足熵不等式則是容許解。

使用了先前章節提到的包裝。容許解是弱解。

小於等於是為了收縮。

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
弱解的其中一個特例：容許解。

然而，有些數值方案，無法顧及平滑函數，僅顧及p次可微函數，導致一致性部分成立。守恆律（微分方程式）、一致守恆方案（遞迴公式）僅p-norm趨近相等。【尚待確認】

consistency with pth-order differentiable test function
‖L(u(t,x)) - L̂(u(t,x))‖ₚ → 0 when Δt→0, Δx→0

先前提到的Oleinik admissibility condition其實是熵不等式其中一種特例：採用Kruzhkov entropy。一致性僅L₁-norm趨近相等。

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₁-convergence。
（用來證明容許解是唯一解。）

Lₚ entropy
η(u) = |u|ᵖ / p
其中一種凸函數。導致Lₚ-convergence。

一致性僅L₁-norm趨近相等，穩定性和收斂性亦然。

discrete entropy inequality：數值解滿足離散熵不等式則是容許解。

conservation
sumᵢ ûᵢ⁽ⁿ⁾ = sumᵢ ûᵢ⁽ⁿ⁺¹⁾

quantity average
      1 ⌠ xᵢ₊₁⸝₂
ûᵢ = —— ⎮        u(tₙ,x) dx   when Δx→0
     Δx ⌡ xᵢ₋₁⸝₂

consistency
f̂(ûᵢ,ûᵢ,...,ûᵢ) = f(ûᵢ)

consistent conservative scheme
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - λ (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)
where λ = Δt / Δx
      f̃(u,u) = f(u)
      fᵢ₊₁⸝₂ = f̃(uᵢ,uᵢ₊₁)

discrete entropy inequality [Tadmor 1984]
η(uᵢ⁽ⁿ⁺¹⁾) ≤ η(uᵢ⁽ⁿ⁾) - λ (qᵢ₊₁⸝₂⁽ⁿ⁾ - qᵢ₋₁⸝₂⁽ⁿ⁾)
where entropy pair (η,q)
      q̃(u,u) = q(u)
      qᵢ₊₁⸝₂ = q̃(uᵢ,uᵢ₊₁)

(η′ʟ - η′ʀ) flux(uʟ, uʀ) ≥ qʟ - qʀ

admissible scheme (entropy satisfying scheme)：一致守恆方案，滿足離散熵不等式、滿足收縮穩定性（連帶滿足CFL condition），則收斂至容許解。目前有兩類遞迴公式滿足離散熵不等式。

entropy satisfying scheme
目前有兩類遞迴公式可以得到容許解。
(1) Osher's E scheme (upwind scheme of finite volume method)
    有限體積法版本的上風方案。滿足離散熵不等式所有凸函數，達成所有收斂性。
    針對線性遞迴公式，至多一階精度。
(2) monotone scheme
    單調方案。滿足離散熵不等式其中一種凸函數Kruzhkov entropy，達成L₁收斂性。
    原始論文沒有證明其他凸函數，不清楚是否達成所有收斂性。
    針對線性遞迴公式，至多一階精度。

[Tadmor 1984]
upwind Lax–Friedrichs method符合此式，得到容許解。

[Osher 1985]
Godunov's method符合此式，得到容許解。
MUSCL符合此式，得到容許解。

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
https://www.uio.no/studier/emner/matnat/math/MAT4380/v15/pensumliste/chapter3.pdf

Osher's E scheme：有限體積法版本的上風方案。通量不逾越Godunov flux，以便正組合、遞減、衝擊波。另外，針對線性遞迴公式，正組合等同保單調。

the Godunov's method is Osher's E scheme under CFL ≤ 1.

Osher's E scheme
sgn(uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) (fᵢ₊₁⸝₂⁽ⁿ⁾ - f(u*)) ≤ 0
for all u* ∈ [uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾]

Osher's E scheme
⎧ fᵢ₊₁⸝₂⁽ⁿ⁾ ≤ min {f(u*) : uᵢ⁽ⁿ⁾ ≤ u* ≤ uᵢ₊₁⁽ⁿ⁾}
⎨                                     , if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
⎪ fᵢ₊₁⸝₂⁽ⁿ⁾ ≥ max {f(u*) : uᵢ⁽ⁿ⁾ ≥ u* ≥ uᵢ₊₁⁽ⁿ⁾}
⎩                                     , if uᵢ⁽ⁿ⁾ > uᵢ₊₁⁽ⁿ⁾

Osher's E scheme
sgn(uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) (fᵢ₊₁⸝₂⁽ⁿ⁾ - Godunov_flux(uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾)) ≤ 0
for all u* ∈ [uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾]

Osher's E scheme
⎰ fᵢ₊₁⸝₂⁽ⁿ⁾ ≤ Godunov_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾) , if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
⎱ fᵢ₊₁⸝₂⁽ⁿ⁾ ≥ Godunov_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾) , if uᵢ⁽ⁿ⁾ > uᵢ₊₁⁽ⁿ⁾

Osher's E scheme => monotone scheme【尚待確認】

when f is convex,
monotone scheme => Osher's E scheme【尚待確認】

Osher's E scheme => L_∞-diminishing and TVD

E系列一覽。標準衝擊波。括號附註是近期名稱，但是名稱也說不上精準。

E condition (Oleinik admissibility condition)
守恆型的約束條件，得到標準衝擊波。

entropy solution (admissibile solution)
守恆型的符號解／數值解，恰是標準衝擊波。

entropy inequality
收縮的標準衝擊波，寫成弱解形式。

E scheme (Osher's E scheme)
守恆型的數值方案，得到收縮的標準衝擊波。

monotone scheme：單調方案。注意到不一定收縮。可以用來處理非收縮的數值解。

   consistent & conservative & monotone
=> Kruzhkov entropy inequality
=> L₁-contraction

[Harten–Hyman–Lax 1976]

   conservative & monotone
=> L₁-contraction

[Crandall–Majda 1980]

跳過中間Kruzhkov entropy inequality，直接得到L₁-contraction。
但是這樣就不能證明得到容許解了。
一、f ≤ max(f,g) => F(f) ≤ F(max(f,g))
二、數值解相減，結果有正有負，拆開正負，分開累計。

   consistent & conservative & TVD
=> L₁-contraction
=> converge to weak solution

[LeVeque 1992] theorem 15.2
基本上胡說八道

https://www.researchgate.net/publication/285599438
https://www.researchgate.net/publication/230676056

一般來說，通量函數是兩個變數（遞迴公式右式是三個變數）。此時通量函數對第一個變數偏微分為正數或零、對第二個變數偏微分為負數或零。

for consistent conservative scheme,
                            ⎧ d                   
                            ⎪ ——— f̃(uᵢ, uᵢ₊₁) ≥ 0 
                            ⎪ duᵢ                 
3-point monotone scheme <=> ⎨ 
                            ⎪ d
                            ⎪ ————— f̃(uᵢ, uᵢ₊₁) ≤ 0
                            ⎩ duᵢ₊₁
[Harten–Hyman–Lax 1976]

此時單調則上風。若數值解收縮（連帶滿足CFL condition），單調即上風。

for consistent conservative scheme,
3-point Osher's E scheme <=> 3-point monotone scheme
under CFL condition

(=>) [Abdi–Hansen–Schroll 2020]
An adaptive viscosity E-scheme for balance laws
http://www.math.ualberta.ca/ijnam/Volume-17-2020/No-3-20/2020-03-08.pdf

(<=)
https://math.tifrbng.res.in/~praveen/slides/gian2017_estable.pdf
since g(u,v) is non-decreasing in u and non-increasing in v
⎰ g(u,v) ≤ g(u,w) ≤ g(w,w) = f(w)   if u ≤ w ≤ v
⎱ g(u,v) ≤ g(w,v) ≤ g(w,w) = f(w)   if u ≤ w ≤ v
sign(v - u)(g(u,v) - f(w)) ≤ 0   for all w ∈ [u,v]

(<=)
https://people.math.ethz.ch/~hiptmair/tmp/NUMHYP_07.pdf
since F(v,w) is non-decreasing in v and non-increasing in w
⎰ F(v,w) - F(u,u) ≤ 0   if v < u < w
⎱ F(v,w) - F(u,u) ≥ 0   if w < u < v

finite element method

有限元素法，其意義隨著時代改變。

一、連續物體，大量取樣，形成離散地點。一個地點，根據物理定律列出一個等式，通常是線性方程式。所有地點，形成線性方程組。線性方程組求解，以便得到特定物理量。線性方程組求解，可以重新改寫成大型稀疏矩陣求解，這就是最初的有限元素法。

二、數學家重新彙整觀念。未知函數的離散近似，重新稱作有限差分法。未知函數的連續近似，重新稱作有限體積法。未知函數的取樣擇鄰，沒有適當的詞彙（當時，網格是計算機繪圖的新興詞彙），只好稱作有限元素法。使用自定義網格，這就是後來的有限元素法。

三、未知函數的近似，也可以採用全部地點的shape function的加權總和。距離太遠的地點的shape function數值幾乎是零，只需累計相鄰地點的shape function。針對這件事，沒有適當的詞彙（當前，RBF內插是機器學習的新興詞彙），只好視作有限元素法。使用shape function，這就是目前的有限元素法。

學術名詞是經年累月風吹日曬而成，其意義通常不是原來的樣子。舉例來說，dynamic programming最初是最佳化演算法，經常跟linear programming相提並論。現在則是演算法設計方法（遞迴公式求值方法），經常跟greedy method相提並論。

在工程領域，有限差分法FDM、有限體積法FVM、有限元素法FEM，大家經常將三者相提並論。

有限差分法：單一樣本，函數值。
有限體積法：相鄰樣本，定積分的平均數。
有限元素法：所有樣本，shape function的加權總和。

finite element method

有限元素法。數值解平滑內插，考慮所有取樣地點。

domain discretization: grid/mesh/particles
codomain discretization: shape function
operator discretization: chain rule

定義域離散化：方格／網格／粒子
對應域離散化：形狀函數
微分運算離散化：形狀函數套用微分連鎖律，求得數學式子。

有限差分法的視角：微分運算離散化考慮所有取樣地點，導數平滑內插。
有限體積法的視角：通量函數的輸入囊括所有取樣地點，通量平滑內插。

shape function：高度為一、總和為一的基底函數。

講義《Computational Magnetohydrodynamics》。

Galerkin finite element method：微分方程式與基底函數的點積都是零。

consistency：符號解與數值解的一致性。

專著《Partial Differential Equations》。

基底內插可以視作線性代數的「垂直投影到子空間」，子空間的基底是各個位置的shape function。保證存在唯一解。

但是我不知道截斷誤差是多少。我不知道如何提升精度。

Sobolev space
https://en.wikipedia.org/wiki/Sobolev_space
Lax—Milgram theorem / Babuška–Lax–Milgram theorem
https://en.wikipedia.org/wiki/Babuška–Lax–Milgram_theorem

numerical scheme：數值解的數學性質、其對應的數值方案，目前沒有定論。這是最接近此主題的台灣學校課程。

discontinuous Galerkin method (DG method)
Nitsche's method
https://relate.cs.illinois.edu/course/cs555-s24/

https://users.oden.utexas.edu/~leszek/classes.html
https://users.oden.utexas.edu/~leszek/classes/EM394H/book.pdf

discontinuous Galerkin method：先前小節的數值方案可以改寫成基底內插函數。

High Order Entropy Stable Discontinuous Galerkin Discretizations
https://www.nas.nasa.gov/pubs/ams/2024/08-13-24.html

High order positivity-preserving entropy stable discontinuous Galerkin discretizations
https://sites.google.com/view/jessechan/talks

monolithic convex limiting / invariant domain preserving
凸的下回合還是凸的。
一個例子是maximum principle，bound在極大極小值之間。

Ladyzhenskaya–Babuška–Brezzi condition (inf-sup condition)
https://en.wikipedia.org/wiki/Ladyzhenskaya–Babuška–Brezzi condition
https://math.stackexchange.com/questions/2886310/
https://web.stanford.edu/class/cme358/notes/cme358_lecture_notes_2.pdf

一些古怪的名詞。

reconstruct–evolve–average method (REA method)
即是semi-Lagrangian method
LeVeque

monotonic upstream-centered scheme for conservation laws (MUSCL)
Godunov's method追加slope limiter
van Leer

discontinuous Galerkin method (DG method)
Godunov's method追加shape function
Cockburn、舒其望

arbitrary high order using derivatives (ADER-DG method)
Godunov's method追加high order shape function
Toro

high-resolution scheme
即是higher-order numerical scheme
大於一階精度的數值方案們。例如MUSCL、DG method。
Harten

spectral method

頻譜法。實施傅立葉轉換，一切挪至頻域處理。

domain discretization: grid
codomain discretization: sinusoid
operator discretization: minimal oscillation trigonometric interpolation

定義域離散化：方格
對應域離散化：弦波
微分運算離散化：最小振盪的三角函數內插

頻譜法：微分方程式實施傅立葉轉換，移項整理成遞迴公式。時間離散化沿用有限差分法，才能製作動畫。

constant-coefficient advection equation:
d             d
—— u(t,x) + c —— u(t,x) = 0
dt            dx

semi-discretization (explicit Euler method):
u(t+Δt,x) - u(t,x)     d
—————————————————— + c —— u(t,x) = 0
         Δt            dx

spectral method:
   ⎛ u(t+Δt,x) - u(t,x)     d         ⎞
DFT⎜ —————————————————— + c —— u(t,x) ⎟ = 0
   ⎝          Δt            dx        ⎠

改寫成陣列：
(DFT(u[n+1]) - DFT(u[n])) / Δt + c (𝑖ωi/Δx) DFT(u[n]) = 0
DFT(u[n+1]) = (1 + C 𝑖ωi) DFT(u[n])
u[n+1] = IDFT((1 + C 𝑖ωi) DFT(u[n]))

where C = c Δt/Δx is CFL number without absolute value
      𝑖 is imaginary unit
      ω = 2π/N is fundamental angular frequency of DFT
      i = (0, 1, ⋯, N-1) is a sequence from 0 to N-1
      1 = (1, 1, ⋯, 1) is a sequence of 1s.
      N is sampling number

也可以全部離散化。優點是adaptive mesh refinement容易收斂。【尚待確認】

constant-coefficient advection equation:
d             d
—— u(t,x) + c —— u(t,x) = 0
dt            dx

full-discretization:
u(t+Δt,x) - u(t,x)     u(t,x+Δx) - u(t,x-Δx)
—————————————————— + c —————————————————————  = 0
         Δt                     2 Δx

spectral method:
   ⎛ u(t+Δt,x) - u(t,x)     u(t,x+Δx) - u(t,x-Δx) ⎞
DFT⎜ —————————————————— + c ————————————————————— ⎟ = 0
   ⎝          Δt                     2 Δx         ⎠

改寫成陣列：
(DFT(u[n+1]) - DFT(u[n])) / Δt
+ c (exp(+𝑖ωi) DFT(u[n]) - exp(-𝑖ωi) DFT(u[n])) / (2Δx) = 0
DFT(u[n+1]) = (1 - ½ C (exp(+𝑖ωi) + exp(-𝑖ωi)) DFT(u[n])
u[n+1] = IDFT((1 - ½ C (exp(+𝑖ωi) + exp(-𝑖ωi)) DFT(u[n]))

傅立葉轉換的數學公式。【尚待確認：索引值標記法】

1. addition <-> addition
   DFT(f[i] + g[i]) = DFT(f[i]) + DFT(g[i])
2. multiplication <-> circular convolution
   DFT(f[i] × g[i]) = DFT(f[i]) ⊛ DFT(g[i])
3. circular shift <-> amplitude modulation
   DFT(f[i⊕1]) = exp(𝑖ωi) DFT(f[i])
4. circular backward difference <-> amplitude modulation
   DFT((f[i] - f[i⊖1]) / Δx) = ((1 - exp(-𝑖ωi)) / Δx) DFT(f[i])
5. differentiation <-> angular acceleration
   DFT(d/dx f(x)) = (𝑖ωi/Δx) DFT(f[i])

https://brianmcfee.net/dstbook-site/content/ch06-dft-properties/Shifting.html
https://dsp.stackexchange.com/questions/70113/
https://math.mit.edu/~stevenj/fft-deriv.pdf

傅立葉轉換的某些數學公式是循環運算，必須小心應對。

(1) 週期邊界條件：無需額外處理，恰好形成循環運算。
(2) 對稱邊界條件：兩種方法。
  (a) method of reflection (DFT)
      追加N-2個鬼點，以便形成循環運算。
      若是實函數，傅立葉轉換函式庫採用實數到複數版本，減少計算時間。
      https://github.com/indutny/fft.js/
  (b) real even-symmetric function (type-II DCT)
      傅立葉轉換替換成餘弦轉換第二型。
      僅適用於實函數。此時等價於第一種方法，但是計算時間更少。
      https://dsp.stackexchange.com/questions/2807/
(3) 開放邊界條件：需要大動土木。我放棄行不行。

original array u[0⋯N-1] under symmetric boundary condition:
⎰ u[-1] = u[1]
⎱ u[N] = u[N-2]

extended array u⃡[0⋯2N-3] under periodic boundary condition:
u⃡[0] = u⃡[2N-2]

method of reflection:
⎰ u⃡[0⋯N-1] = u[0⋯N-1]
⎱ u⃡[N⋯2N-3] = u[N-2⋯1]

method            | number  | naive boundary condition
------------------|---------|-------------------------------
Fourier transform | complex | periodic
sine transform    | real    | odd-symmetric (zero Dirichlet)
cosine transform  | real    | even-symmetric (zero Neumann)

consistency：頻譜法不需要自行調整截斷誤差。

傅立葉轉換：截斷誤差是零。雙射函數，精準無誤。
加法：截斷誤差是零。數學公式，精準無誤。
乘法：同上。
循環移位：同上。
循環前向差分：同上。
微分：截斷誤差源自三角函數內插。
　　　已經保證誤差盡量小，故不需要自行調整誤差數量級。

stability：von Neumann stability analysis。

Schwartz space

numerical scheme：先前小節的數值方案可以套用頻譜法。

(1) spectral difference method
(2) spectral volume method
(3) spectral element method

其中隱方法有極大優勢，時間複雜度較低。以一維空間為例：

有限差分法　隱方法　每回合O(NNT)  　取決於矩陣求解演算法
　　　　　　顯方法　每回合O(N)    　達到理論下限
頻譜法　　　隱方法　每回合O(NlogN)　取決於傅立葉轉換演算法
　　　　　　顯方法　每回合O(NlogN)　取決於傅立葉轉換演算法

Galerkin method：頻譜法其實可以視作Galerkin method。

(1) spectral method
(2) pseudo-spectral method

spectral method：傅立葉轉換是輸入數列與N種弦波分別點積。頻譜法是微分方程式與N種弦波分別點積等於零。因此頻譜法是Galerkin method。

spectral method is Galerkin method:

⎧ ⎛ u(t+Δt,x) - u(t,x)     d         ⎞
⎪ ⎜ —————————————————— + c —— u(t,x) ⎟ ∙ exp(𝑖ω0i) = 0
⎪ ⎝          Δt            dx        ⎠
⎪
⎪ ⎛ u(t+Δt,x) - u(t,x)     d         ⎞
⎨ ⎜ —————————————————— + c —— u(t,x) ⎟ ∙ exp(𝑖ω1i) = 0
⎪ ⎝          Δt            dx        ⎠
⎪              :
⎪              :
⎪ ⎛ u(t+Δt,x) - u(t,x)     d         ⎞
⎪ ⎜ —————————————————— + c —— u(t,x) ⎟ ∙ exp(𝑖ω(N-1)i) = 0
⎩ ⎝          Δt            dx        ⎠

https://www.mech.kth.se/~ardeshir/courses/literature/Notes_Spectral_Methods.pdf

pseudo-spectral method：傅立葉轉換改成小波轉換，N種弦波改成正規正交多項式。因此虛擬頻譜法也是Galerkin method。

method of characteristics

特徵線法是一個非常具有特色的方法，但是被易容。

特徵即是速度向量場。特徵曲線即是粒子移動軌跡。

characteristics = velocity vector field
characteristic curve = pathline

特徵線法：一、手工計算速度向量場的數學公式。二、整個空間均勻散布大量粒子。三、每回合每個粒子依照速度向量場移動一步。

numerical scheme：本文介紹其中三類。

(1) Lagrangian particle method：僅使用粒子。
(2) semi-Lagrangian method：粒子搭配靜止網格。
(3) arbitrary Lagrangian-Eulerian method：粒子搭配動態網格。

物理學當中，特徵線法之粒子移動，稱作「拉格朗日描述Lagrange description」，有限體積法之格子流動，稱作「歐拉描述Euler description」。因此這三類方法取名為拉格朗日。

來張圖片

大家也針對特定的微分方程式，發明特定的數值方案。

(1) smoothed particle hydrodynamics (SPH):
    針對流體方程組。粒子套用shape function，變成有限元素法。
(2) vortex method:
    針對流體方程組。
(3) material point method (MPM):
    針對柔體方程組。速度用粒子來算，應力用格子來算。
(4) particle mesh method:
    針對N-body simulation。粒化格，解Poisson equation，採用頻譜法。
(5) particle–particle–particle–mesh (P³M):
    針對N-body simulation。

consistency：Courant–Friedrichs–Lewy condition。

symplectic method (semi-implicit method)：全域誤差受限且振盪。常見做法是operator splitting method，重新安排式子順序。

專著《Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations》。

     d/dt f(t) = f(t) + g(t)
---> 1. (f[n+1] - f[n]) / Δt = g[n+1]
     2. (f[n+1] - f[n]) / Δt = f[n+1]

Verlet method
leapfrog method: kick-drift-kick / drift-kick-drift
Hermite method

https://en.wikipedia.org/wiki/Leapfrog_integration
https://arxiv.org/abs/0708.0738

level set method

等高線法是一個非常具有特色的方法，但是被雪藏。

一、數學家至今仍未發掘等高線的性質和運算。

二、計算學家至今仍未發明高效的資料結構和演算法。

導致大家在實務上使用其他方法。

等高線法：未知函數的等高線，每回合沿著表面垂直方向移動（梯度方向）。等高線每一點的移動速度取決於守恆律的倍率大小。所有等高線可以改寫成一個有向距離場，其資料結構是三維陣列。

數值模擬過程當中，相同高度的等高線可能融合或分離（宛如水珠融合和分離）。這件事的處理方法，以及收斂性證明，我沒有研究。

Hamilton–Jacobi equation  d/dt f + H(x, df/dx, t) = 0
                          where f has the form f(x(t), t)
eikonal equation          |∇f| = n
level set equation        d/dt φ + v|∇φ| = 0

https://math.mit.edu/classes/18.086/2007/levelsetpres.pdf
https://math.mit.edu/classes/18.086/2007/levelsetnotes.pdf
https://math.mit.edu/classes/18.376/Notes/LectureTopics18376.pdf

https://en.wikipedia.org/wiki/Viscosity_solution
https://benjaminmoll.com/wp-content/uploads/2019/07/viscosity_slides.pdf   

In the case of hyperbolic systems, the notion of weak
solution based on distributions does not guarantee
uniqueness, and it is necessary to supplement it with
entropy conditions or some other selection criterion.

In fully nonlinear PDE such as the Hamilton–Jacobi
equation, there is a very different definition of weak
solution called viscosity solution.

level set ghost point
https://www.researchgate.net/publication/334643942

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

Fronts propagating with curvature-dependent speed
https://math.berkeley.edu/~sethian/Papers/sethian.osher.88.pdf

capturing multivalued solution

Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations
https://math.berkeley.edu/~sethian/Papers/sethian.osher.88.pdf

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
https://www.researchgate.net/publication/2855502

differential equation - in mathematics🚧

微分方程式的分類：多變數函數

微分方程式分為兩類。一、常微分方程式：對同一種變數進行微分。二、偏微分方程式：對多種變數進行微分。

ordinary differential equation, ODE:

             ∂             ∂  ∂ 
f(x,y) + 2 = —— f(x,y) + 3 —— —— f(x,y) + 2 g(x,y)
             ∂x            ∂x ∂x

        ∂f     ∂²f
f + 2 = —— + 3 ——— + 2g     省略括號的部分、合併多次微分的部分
        ∂x     ∂x²          (Leibniz)

f + 2 = Dₓf + 3Dₓₓf + 2g    微分簡寫成 Dₓf Dₓₓf (Arbogast)

f + 2 = fₓ + 3fₓₓ + 2g      微分簡寫成 fₓ fₓₓ   (???)

f + 2 = f′ + 3f″ + 2g       微分簡寫成 f′ f″    (Lagrange)

f + 2 = ḟ + 3f̈ + 2g         微分簡寫成 ḟ f̈      (Newton)

partial differential equation, PDE:

             ∂             ∂  ∂ 
f(x,y) + 2 = —— f(x,y) + 3 —— —— f(x,y) + 2 g(x,y)
             ∂x            ∂x ∂y

        ∂f     ∂²f
f + 2 = —— + 3 ——— + 2g     省略括號的部分、合併多次微分的部分
        ∂x     ∂xy

f + 2 = ∂ₓf + 3∂ₓ₝f + 2g    微分簡寫成 ∂ₓf ∂ₓ₝f

f + 2 = fₓ + 3fₓ₝ + 2g      微分簡寫成 fₓ fₓ₝

微分方程式的分類：多變量函數

微分方程式分為兩類。一、純量微分方程式：純量函數的微分方程式。二、向量微分方程式：向量函數的微分方程式，得以引入梯度、散度、旋度。大家經常省略字尾-valued。

scalar-valued differential equation:

      ∂f     ∂f     ∂²f      ∂²f     ∂²f
f + 2 —— + 3 —— + 5 ———— + g ——— + h ——— = 0
      ∂x     ∂y     ∂x∂y     ∂x²     ∂y²

vector-valued differential equation:

  ∂Fₓ   ∂F₝          ∂F₝   ∂Fₓ 
( ——— + ——— ) Fₓ + ( ——— - ——— ) F₝ = 0 
  ∂x    ∂y           ∂x    ∂y

微分方程式的分類：加法分離

微分方程式分為兩類。一、一次微分方程式：函數微分視作變數，形成一次函數。二、非一次微分方程式。

一次微分方程式，可以推導符號解（分離變數法、格林函數）、演算數值解（時域一次方程組、頻域傅立葉轉換）。非一次微分方程式，則是數學界的大難題，至今只有少數特例找到了符號解。

linear differential equation:

      ∂f     ∂f     ∂²f      ∂²f     ∂²f
f + 2 —— + 3 —— + 5 ———— + g ——— + h ——— + c = 0
      ∂x     ∂y     ∂x∂y     ∂x²     ∂y²

nonlinear differential equation:

  ∂f   ∂f ∂f   ⎛∂f⎞²     ⎛∂f⎞
f —— + —— —— + ⎜——⎟ + sin⎜——⎟ = 0
  ∂x   ∂x ∂y   ⎝∂x⎠      ⎝∂x⎠

微分方程式還有兩類：一、齊次微分方程式：函數微分視作變數，沒有常數項。二、擬一次微分方程式：函數微分視作變數，係數不是常數函數而是泛函數，泛函數的輸入變數不含最高次微分。

齊次微分方程式、擬一次微分方程式，容易推導符號解。

homogeneous differential equation:

      ∂f     ∂f     ∂²f      ∂²f     ∂²f
f + 2 —— + 3 —— + 5 ———— + g ——— + h ——— = 0
      ∂x     ∂y     ∂x∂y     ∂x²     ∂y²

quasilinear differential equation:

     ∂f        ∂f
A(f) —— + B(f) —— + C(f) = 0
     ∂x        ∂y

quasilinear differential equation:

 ⎛   ∂f   ∂f⎞ ∂f    ⎛   ∂f   ∂f⎞ ∂f    ⎛   ∂f   ∂f⎞ ∂²f 
A⎜f  ——   ——⎟ —— + B⎜f  ——   ——⎟ —— + C⎜f  ——   ——⎟ ————
 ⎝ , ∂x , ∂y⎠ ∂x    ⎝ , ∂x , ∂y⎠ ∂y    ⎝ , ∂x , ∂y⎠ ∂x∂y

   ⎛   ∂f   ∂f⎞ ∂²f    ⎛   ∂f   ∂f⎞ ∂²f    ⎛   ∂f   ∂f⎞
+ D⎜f  ——   ——⎟ ——— + E⎜f  ——   ——⎟ ——— + F⎜f  ——   ——⎟ = 0
   ⎝ , ∂x , ∂y⎠ ∂x²    ⎝ , ∂x , ∂y⎠ ∂y²    ⎝ , ∂x , ∂y⎠

微分方程式的分類：乘法分離

微分方程式分為兩類。一、可分離微分方程式：微分方程式／解恰是多個函數相乘，每個函數輸入變數不重複。二、不可分離微分方程式。

可分離微分方程式，容易推導符號解（分離變數法）、演算數值解（顯式方法）。

separable ordinary differential equation:

      ⎛ ∂f ⎞
A(f) B⎜ —— ⎟ = c     where B is identity functional
      ⎝ ∂x ⎠

separable partial differential equation:

equation L( f(t,x,y) ) = 0
solution f(t, x, y) = f₁(t) f₂(x) f₃(y)

non-separable differential equation:

∂f   f
—— = — + f x
∂x   x

微分方程式的分類：相依

常微分方程式分為三類。一、純時間常微分方程式：左式對時間微分，右式只有時間變數。二、自主常微分方程式：左式對時間微分，右式沒有時間變數。三、以上皆非。

pure-time ODE:

df
—— = 3x² + 2x + 1
dx

pure-time ODE:

  d²f     df
5 ——— + 4 —— = 3x² + 2x + 1
  dx²     dx

autonomous ODE:

df
—— = 3f² + 2f + 1
dx

autonomous ODE:

  d²f     df
5 ——— + 4 —— = 3f² + 2f + 1
  dx²     dx

微分方程式的分類：正定

二階偏微分方程式分為三類：橢圓、拋物、雙曲。圓錐曲線。

  ∂²f     ∂²f      ∂²f     ∂f     ∂f   
a ——— + b ———— + c ——— + d —— + e —— + g f + h = 0
  ∂x²     ∂x∂y     ∂y²     ∂x     ∂y   

elliptic PDE     b² - 4ac < 0
parabolic PDE    b² - 4ac = 0
hyperbolic PDE   b² - 4ac > 0

elliptic PDE     e.g. Poisson equation   fₓₓ + f₝₝ = 0
parabolic PDE    e.g. heat equation      fₜ = fₓₓ
hyperbolic PDE   e.g. wave equation      fₜₜ = fₓₓ

微分方程式的分類：函數運算

algebraic equation：代數方程式。不做微分積分運算。
differential equation：微分方程式。只做微分運算。
integral equation：積分方程式。只做積分運算。
difference equation：差分方程式。只做差分運算。
differential-algebraic equation：微分代數方程式。兩個一起上。
differential-difference equation：微分差分方程式。兩個一起上。
integro-differential equation：積分微分方程式。兩個一起上。

微分方程式的分類：函數輸入運算

functional differential equation：泛函微分方程式。函數輸入做各種運算。
delay differential equation：延遲微分方程式。只做減法運算（減去常數）。
ordinary differential equation：常微分方程式。函數輸入不做運算。

微分方程式求解

微分方程式有多解。舉例來說：

一、常數函數，微分通通是零，答案很多種。
二、散度旋度運算，好比a + b = 1，答案很多種。

微分方程式求解，必須事先確保唯一解。

unisolvence 確保唯一解
convergence 確保求解過程可以得到唯一解

追加限制條件，以確保唯一解。

initial condition  指定解的某處的函數值、導數值
boundary condition 指定解的邊界的位置、函數值、導數值

函數值：確保解的零次項（常數項）、一次項、二次項、……。
導數值：確保解的一次項、二次項、……。

微分方程式的標準格式。

explicit ODE：常微分方程ODE，最高次偏微分，挪到左式。
Cauchy problem：偏微分方程PDE，對時間最高次偏微分，挪到左式。

微分方程式的求解套路。

(1) steady state system
  (a) algebraic equation (direct stiffness method)
  (b) differential equation at n → ∞ (relaxation method)
(2) dynamical system
  (a) differential equation

(1) IVP quadratic differential equation
    --> complex number system
    https://www.ryanjuckett.com/damped-springs/
    https://tutorial.math.lamar.edu/classes/de/FundamentalSetsofSolutions.aspx
(2) BVP eigenvalue problem f′ = k f
    --> eigenfunction/eigenmode
    https://tutorial.math.lamar.edu/classes/de/bvpevals.aspx
(3) BVP differential equation
    --> Laplace transform (spectral method)
    https://tutorial.math.lamar.edu/Classes/DE/IVPWithLaplace.aspx
    https://tutorial.math.lamar.edu/Classes/DE/Laplace_Table.aspx

differential equation           --> fundamental solution
system of differential equation --> eigendecomposition

動態系統的解可以畫成圖形。

time series graph   函數曲線。以時間為主軸。
gradient field plot 移動速度。(dx/dt, dy/dt)向量場。
phase portrait      移動路線。初始條件(x₀,y₀)沿著向量場跑。
bifucation diagram  穩態位置。每種參數的穩態位置。

微分方程式的調整

介紹幾個手法。

linearization：線性化。硬是分開變數，成為穩態解與差異量，成為無常數項一次微分方程式。
method of lines：線法。硬是分開時間變數與空間變數，成為時間離散化、空間離散化，依序求解。
operator splitting：算子分裂法。硬是分開每個運算，成為多道式子，依序求解。
domain decomposition：領域分解法。硬是分開邊界條件，成為多個問題，分別求解、整合答案。

nonlinear differential equation
df     df
—— + f ——  = k
dt     dx

linearization f = f₀ + δf
d                        d
—— (f₀ + δf) + (f₀ + δf) —— (f₀ + δf) = k
dt                       dx

                      df₀      df₀
when f₀ is a solution ——— + f₀ ——— = k
                      dt       dx
dδf       dδf       df₀      dδf
———  + f₀ ———  + δf ——— + δf ——— = 0
dt        dx        dx       dx 

when δf is small
dδf       dδf      df₀
———  + f₀ ——— + δf ——— = 0
dt        dx       dx

linear differential equation with variable δf
⎛ dδf ⎞      ⎛ dδf ⎞   df₀ ⎛    ⎞ 
⎜ ——— ⎟ + f₀ ⎜ ——— ⎟ + ——— ⎜ δf ⎟ = 0
⎝ dt  ⎠      ⎝ dx  ⎠   dx  ⎝    ⎠

https://www.math.arizona.edu/~kglasner/math456/fixedpointlin.pdf

高階微分方程式，化作一階微分方程組。

d³           d²           d                   
——— f(t) = 6 ——— f(t) + 7 —— f(t) + 8 f(t) + 9
dt³          dt²          dt

錯誤方式
⎧ f  = f                           ⎡ f  ⎤   ⎡ 0 1 0 0 ⎤ ⎡ 1  ⎤
⎨ f′ = f′                          ⎢ f′ ⎥   ⎢ 0 0 1 0 ⎥ ⎢ f  ⎥
⎪ f″ = f″                     -->  ⎢ f″ ⎥ = ⎢ 0 0 0 1 ⎥ ⎢ f′ ⎥
⎩ f‴ = 6f″ + 7f′ + 8f + 9          ⎣ f‴ ⎦   ⎣ 9 8 7 6 ⎦ ⎣ f″ ⎦

正確方式                          
⎧ (f )′ = f′                       ⎧ f₀′ = f₁
⎨ (f′)′ = f″                  -->  ⎨ f₁′ = f₂
⎩ (f″)′ = 6f″ + 7f′ + 8f + 9       ⎩ f₂′ = 6f₂ + 7f₁ + 8f₀ + 9

微分方程式，化作代數方程式。

differential equation - in physics🚧

經典的微分方程式：Euler's Identity

引入了歐拉恆等式，作為微分方程式的解。

exponential equation  ḟ = f
oscillation equation  ḟ = 𝑖f

exponential equation

「指數方程式」。

ḟ = f，微分之後還是一樣。符號解是指數函數。微分不動點。

ḟ = γf ，添上倍率γ。γ>1是成長率，γ<1是衰減率。

oscillation equation

「振盪方程式」。指數方程式，令實部為零。

ḟ = 𝑖f，只考慮虛部。

ḟ = 𝑖ωf ，添上角速度ω。

經典的微分方程式：Laplacian

引入了梯散運算。全是人名，紀念古人。

Laplace equation    ∆f = 0      fₓₓ + f₝₝ + f₞₞ = 0
Poisson equation    ∆f = g      fₓₓ + f₝₝ + f₞₞ = g
Helmholtz equation  ∆f = λf     fₓₓ + f₝₝ + f₞₞ = λf

人名後面到底要不要加apostrophe s？我不知道。

Laplace equation

「拉普拉斯方程式」。處處梯散為零。

一坨東西的勢力均衡。

Poisson equation

「泊松方程式」。處處梯散已知。

常見用法是兩邊梯散相同∆f₁ = ∆f₂，f₂已知或∆f₂已知。

兩坨東西的勢力分布相等。其中一坨已知。

梯散反運算。已知梯散g，求得原函數f。

符號解：格林函數疊加而得。

數值解：擁有特殊演算法。時域一次方程組、頻域傅立葉轉換。

Helmholtz equation

「亥姆霍玆方程式」。梯散運算的特徵函數。

梯散運算的特徵函數f是各種頻率的cos和sin波。

符號解：f(x) = e^{𝑖√λ‖x‖} / 4π‖x‖。

經典的微分方程式：spacetime

引入了時間變數、空間變數，形成了守恆定律。

Helmholtz equation  f = ∆f     f   = fₓₓ + f₝₝ + f₞₞
heat equation       ḟ = ∆f     fₜ  = fₓₓ + f₝₝ + f₞₞
wave equation       f̈ = ∆f     fₜₜ = fₓₓ + f₝₝ + f₞₞

注意到，數學家使用牛頓微分符號、梯度／散度／旋度／梯散符號，微分對象涵蓋所有變數。物理學家異於常人，牛頓微分符號只針對時間變數，梯度／散度／旋度／梯散符號只針對空間變數。

物理學家另行發明梯度運算spacetime gradient、梯散運算d'Alembertian，以便針對所有變數。食飽換枵。

Helmholtz equation

「亥姆霍玆方程式」。

f = ∆f，函數等於位勢差。

f = -k²∆f，添上波數k。

heat equation

「熱傳導方程式」。硬生生多出傳導二字。

ḟ = ∆f，變化速度等於位勢差。

ḟ = v∆f，添上擴散速度v。

符號解：f(x) = e^-x²/4t / √t。稱作heat kernel。

wave equation

「波動方程式」。硬生生多出動字。

f̈ = ∆f，加速度等於位勢差。位勢差產生彈簧力，彈簧力產生加速度。

f̈ = v²∆f，添上傳播速度v，有如彈性係數。

方程式重新整理成：時間微分的特徵函數、空間二次微分的特徵函數。解是兩者聯立，稱為特徵模態（eigenmode）。物理意義：駐波。

一維空間的符號解：特徵模態f(x,t) = a sin(x+vt) + b sin(x-vt)，其中a與b取決於初始條件。物理意義：兩個波，振幅a與b，往反方向傳播，疊加之後形成駐波。

二維空間的符號解：取決於邊界形狀。特殊的邊界形狀，才能順利推導符號解。

一維琴弦振動：直線線段
二維薄膜振動：方形、圓形、L形
三維固體振動：懸臂梁、H型鋼、平板

經典的微分方程式：higher-order derivative

引入了合體技：零次微分項（泛函數R(f)）、一次微分項（散度∇∙f）、二次微分項（梯散∆f）、常數函數項（已知函數g）。

分別稱作反應、傳輸、擴散、源／匯。

transport equation                      ḟ = ∇∙f
diffusion equation                      ḟ = ∆f
transport–diffusion equation            ḟ = ∆f + ∇∙f
reaction–diffusion equation             ḟ = R(f) + ∆f
reaction–transport–diffusion equation   ḟ = R(f) + ∇∙f + ∆f

合體技導致難以推導符號解。形成數學領域：動態系統。

經典的微分方程式：transport

引入了物理量、速度。

continuity equation  q̇ + ∇∙(qv) = 0   ∂ₜq + ∂ₓ(qv) + ∂₝(qv) = 0
advection equation   q̇ + v∙∇q = 0     ∂ₜq + vₓ∂ₓq + v₝∂₝q = 0
incompressibility    ∇∙v = 0          ∂ₓvₓ + ∂₝v₝ = 0

transport–diffusion equation  q̇ + ∇∙(qv) - ∇∙(D∇q) = S

注意到，數學偏微分vₓ跟物理學分量vₓ無法辨別。一種解法是更換偏微分表示法，偏微分寫成∂/∂x q或者∂ₓq。另一種解法是更換分量表示法，速度分量寫成u v w。實務上的解法是一起更換。

continuity equation

「連續方程式」。

q̇ + ∇∙q = 0，時間變化量等於空間轉移量。

q̇ + ∇∙(qv) = 0，添上每一處的轉移速度v。

advection equation

「平流方程式」。

q̇ + ∇∙q = 0，當速度是常數函數，恰是連續方程式。

q̇ + v∙∇q = 0，添上每一處的轉移速度v。

連續方程式可以拆成兩項：平流與壓縮。

微分的乘法公式：前微後不微、前不微後微。

q̇ + ∇∙(qv) = 0

q̇ + v∙(∇q) + q(∇∙v) = 0
    ^^^^^^   ^^^^^^^
 advection   compression

∂        ∂        ∂           
—— q + ( —— qvₓ + —— qv₝ ) = 0
∂t       ∂x       ∂y

∂           ∂         ∂            ∂       ∂          
—— q + ( vₓ —— q + v₝ —— q ) + q ( —— vₓ + —— v₝ ) = 0
∂t          ∂x        ∂y           ∂x      ∂y         
       ^^^^^^^^^^^^^^^^^^^^^   ^^^^^^^^^^^^^^^^^^^
       advection               compression

平流方程式只有時間變化項、平流項。

有人簡寫成大寫D，稱作「物質導數」。

∂           ∂         ∂                   D       
—— q + ( vₓ —— q + v₝ —— q ) = 0   --->   —— q = 0
∂t          ∂x        ∂y                  Dt

平流方程式追加∇∙v = 0，則滿足連續方程式，則物理量守恆。

incompressibility

「不可壓縮」。

∇∙v = 0。出入速度，總和相等。流動順暢，不疾不徐。

transport–diffusion equation

「傳輸－擴散方程式」。

D是擴散係數。S是源／匯，根據正負號。

q̇ + ∇∙(qv) - ∇∙(D∇q) = S
    ^^^^^^   ^^^^^^^   ^^
 transport   diffusion  source/sink

 ^^^^^^^^^   ^^^^^^^^^  ^^^^^^
 1st-order   2nd-order  constant
 derivative  derivative function

經典的微分方程式：hydrodynamics

引入了質量、動量、能量。

continuity equation    ∂/∂t ρv + ∇∙(ρv⊗v) = 0
Euler equation         ∂/∂t ρv + ∇∙(ρv⊗v) + ρg + ∇p = 0
Cauchy equation        ∂/∂t ρv + ∇∙(ρv⊗v) + ρg + ∇∙τ = 0
Navier–Stokes equation ∂/∂t ρv + ∇∙(ρv⊗v) + ρg + ∇p + ∇∙τ = 0

continuity equation for momentum conservation

「連續方程式之動量守恆」。物理量是動量密度q = ρv。

質量m。體積V。密度ρ = m/V。速度v。動量mv。動量密度ρv = mv/V。

二維動量有兩個方向，總共兩道連續方程式，三維動量有三個方向，總共三道連續方程式。有人利用Kronecker product ⊗，簡寫成一道方程式。

連續方程式的由來是泰勒近似。泰勒級數，取零次項與一次項。二次項之後數值極小，物理學家選擇忽略不計。事實上，所有物理公式，其本質都是泰勒級數取到一次項或者取到二次項。

符號解沒有特別取名。

Euler momentum equation

「歐拉動量方程式」。重力密度ρg。壓力p。

動量相聚，產生壓力。

重力加速度g。重力mg。重力密度ρg = mg/V。

力f。表面積A。壓力p = f/A。

符號解是Arnold–Beltrami–Childress flow。

Cauchy momentum equation

「柯西動量方程式」。重力密度ρg。應力τ。

動量相撞、動量相擦，產生應力。

重力加速度g。重力mg。重力密度ρg = mg/V。

力f。表面積A。應力τ = f/A。

符號解我不清楚。

Navier–Stokes equation for momentum conservation

「流體方程組之動量守恆」。大家一起上。

重力密度可以視作「每單位面積的重力之空間導數」，使得右式每一項都有導數。然而重力密度是零階張量，不適合改寫成導數。

流體方程組是慣量、動量、能量三種式子聯立。此處只講動量。

流體方程組的原始版本只針對牛頓流體。牛頓流體是流體的特例。此處版本是後人重新歸納而得。

符號解是千禧年大獎難題，跟P/NP問題齊名。

順帶一提，多維空間的情況下，也許可以把流場調整為無旋場。例如最短路徑演算法Johnson's algorithm的調整路徑手法，把邊權重調整為無負環。然而我們對調整函數知之甚少。等你發表論文。

備忘

牛頓第一方程
1. 質量沿著速度方向移動
2. 速度是質量的附庸品，速度也沿著速度方向移動
3. 動量沿著速度方向移動

PV=NRT方程
空間變小，但是動量總和不變
---> 動量撞牆次數變大，壓力變大 (跟表面積有啥關係?)

Bernoulli equation 白努力方程
1. 1/2 dvv + dgh + p = 常數
2. 1/2 mvv + mgh + pV = 常數
   動能      位能  nrt
3. Euler equation的特例：無旋(有勢)無散(不可壓縮)純粹諧場
   ∂q/∂t + 1/2 (∂vₓ²/∂x + ∂v₝²/∂x) + 1/ρ ∂p/∂x = 0
   勢變化

經典的微分方程式：discontinuity

符號解不連續。

Burgers equation                  衝擊波
traffic-flow equation             衝擊波
Friedlander equation              爆炸波

Burgers equation

wave breaking time。

https://en.wikipedia.org/wiki/Breaking_wave

https://physics.stackexchange.com/questions/76852/

Euler–Poisson equation。組合拳。

traffic-flow equation

http://faculty.washington.edu/rjl/riemann_book/Traffic_flow.html

交通車流問題。車輛不會互相穿越，整個過程都是函數。

經典的微分方程式：soliton

符號解不動點。

Korteweg–de Vries equation        孤立波
shallow water equations           淺水波
Boussinesq's equations            淺水波
de Saint-Venant equations         淺水波

Korteweg–de Vries equation

符號解是solitary wave與cnoidal wave。

記得介紹一下dispersion和group velocity。

Huygens–Fresnel principle。

shallow water equations

專著《Computational Algorithms for Shallow Water Equations》。

Boussinesq's equations

https://commons.wikimedia.org/wiki/File:Celeris_Wave_Model.gif

經典的微分方程式：fractal

符號解自遞迴。

Turing model
Gray–Scott model
FitzHugh–Nagumo model
Hodgkin–Huxley model

Turing model

各式各樣的參數，各式各樣的造型，造就大自然。

Turing pattern

Gray–Scott model

https://pmneila.github.io/jsexp/grayscott/

經典的微分方程式：chaos

符號解無規律。亂繞圓圈的路線。

Lorenz equation                   大氣對流
Lorenz–Emanuel system             大氣對流
Rabinovich–Fabrikant equation     雙角

經典的微分方程式：static equilibrium

萬物靜止，保持均衡。

nodal equilibrium equation

nodal equilibrium equation

方程式清單請見課程講義：

土木工程之結構分析，已經發展一套分析靜態結構的方法。其中跟本文有關的方法是direct stiffness method：運用虎克定律F = kx計算位移量。

force
(1) flexibility method (method of consistency deformation)

displacement
(1) slope-deflection method
(2) moment distribution method
(3) direct stiffness method

所有元件／所有取樣地點，皆滿足虎克定律。列出所有等式，化作線性方程組，得到大型稀疏矩陣。矩陣求解，得到每個元件／每個取樣地點的位移量。整件事情可以想成是：給定常微分方程式，求得數值解。這個方法最初稱作有限元素法，後來又稱作矩陣法。

土木工程之結構分析之有限元素法，應用數學之微分方程之有限元素法，後者源自前者，不過後者已經發展為完全不同的事物。

經典的微分方程式：dynamic equilibrium

萬物消長，走向均衡。

equation of motion in vibration   物理運動
Lotka–Volterra equations          生態消長
Michaelis–Menten equation         化學反應
Hill equation                     化學反應

equation of motion in vibration

F(t) = m d²/dt² x(t) + c d/dt x(t) - k x(t)
where m c k are constants

F = mẍ + cẋ + kx

F: force
x: position (x is displacement when x₀ = 0. now let x₀ = 0.)
m: mass
c: damping coefficient
k: spring constant

彈簧位置符號解：

Lotka–Volterra equations

⎰ d/dt x(t) = a x(t) - b x(t) y(t)
⎱ d/dt y(t) = c x(t) y(t) - d y(t)
where a b c d are constants

⎰ ẋ = ax - bxy
⎱ ẏ = cxy - dy

x: population density of prey (e.g. rabbit)
y: population density of predator (e.g. fox)

此微分方程組可以改造成微分方程式，方便推導符號解。【尚待確認】

dx(t)   a x(t) - b x(t) y(t)
————— = ————————————————————
dy(t)   c x(t) y(t) - d y(t)

三個數學技巧：函數對函數微分、微分連鎖律、反函數定理（函數與反函數在同一點的斜率互為倒數）。三者形成等量微分公理。

                       dx(t)
                       —————
dx(t)   dx(t) dt       dt       a x(t) - b x(t) y(t)
————— = ————— ————— = ——————— = ————————————————————
dy(t)   dt    dy(t)    dy(t)    c x(t) y(t) - d y(t)
                       —————
                       dt

^^^^^^  ^^^^^^^^^^^   ^^^^^^^
1.      2.            3.

1. rate of two differentials
2. chain rule
3. inverse function theorem
1+2+3. axioms of equality

注意到此式無法用來製作動畫。

經典的微分方程式：instability

多種物質的互動。

vorticity equation                渦流
Ginzburg–Landau equation          剪力流
thin-film equation                摩擦流
Buckley–Leverett equation         二相流
Kuramoto–Sivashinsky equation     冷熱交融

Ginzburg–Landau equation

專著《Stability and Transition in Shear Flows》。

thin-film equation

Hele–Shaw equation
https://www.nature.com/articles/ncomms1289

Kuramoto–Sivashinsky equation

經典的微分方程式：relativity

物理元件的互動。

Maxwell equations                 電磁對偶
Schrödinger equation              波粒對偶

Maxwell equations

「電磁場方程組」。科學家們先後發現各式各樣的電磁關係式，由Maxwell集結成一套完整理論，由Heaviside精煉成散度與旋度。

沒有電磁感應的時候，電場靜止。靜電場散度為零。靜電場旋度為零，故靜電場是梯度場／守恆場，故靜電場可以改寫成電勢。

發生電磁感應的時候，電場運動。動電場散度依然為零。動電場旋度不再為零，故動電場不再是梯度場／守恆場。電場空間旋度恰等於磁場時間導數。電場空間變化導致磁場時間變化。電磁對調亦然。

電磁場方程組用來描述動電場。

⎧ ∇∙E = 0
⎨ ∇∙B = 0
⎪ ∇×E = - d/dt B
⎩ ∇×B = - d/dt E

electromagnetic wave equation

「電磁波方程式」。電場振盪產生波。磁場隨之振盪。

電／磁／波，其方向是右手四指／手掌／拇指。形成橫波。

假設只往X軸方向傳輸，導致旋度運算只剩一項。

⎰ d/dx E = - d/dt B
⎱ d/dx B = - d/dt E

一式對空間微分，另一式對時間微分，兩式合併，形成波動方程式。電場版本和磁場版本等價，實際上只有一道方程式。

∂²/∂t² E = ∂²/∂x² E

∂²/∂t² B = ∂²/∂x² B

Schrödinger equation

oscillation equation    ḟ = 𝑖f      fₜ  = 𝑖f
Schrödinger equation    ḟ = 𝑖∆f     fₜ  = 𝑖(fₓₓ + f₝₝ + f₞₞)

heat equation額外乘上虛數𝑖。