引言

真實世界的物理現象，物理學家習慣寫成函數方程式。想要用電腦模擬真實世界，設計函數方程式、解函數方程式是必備技能。

比方說，記錄物體所在位置。根據人類目前所知，物體不會分身，不會同時出現在兩個位置，符合函數的概念；物體不會瞬移，不會瞬時出現在遙遠位置，符合連續的概念。因此物體所在位置可以表示成一個連續函數x(t)，簡寫成x。

數學家創造了函數、連續，主要是為了符合人類認知。如果影分身之術、飛雷神之術成真，那麼數學家勢必要創造其他數學元件，以符合新認知。

位置的變化多寡，稱作速度，符合微分的概念。因此物體當前速度可以表示成d/dt x(t)，簡寫成x′。

物理課教過等速運動。物體速度是40，可以列出等式x′ = 40。大家把40視作一個函數，而非一個實數。

速度也可以忽快忽慢。自訂速度v，可以列出等式x′ = v。

速度的變化多寡，稱作加速度，符合二次微分的概念。因此物體當前加速度可以表示成 d²/dt² x(t)，簡寫成x″。

物理課教過等加速度運動。自由落體，加速度是重力加速度g，g是一個常數約9.8，可以列出等式x″ = -g。如果又有空氣阻力f，得到加速度a = f/m，可以列出等式x″ = -g + f/m。

加速度也可以不斷變化。彗星撞地球，加速度是引力加速度g = -Gm₁m₂ / r²，G是萬有引力常數約6.7e-11，m₁ m₂是質量，r是距離。地心座標定成0，可以列出等式x″ = -Gm₁m₂ / x²。

我們的目標就是解x，知道物體的所在位置。

functional equation

方程式：已知數、未知數全是實數。實數運算有加減乘除。

x² + 2xy + 2y = 1

函數方程式：已知數、未知數全是函數。函數運算有加減乘除微積代入複合。

                d            g(x+2)
∫ f(x) dx + 2 = —— f(g(x)) + ——————
                dx            f(x)

函數方程式當中的實數，其實是函數，稱作「常數函數」。

functional

函數：輸入、輸出全是實數。

f(x, y) = x² + 2xy + 2y - 1

泛函數：輸入、輸出全是函數。函數的函數。

                                d            g(x+2)
L(f(x), g(x)) = ∫ f(x) dx + 2 - —— f(g(x)) - ——————
                                dx            f(x)

我不清楚數學家為何故意讓「泛函數（名詞）」跟「函數的（形容詞）」撞名。

函數方程式的各種要素

一、未知函數：至少一個。
二、函數輸入：至少一個。
三、函數輸出：至少一個。
四、函數輸入運算：變數變換。
五、函數輸出運算：變換。
六、方程式：也包含函數輸入變數。
七、方程式：至少一道。

1. number of unknown function.

   ∂         ∂
   —— f(t) = —— g(t)
   ∂t        ∂t

2. number of input.

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

3. number of output.

                 ⎡ ∂              ∂            ⎤
                 ⎢ —— Fₓ(t,x,y) + —— F₝(t,x,y) ⎥
   ∂             ⎢ ∂x             ∂y           ⎥
   —— F(t,x,y) = ⎢                             ⎥
   ∂t            ⎢ ∂              ∂            ⎥
                 ⎢ —— F₝(t,x,y) - —— Fₓ(t,x,y) ⎥
                 ⎣ ∂x             ∂y           ⎦

4. transform of input.

   ∂              ∂
   —— f(2t,x+1) = —— f(t+1,2x)
   ∂t             ∂x

5. transform of output.

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

6. equation with input.

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

7. number of equations.

   ⎧ ∂             ∂             ∂
   ⎪ —— f(t,x,y) = —— f(t,x,y) + —— f(t,x,y)
   ⎪ ∂t            ∂x            ∂y
   ⎨ 
   ⎪ ∂
   ⎪ —— f(t,x,y) = f(t,x,y)
   ⎩ ∂t

symbolic solution / numeric solution

方程式的解，區分為符號解、數值解。

符號解外觀宛如公式、宛如定律，深受數學家喜愛。

操作科學計算軟體、查閱工程數學書籍，推導符號解。

solve ax² + bx + c = 0   |  solve f(x) + f′(x) = 0
          _________      |
    -b ± V b² - 4ac      |          -x
x = ———————————————      |  f(x) = e
          2a             |

數值解外觀宛如陣列、宛如資料，深受計算學家喜愛。

操作數值模擬軟體、查閱數值分析書籍，演算數值解。

solve 3x² + 12x + 8 = 0  |  solve f(x) + f′(x) = 0
x = -0.7790              |  f(x) = [7.4, 2.7, 1.0, 0.4, 0.1]

古代沒有計算機，古人只找符號解。導致數學課本只教符號解，不教數值解。大部分的方程式，沒人知道符號解，只好改找數值解。

符號解計算工具如Symbolab。數值解計算工具如sundials。

符號解的詳細分類

符號組成式子，式子細分為代數式、閉合式、解析式、……。大家口語上不說符號解，而是說代數解、閉合解、解析解、……。

algebraic solution：只有四則運算。
closed-form solution：包含以上項目，而且多了三角函數。
analytic solution：包含以上項目，而且多了無窮級數。

許多方程式已被證明不存在代數解、閉合解、解析解、……。

Kepler's equation：橢圓軌道運動，已知座標位置(x,y)，求時刻t，沒有閉合解。
Bézier curve：加權平均數曲線，已知座標位置(x,y)，求參數t，沒有閉合解。

延伸閱讀：numerical vs. numeric

numerical：使用到數字的。numeric：本體是數字的。

numerical function：函數的輸入與輸出是數字。

例如地板函數（實數到整數）、絕對值函數（實數到實數）。

numeric function：由數字構成的函數。

例如平方函數，寫成symbolic function是f(x) = x²，寫成numeric function是f([-3,-2,-1,0,1,2,3]) = [9, 4, 1, 0, 1, 4, 9]。

differential equation

「微分方程式」。函數方程式，只有四則運算、微分運算。

           d
f(x) + 2 = —— f(x) + 2 g(x)
           dx

符號解數學公式

微分等於零。解是任意的常數函數。

d
—— f(x) = 0     --->   f(x) = c
dx

微分等於原本函數。解是e^x。原理是e^x微分之後還是e^x。

d                              x
—— f(x) = f(x)  --->   f(x) = e
dx

數學家以此為基礎，讓方程式漸趨複雜，發展各種求解技巧。由於這不是我的專長，就不多提了。

符號解演算法（least squares method）（regression）

函數方程式重新整理成「泛函數等於零」的格式。

函數方程式。
d
—— f(x) = f(x) + 1
dx

等量減法公理，兩邊同減一樣多的東西。（移項）
d
—— f(x) - f(x) - 1 = 0
dx

整理成泛函數的模樣。
          d
L(f(x)) = —— f(x) - f(x) - 1
          dx         

省略括號，簡化符號，比較好讀。
L(f) = f′ - f - 1

推定解是自訂函數，係數待求。

assume f(x) = ax² + bx + c     find a,b,c

自訂函數代入泛函數。通常不是零，通常有誤差。

L(f(x)) = (2a + b) - (ax² + bx + c) - 1

平方誤差盡量小：函數每一處的平方，總和盡量小。

-∞累計到+∞，誤差無限大。因此大家只累計一小段區間。

argmin L(f(x)) dot L(f(x))     點積風格
 a,b,c
argmin ∫ L(f(x))² dx           積分風格
 a,b,c

最小值位於一次微分等於零的地方。對各個係數偏微分。

⎧ ∂/∂a ∫ L(f(x))² dx = 0
⎨ ∂/∂b ∫ L(f(x))² dx = 0
⎩ ∂/∂c ∫ L(f(x))² dx = 0

如果積分很難算，可以偷工減料，對輸入取樣。

let x = -2, -1, 0, 1, 2
⎧ ∂/∂a L(f(-2))² + L(f(-1))² + L(f(0))² L(f(1))² + L(f(2))² = 0
⎨ ∂/∂b L(f(-2))² + L(f(-1))² + L(f(0))² L(f(1))² + L(f(2))² = 0
⎩ ∂/∂c L(f(-2))² + L(f(-1))² + L(f(0))² L(f(1))² + L(f(2))² = 0

整件事情宛如迴歸。一筆資料：輸入值、輸出值（零）。迴歸函數：自訂函數代入泛函數。

data: { (-2,0), (-1,0), (0,0), (1,0), (2,0) }
regression function: L(f(x)) = (2a + b) - (ax² + bx + c) - 1

自訂函數可以推廣成自訂函數網路。函數網路末端接上泛函數，即是代入。大量資料拿來訓練函數網路。

符號解演算法（Galerkin method）（linearization）

推定解是自訂函數的加權總和（基底的線性組合），係數待求。

assume f(x) = w₀ φ₀(x) + w₁ φ₁(x) + w₂ φ₂(x)

再推定泛函數與自訂函數的點積等於零。

可以想成施測。挑選數個函數，進行投影，必須為零。

⎧ L(f(x)) dot φ₀(x) = 0
⎨ L(f(x)) dot φ₁(x) = 0
⎩ L(f(x)) dot φ₂(x) = 0

值得一提的是，一次微分方程式，得到一次方程組，容易求解。

符號解演算法（Parker–Sochacki method）（approximation）

推定解是多項式函數（例如泰勒級數），係數待求。

多項式函數容易微分、容易移項。

數值解演算法（backward Euler method）（implicit Euler method）（dynamic programming）

函數方程式離散化，移項整理，得到遞迴公式，依序填值。

d
—— f(x) = f(x) + 1
dx

(f(x) - f(x - Δx)) / Δx = f(x) + 1     backward difference
(f[i] - f[i-1]   ) / Δx = f[i] + 1     array indexing
f[i] - f[i-1]  = f[i] Δx + Δx
f[i] - f[i] Δx = f[i-1] + Δx
f[i] (1 - Δx)  = f[i-1] + Δx
f[i] = (f[i-1] + Δx) / (1 - Δx)        recurrence

float f[1000]; float Δx = 0.01; float x0 = 0, f0 = 2; // initial condition void backward_euler() { f[0] = f0; // initial condition for (int i=1; i<1000; ++i) f[i] = (f[i-1] + Δx) / (1 - Δx); } void print_solution() { for (int i=0; i<1000; ++i) { float x = x0 + Δx * i; float ans = exp(x) * 3 - 1; // 印出x值、數值解、符號解（事先計算而得） cout << x << ' ' << f[i] << ' ' << ans << '\n'; } }

數值解演算法（forward Euler method）（explicit Euler method）（extrapolation）

難以移項整理的情況下，不得不修改索引值。恰是前向差分。

d
—— f(x) = x sqrt(f(x))
dx

(f(x) - f(x - Δx)) / Δx = x     sqrt(f(x))     backward difference
(f[i] - f[i-1]   ) / Δx = (i-1) sqrt(f[i-1])   array indexing
f[i] = f[i-1] + (i-1) sqrt(f[i-1]) Δx          recurrence

(f(x + Δx) - f(x)) / Δx = x sqrt(f(x))         forward difference
(f[i+1]    - f[i]) / Δx = i sqrt(f[i])         array indexing
f[i+1] = f[i] + Δx i sqrt(f[i])                same recurrence

畫成函數圖形，宛如外插。利用原處斜率，求得下一處函數值。

斜率取自微分方程式的右式，完全不需要移項整理。

f(xₙ₊₁) = f(xₙ) + Δx f′(xₙ)

// d/dx f(x) = x sqrt(f(x)) float f′(float x, float fx) {return x * sqrt(fx);} float f[1000]; float Δx = 1.0; float x0 = 0, f0 = 1; // initial condition void forward_euler() { f[0] = f0; // initial condition for (int i=0; i<1000-1; ++i) f[i+1] = f[i] + Δx * i * sqrt(f[i]); } void print_solution() { for (int i=0; i<1000; ++i) { float x = x0 + Δx * i; float ans = x*x*x*x/16 + x*x/2 + 1; // 印出x值、數值解、符號解（事先計算而得） cout << x << ' ' << f[i] << ' ' << ans << '\n'; } }

數值解演算法（multistep method）（Adams–Bashforth method）

指數函數、三角函數、非線性函數，函數曲線斜率變化極大。只取一處斜率，往往不夠精準。因此取多處斜率，求加權平均數。取先前計算過的斜率，節省計算時間。宛如線性遞迴公式。

斜率數量、斜率地點、權重大小，可以自由設定，但是必須滿足一致性、滿足泰勒近似。一致性留待後面章節介紹。

經典的設定方式是AB2。取兩處斜率：原處、上一處。其兩個權重：3/2、-1/2。

另外還有很多花招，例如AM4、BDF2，這裡就不介紹了。

⎰ f(xₙ₊₁) = f(xₙ) + Δx sₙ
⎱ sₙ = 3/2 f′(xₙ) - 1/2 f′(xₙ₋₁)

f(xₙ) = f(xₙ₋₁) + 3/2 Δx f′(xₙ₋₁) - 1/2 Δx f′(xₙ₋₂)

// d/dx f(x) = x sqrt(f(x)) float f′(float x, float fx) {return x * sqrt(fx);} float f[1000]; float Δx = 1.0; float x0 = 0, f0 = 1; // initial condition float Δf[1000]; void AB2() { f[0] = f0; // initial condition Δf[0] = Δx * f′(x0, f0); // 初始位置斜率（順手乘上Δx） for (int i=1; i<2; ++i) { float x = x0 + Δx * i; // 得到下一處函數值（不用再乘Δx） // 總共兩處斜率：原處、上一處 // 初始位置沒有上一處，姑且用原處代替上一處。 f[i] = f[i-1] + Δf[i-1] * 3/2 - Δf[i-1] * 1/2; // 得到下一處斜率（順手乘上Δx） Δf[i] = Δx * f′(x, f[i]); } for (int i=2; i<1000; ++i) { float x = x0 + Δx * i; // 得到下一處函數值（不用再乘Δx） // 總共兩處斜率：原處、上一處 f[i] = f[i-1] + Δf[i-1] * 3/2 - Δf[i-2] * 1/2; // 得到下一處斜率（順手乘上Δx） Δf[i] = Δx * f′(x, f[i]); } }

數值解演算法（multistage method）（Runge–Kutta method）

方才取左方多處斜率。現在取右方多處斜率。畢竟目的是右方函數值。

右方各處斜率未知，需要一一計算，甚至反覆計算。遞增法改成遞推法。從原處斜率開始，遞推下一處斜率。即時採用最新斜率，更加精準。

斜率數量、斜率地點、權重大小，可以自由設定，但是必須滿足一致性、滿足泰勒近似。一致性留待後面章節介紹。

經典的設定方式是RK4。取四處斜率：原處、半步、半步、一步。其四個權重：1/6、1/3、1/3、1/6。

另外還有很多花招，例如BS23、RKF45，這裡就不介紹了。

// d/dx f(x) = x sqrt(f(x)) float f′(float x, float fx) {return x * sqrt(fx);} float f[1000]; float Δx = 1.0; float x0 = 0, f0 = 1; // initial condition void RK4() { f[0] = f0; for (int i=0; i<1000-1; i++) { float x = x0 + Δx * i; float fx = f[i]; // 遞推得到下一處斜率（順手乘上Δx） // 總共四處斜率：原處、半步、半步、一步 float k1 = Δx * f′(x, fx); float k2 = Δx * f′(x + Δx / 2, fx + k1 / 2); float k3 = Δx * f′(x + Δx / 2, fx + k2 / 2); float k4 = Δx * f′(x + Δx, fx + k3); // 斜率的加權平均數（順手乘上Δx） // 四個權重：1/6、1/3、1/3、1/6 float Δf = (k1 + k2 * 2 + k3 * 2 + k4) / 6; // 得到下一處函數值（不用再乘Δx） f[i+1] = f[i] + Δf; } }

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) basis interpolation (FEM)     基底內插（有限元素法）
  (d) contour (LSM)                 等高線（等高線法）
(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, ⋅));
	}
}

const float t0 = 0.0, tn = 100.0; // interval const float x0 = 0.0, xn = 10.0; // (sampling range) const float Δt = 1.0; // step size const float Δx = 1.0; // (sampling gap) const int ♯t = (tn - t0) / Δt + 1; // number of steps+1 const int ♯x = (xn - x0) / Δx + 1; // (sampling number) float f[♯t][♯x]; // numeric solution void numerical_simulation() { // initialize first frame int n = 0; for (int i=0; i<♯x; i++) { float t = t0 + Δt * n; float x = x0 + Δx * i; f[n][i] = initial_value(t, x); } // draw first frame draw(f[n]); // time-stepping loop for (int n=0; n<♯t-1; n++) { // compute next frame for (int i=0; i<♯x; i++) { float t = t0 + Δt * n; float x = x0 + Δx * i; f[n+1][i] = f[n][i] + F(f[n][i], t, x) * Δt; } // draw next frame draw(f[n+1]); } }

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 value / boundary value：邊界數值。遞迴公式已知數值。
(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 value / boundary value

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

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

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

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

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

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

最難纏的初始值是隨機函數。必須謹慎內插適當數值，以確保數值連續平滑。一種替代方案是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)

float f[1000]; float Δt = 0.01; float t0 = 0, f0 = 2; void exponential_growth() { f[0] = f0; for (int n=0; n<1000-1; ++n) f[n+1] = (f[n] + Δt) / (1 - Δt); // 印出t、f(t)、正確解答 for (int n=0; n<1000; ++n) { float t = t0 + Δt * n; float ans = exp(t) * 3 - 1; cout << t << ' ' << f[n] << ' ' << ans << '\n'; } }

零維空間微分方程式。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

float f[1000]; // 解 float Δt = 0.01; // 取樣間距 float t0 = 0, f0 = 2; // 初始條件 f(t0) = f0 void exponential_growth() { f[0] = f0; for (int n=0; n<1000-1; ++n) f[n+1] = f[n] + (f[n] + 1) * Δt; // 印出t、f(t)、正確解答 for (int n=0; n<1000; ++n) { float t = t0 + Δt * n; float ans = exp(t) * 3 - 1; cout << t << ' ' << f[n] << ' ' << ans << '\n'; } }

零維空間微分方程組。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]

float x[1000]; float y[1000]; float Δt = 0.02; // float t0 = 0; // 初始時刻不重要 float x0 = 10; float y0 = 10; void Lotka_Volterra_equation() { x[0] = x0; y[0] = y0; for (int n=0; n<1000-1; ++n) { x[n+1] = x[n] + (1.1*x[n] - 0.4*x[n]*y[n]) * Δt; y[n+1] = y[n] + (0.1*x[n]*y[n] - 0.4*y[n]) * Δt; } }

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

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

float f[100][10]; float Δt = 0.01, Δx = 0.01; float v = -0.4; // 平流速度 void advection_equation() { // initial condition // 左端有一水源，供水量1。 for (int n=0; n<100; ++n) f[n][0] = 1; // 最初一片虛空 for (int i=1; i<10; ++i) f[0][i] = 0; float C = v * Δt / Δx; for (int n=0; n<100-1; ++n) for (int i=1; i<10; ++i) f[n+1][i] = f[n][i] + C * (f[n][i] - f[n][i-1]); // 印出f(t,x) for (int n=0; n<100; ++n) { for (int i=0; i<10; ++i) cout << f[n][i] << ' '; cout << '\n'; } }

二維空間微分方程式。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 // 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) { 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); } }

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̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖₂ = √ sum (f̂⁽ⁿ⁾(xᵢ) - ĝ⁽ⁿ⁾(xᵢ))² Δx
               = √ sum (f[i] - g[i])² Δx

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

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

float global_error(int n) { float sum = 0; for (int i=0; i<sample_number_x; ++i) { float t = t0 + n * Δt; // 可以挪到迴圈外面 float x = x0 + i * Δx; float diff = f(t, x) - f[n][i]; sum += diff * diff * Δx; // 乘Δx可以挪到迴圈外面 } return sqrt(sum); }

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̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ when n≥0

           with consistency,
           distance of corresponding iteration of
           a symbolic solution and a numeric solution with
           same initial condition is bounded.
           ‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ ≤ ‖f⁽⁰⁾ - f̂⁽⁰⁾‖ when n≥0

convergence: distance of corresponding iteration of
             two numeric solutions with
             different initial conditions leads to zero.
             ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ → 0 when n→∞

             with consistency,
             distance of corresponding iteration of
             a symbolic solution and a numeric solution with
             same initial condition leads to zero.
             ‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖ → 0 when n≥0

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

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

有些文獻又進一步修改定義，只討論動態系統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

數值模擬當中的一致性：當取樣間距Δ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

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

穩定性細分為許多類型。

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₁收斂性

absolute stability

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

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

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

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

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

contractive stability

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

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

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

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

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

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

exponential stability

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

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

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

‖f̂⁽ⁿ⁺¹⁾ - ĝ⁽ⁿ⁺¹⁾‖ ≤ k ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖   time-invariant
k is a constant and 0 ≤ k ≤ 1

當遞迴公式恰是線性遞迴函數，兩種數值解相加減仍是數值解。指數穩定性變成：數值解每回合範數乘上固定倍率，小於等於1。

‖f̂⁽ⁿ⁺¹⁾‖ ≤ k ‖f̂⁽ⁿ⁾‖                   linearity
k is a constant and 0 ≤ k ≤ 1

stability analysis with matrix
（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̂⁽ⁿ⁾

當遞迴公式恰是線性遞迴函數，放大因子的地位宛如先前小節的矩陣。放大因子的平方範數小於等於1，才會滿足收縮穩定性。

       ‖f̂⁽ⁿ⁺¹⁾‖₂
‖a‖₂ = ————————— ≤ ‖A‖₂ ≤ 1           contractive stability
        ‖f̂⁽ⁿ⁾‖₂

當遞迴公式恰是線性遞迴函數，而且，微分方程式恰是零維空間（常微分方程式ODE），此時放大因子與矩陣皆退化成相同數值a。放大因子的絕對值小於等於1，才會滿足收縮穩定性。

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

此時放大因子是複數。複數可以拆解成長度和角度，分別稱作強度magnitude和相位phase，進一步討論強度變化快慢、相位變化快慢。若有興趣請自行深究。

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

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

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

forward difference of t
central difference of x

傅立葉轉換（甚至是拉普拉斯轉換）：

𝓕(f[n+1]) = 𝓕(f[n]) - ½ C (𝓕(f[n]) exp(+𝑖wΔx) - 𝓕(f[n]) exp(-𝑖wΔx))
𝓕(f[n+1]) = (1 - ½ C (exp(+𝑖wΔx) + exp(-𝑖wΔx)) 𝓕(f[n])

放大因子隨之轉換到頻域。此遞迴公式無論如何不可能穩定。

a(w) = 1 - ½ C (exp(𝑖wΔx) - exp(-𝑖wΔx))     amplification factor
a(w) = 1 - C 𝑖 sin(wΔx)                     amplification factor
a numerical scheme is stable iff |a(w)| ≤ 1
this numerical scheme is unstable since |a(w)| > 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(w) = ½ (exp(𝑖wΔx) + exp(-𝑖wΔx)) - ½ C (exp(𝑖wΔx) - exp(-𝑖wΔx))
a(w) = cos(wΔx) - C 𝑖 sin(wΔx)

let |C| ≤ 1 is CFL condition
thus |a(w)|² = cos²(wΔx) + C² sin²(wΔx)
             ≤ cos²(wΔx) + sin²(wΔx)
             = 1
a numerical scheme is stable iff |a(w)| ≤ 1
this numerical scheme is stable under CFL condition |C| ≤ 1

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

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的實部虛部。

【圖片尚待修正】

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。

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

stability analysis

穩定性分析。使用數學工具推導穩定性。業已介紹，茲此總結。

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

(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。

根的數學性質。總之都是單位圓裡面。

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

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

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

專論《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]

numerical scheme

先前介紹了三個項目：

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

以下補充更多詞彙，但是不詳細講解。你們自己通靈吧。

temporal discretization

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

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

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

時間離散化也可以採用高階差分公式。雞肋。

(1) multistep method：多步法。取多處斜率。都是左方斜率。
(2) multistage method：多階段法。求多次斜率。都是右方斜率。
(3) 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

spatial discretization

差分公式，故意增加精度。

泰勒近似，保留更多低次方項。然後想辦法將高階導數替換成低階導數。推導過程在此。

針對中央差分公式，可以使用泰勒近似然後替換高階導數，也可以使用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

差分公式，故意改變項數，並且維持精度。

求得差分公式，可以使用泰勒近似。簡化差分公式並且維持精度，可以進一步使用Padé approximation。

(1) compact finite difference：用更少的項數，達到更高的精度。
(2) degraded finite difference：達到更低的精度。自虐。

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

implicit method / explicit method

(1) 決定時刻。
  (a) implicit method：右式取新值。
  (b) explicit method：右式取舊值。
  (c) implicit–explicit method：有些取新值、有些取舊值。
  (d) Crank–Nicolson method：新值加舊值然後除以二。理都懂然并卵。
(2) 分割步驟。
  (a) operator splitting
  (b) Strang splitting
(3) 更動步驟。
  (a) semi-implicit method (symplectic 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
                            ^^^^
Crank–Nicolson method。新值加舊值然後除以二。
     d/dt f(t) = f(t) + 1
---> (f[n+1] - f[n]) / Δt = (f[n+1] + f[n]) / 2 + 1
                            ^^^^^^^^^^^^^^^^^^^
operator splitting。項次太多，逐步更新。
     d/dt f(t) = f(t) + g(t)
---> 1. (f[n+1] - f[n]) / Δt = f[n+1]
     2. (f[n+1] - f[n]) / Δt = g[n+1]

semi-implicit method。式子太多，安排順序。
     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]

https://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book018.html
https://en.wikipedia.org/wiki/Strang_splitting

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

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

domain discretization

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

(1) mesh-based method：網格法。取樣擇鄰：隨機點／網格／基底內插。
(2) meshfree method：無網格法。隨機點／半徑範圍／基底內插。
(3) fixed mesh refinement：固定網格細緻化。事先指定空間各處網格精細程度。
(4) adaptive mesh refinement：自適應網格細緻化。各處根據誤差調整網格精細程度。

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

numerical iteration

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

(1) domain decomposition method：領域分解法。矩陣求解切成多個區塊輪流計算。
(2) multigrid method：多格法。矩陣求解採用scaling method。
(3) spectral method：頻域法。矩陣求解轉換至頻域計算。

matrix solver

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

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

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

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

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

domain decomposition method

spectral method

頻域法即是傅立葉轉換。FFT。

multigrid method

numerical error

error suppression：誤差抑制。調整數值解，減少數值誤差。

(1) smooth transition：不連續處強行平滑化，避免Runge phenomenon。
(2) numerical smoothing：整體強行平滑化，避免不穩定。
(3) lake at rest：零、負數強行改成machine epsilon，避免歸零。

smooth transition

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

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

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

波浪撞牆反彈demo

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

numerical 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。

differential equation: linear 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

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/diffusion/viscosity。牛頓流體當中，這些物理現象都是對空間二次微分，即是Laplacian。

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

lake at rest

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

numerical scheme

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

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

(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̂⁽ⁿ⁾ ≤ ĝ⁽ⁿ⁾ => f̂⁽ⁿ⁺¹⁾ ≤ ĝ⁽ⁿ⁺¹⁾

[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/

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，畢竟沒有二階偏微分。

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。
純量：單變量。函數輸出是數值、而不是向量。

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

linear advection equation           線性平流方程式
linear transport equation           線性傳輸方程式
u ≔ u , f(u) ≔ au                   u: physical quantity
                                    a: constant velocity

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

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

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

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

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

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

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

solution of scalar conservation law in one-dimensional space

線性守恆律，初始條件是函數，那麼時時刻刻都是函數，形成單值解。例如linear 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》。

線性守恆律。

linear 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)：抱歉我學不會。我選擇放棄。

numerical method

數值模擬，總共四種方法：

(1) finite difference method：有限差分法。使用原本的微分方程式。
(2) finite volume method：有限體積法。改成通量。
(3) finite element method：有限元素法。改成引力。
(4) level set method：等高線法。改成勢。

finite difference method

有限差分法，可以求得守恆律的解析解。

有限差分法是一類數值方案的泛稱。先前的動態系統章節，每一道遞迴公式都是有限差分法。

domain discretization: uniform grid
codmain discretization: function point value
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|

Courant–Friedrichs–Lewy condition

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

守恆律，一旦滿足收縮穩定性，也會連帶滿足「CFL條件」。

for linear scheme of conservation law,
contractive stability => CFL condition

【待補範例】

線性守恆律，CFL condition是取樣間距Δt與Δx、倍率a，三者構成的不等式。

linear conservation law
du     du
—— + a —— = 0
dt     dx

CFL condition for linear conservation law
a Δt
———— ≤ 1
 Δx

一般的守恆律，CFL condition修改成數值模擬過程當中的倍率最大值。硬撐過去。

conservation law
du        du
—— + a(u) —— = 0
dt        dx

CFL condition for conservation law
aₘₐₓ Δt
——————— ≤ 1    where aₘₐₓ is the maximum value of a(u)
   Δx          in numerical simulation

順帶一提，當守恆律是方程組，而且矩陣A(u)可以對角化，透過特徵分解和座標轉換，每道等式得以各自討論CFL condition。

順帶一提，守恆律的某些衍生版本，也有類似性質。舉例來說，波動方程式是a Δt / Δx ≤ 1。熱傳導方程式是2 a Δt / (Δx)² ≤ 1。這裡就不細講了，推導過程請見講義：

倍率a可以視作移動速度。教科書習慣將a改寫成v。

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

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

【待補圖片】

CFL condition原始定義：characteristic curve當中，符號解的domain of dependence包含於數值解的domain of dependence。

線性守恆律，符號解的domain of dependence是一條直線，斜率1 / a。數值解必須遵循這條直線，斜率Δt / Δx。因此a Δt / Δx = 1。【尚待確認】

一般的守恆律，符號解的domain of dependence是一塊範圍。因此aₘₐₓ Δt / Δx ≤ 1。

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

CFL condition是必要條件：滿足收縮穩定性，導致滿足CFL condition。CFL condition不是充分條件：滿足CFL condition，可能穩定、可能不穩定。

CFL condition的主要功用是快篩試劑。逆否命題：不滿足CFL condition，導致不滿足收縮穩定性。

for linear scheme of continuity equation,
contractive stability ⇍ CFL condition

https://scicomp.stackexchange.com/questions/2927/
https://scicomp.stackexchange.com/questions/25398/

換個說法。各種數值方案，經過推導，得到各種穩定性條件。CFL condition只不過是數學式子最漂亮的那一個穩定性條件，不等式右側恰是1。大部分穩定性條件更加嚴格，不等式右側低於1。

至少要滿足CFL condition，才有可能滿足收縮穩定性。如果連CFL condition都達不到，那麼不可能滿足收縮穩定性。

例如後向差分（一階精度）是1。
例如中央差分（四階精度）是0.728...。

Courant number

CFL condition的左式，代入a、Δt、Δx的實際數值，求得左式數值，此數值稱作「Courant number」或「CFL number」。

CFL number
    aₘₐₓ Δt
C = ———————
       Δx

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

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

一、設計數值方案。
二、設定初始條件、邊界條件。
三、CFL number調整為1，以便滿足CFL condition。
四、逐步降低CFL 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 ≤ ΔtCFL,
  then explicit multistage method reaches contractive stability
  under CFL condition Δt ≤ min (aᵢⱼ/bᵢⱼ) ΔtCFL
                           ⁱʲ
  when a and b are non-negative
  and a satisfies convex combination property.

proof:

  ‖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ᵢⱼ) ≤ ΔtCFL

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

「保強穩定性」內容很雲。CFL condition只是必要條件，此處卻以充分條件作為前提。前提是一件未必成立的事情。

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

finite volume method

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

有限體積法是一類數值方案的泛稱：未知函數切成小段，每一段的定積分的平均數，作為遞迴公式的變數。

domain discretization: uniform grid
codmain discretization: numerical flux
                       (average of definite integral)
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的正解。

順帶一提，當守恆律是linear 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:
  ũᵢ₊₁⸝₂(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ᵢ ûᵢ⁽ⁿ⁾ = sumᵢ ûᵢ⁽ⁿ⁺¹⁾

quantity average
      1 ⌠ xᵢ₊₁⸝₂
ûᵢ = —— ⎮        u(tₙ,x) dx   when Δx→0
     Δx ⌡ xᵢ₋₁⸝₂

consistency
f̂(ûᵢ,ûᵢ,...,ûᵢ) = f(ûᵢ)

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

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

我還不確定要不要介紹

artificial viscosity：人工黏性。姑且先放一篇Lax的文章，撐一下場面。

Cauchy problem：微分方程式，指定邊界。

viscous Cauchy problem：追加二次微分項。

entropy conservative scheme：僅等號成立。不收縮。

entropy conservative scheme [Tadmor 1987]
η(uᵢ⁽ⁿ⁺¹⁾) = η(uᵢ⁽ⁿ⁾) - λ (qᵢ₊₁⸝₂⁽ⁿ⁾ - qᵢ₋₁⸝₂⁽ⁿ⁾)
where entropy pair (η,q)
      q̃(u,u) = q(u)
      qᵢ₊₁⸝₂ = q̃(uᵢ,uᵢ₊₁)

entropy conservation
(η′ʟ - η′ʀ) flux(uʟ, uʀ) = qʟ - qʀ

https://www.ams.org/journals/mcom/1987-49-179/S0025-5718-1987-0890255-3/

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: triangular mesh
codmain discretization: shape function
                       (basis interpolation by shape functions)
operator discretization: chain rule

定義域離散化：方格／網格／粒子。
函數離散化：即是線性代數的「垂直投影到子空間」。
　　　　　　子空間的基底是shape function。
　　　　　　保證存在唯一解。但是我還不知道截斷誤差是多少。
微分運算離散化：其實就是基底函數套用微分連鎖律，求得數學式子。

consistency：符號解與數值解的一致性。

專著《Partial Differential Equations》。

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

shape function：基底函數。我不知道如何提升精度。

講義《Computational Magnetohydrodynamics》。

Galerkin finite element method：微分方程式與基底函數的點積都是零。
pseudo-spectral method：基底函數是正交多項式，例如弦波。

numerical scheme：數值解的數學性質、其對應的數值方案，目前沒有定論。這是最接近此主題的台灣學校課程。

discontinuous Galerkin method (DG method)
Nitsche's method

spectral difference method
spectral volume method

https://users.oden.utexas.edu/~leszek/classes.html
https://users.oden.utexas.edu/~leszek/classes/EM394H/book.pdf

finite volume method：先前小節的數值方案可以改寫成基底內插函數。

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追加basis interpolation
Cockburn、舒其望

arbitrary high order using derivatives (ADER-DG method)
Godunov's method追加high order basis interpolation
Toro

high-resolution scheme
即是higher order numerical scheme
大於一階精度的數值方案們。例如MUSCL、DG method。
Harten

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

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

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

微分方程式的分類：多變數函數

微分方程式分為兩類。一、常微分方程式：對同一種變數進行微分。二、偏微分方程式：對多種變數進行微分。

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 = ḟ + 3f̈ + 2g         微分簡寫成 ḟ 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 = Dₓf + 3Dₓ₝f + 2g    微分簡寫成 Dₓf Dₓ₝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 指定解的邊界的位置、函數值、導數值

函數值：確保解的零次項（常數項）、一次項、二次項、……。
導數值：確保解的一次項、二次項、……。

解可以畫成圖形。

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

微分方程式，化作代數方程式。

經典的微分方程式：Euler's Identity

引入了歐拉恆等式，作為微分方程式的解。

exponential equation  ḟ = f
oscillation equation  ḟ = 𝑖 f

exponential equation

「指數方程式」。

ḟ = f，微分之後還是一樣。符號解是指數函數。微分不動點。

ḟ = γf ，添上倍率γ。γ>1是成長率，γ<1是衰減率。

oscillation equation

「振盪方程式」。指數方程式，令實部為零。

ḟ = 𝑖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₞₞
Wave equation       f̈ = ∆f     fₜₜ = fₓₓ + f₝₝ + f₞₞

注意到，數學家使用牛頓微分符號、梯度／散度／旋度／梯散符號，微分對象涵蓋所有變數。物理學家異於常人，牛頓微分符號只針對時間變數，梯度／散度／旋度／梯散符號只針對空間變數。

物理學家另行發明梯度運算spacetime gradient、梯散運算d'Alembertian，以便針對所有變數。食飽換枵。

Helmholtz equation

「亥姆霍玆方程式」。

f = ∆f，函數等於位勢差。

f = -k²∆f，添上波數k。

heat equation

「熱傳導方程式」。硬生生多出傳導二字。

ḟ = ∆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型鋼、平板

經典的微分方程式：transport

引入了散度運算。

continuity equation  ψ̇ + ∇∙(ψu) = 0     ψₜ + (ψu)ₓ + (ψu)₝ = 0
advection equation   ψ̇ + u∙∇ψ = 0       ψₜ + u₁ψₓ + u₂ψ₝ = 0
incompressibility    ∇∙u = 0            (u₁)ₓ + (u₂)₝ = 0

continuity equation

「連續方程式」。

ψ̇ + ∇∙ψ = 0，時間變化量等於空間轉移量。

ψ̇ + ∇∙(ψu) = 0，添上每一處的轉移速度u。

advection equation

「平流方程式」。

ψ̇ + ∇ψ = 0，時間變化量等於相鄰空間差距。

ψ̇ + u∙∇ψ = 0，添上每一處的轉移速度u。

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

連續方程式可以拆成兩項：平流與壓縮。

微分的乘法公式：前微後不微、前不微後微。

ψ̇ + ∇∙(ψu) = 0
ψ̇ + u∙(∇ψ) + ψ(∇∙u) = 0
    ^^^^^^     ^^^^
    advection  compression

∂        ∂        ∂           
—— ψ + ( —— ψu₁ + —— ψu₂ ) = 0   u₁是X速度
∂t       ∂x       ∂y             u₝是Y速度

∂           ∂         ∂            ∂       ∂          
—— ψ + ( u₁ —— ψ + u₂ —— ψ ) + ψ ( —— u₁ + —— u₂ ) = 0
∂t          ∂x        ∂y           ∂x      ∂y         
       ^^^^^^^^^^^^^^^^^^^^^     ^^^^^^^^^^^^^^^^^
       advection                 compression

平流方程式只有時間變化項、平流項。

有人簡寫成大寫D，稱作「物質導數」。

∂           ∂         ∂                   D       
—— ψ + ( u₁ —— ψ + u₂ —— ψ ) = 0   --->   —— ψ = 0
∂t          ∂x        ∂y                  Dt      

平流方程式追加∇∙ψ = 0，則滿足連續方程式，物理量守恆。

平流方程式不追加∇∙ψ = 0，則不一定滿足連續方程式，物理量不一定守恆。

符號解：f(x,t) = e^-(x+t)²。【尚待確認】

incompressibility

「不可壓縮」。

∇∙u = 0。出入速度，總和相等。流動順暢，不疾不徐。

經典的微分方程式：flow

引入了梯度運算，形成了勢、力。

continuity equation    ∂/∂t ρu + ∇∙(ρu⊗u) = 0
Euler equation         ∂/∂t ρu + ∇∙(ρu⊗u) + ρg + ∇p = 0
Cauchy equation        ∂/∂t ρu + ∇∙(ρu⊗u) + ρg + ∇∙τ = 0
Navier–Stokes equation ∂/∂t ρu + ∇∙(ρu⊗u) + ρg + ∇p + ∇∙τ = 0

hydrodynamics equations

continuity equation for momentum conservation

「連續方程式之動量守恆」。物理量是動量密度ψ = ρu。

質量m。體積V。密度ρ = m/V。速度u。動量mu。動量密度ρu = mu/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的特例：無旋(有勢)無散(不可壓縮)純粹諧場
   ∂ψ/∂t + 1/2 (∂u²/∂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。

記得介紹一下diffusion和dispersion。

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

符號解自遞迴。

reaction–diffusion–advection equation
Turing pattern
Gray–Scott model

reaction–diffusion–advection equation

「反應－擴散－平流方程式」。各式各樣的參數，各式各樣的造型，造就大自然。

引入了合體技的概念。

advection equation                      ḟ = ∇∙f
diffusion equation                      ḟ = ∆f
advection–diffusion equation            ḟ = ∆f + ∇∙f
reaction–diffusion equation             ḟ = ∆f + R(f)
reaction–diffusion–advection equation   ḟ = ∇∙f + ∆f + R(f)

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 = mẍ + cẋ + 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

⎰ ẋ = ax - bxy
⎱ ẏ = 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
Schrödinger equation    ḟ = 𝑖 ∆f     fₜ  = 𝑖(fₓₓ + f₝₝ + f₞₞)

heat equation額外乘上虛數𝑖。