引言

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

比方說，記錄物體所在位置。根據人類目前所知，物體不會分身，不會同時出現在兩個位置，符合函數的概念；物體不會瞬移，不會瞬時出現在遙遠位置，符合連續的概念。因此物體所在位置可以表示成一個連續函數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. change of variable.

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

5. transformation.

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

6. transformation by input. (if scalar-valued)

   ∂                  ∂
   —— 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）

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

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

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

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

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

重新整理方程式，以便形成動畫。左式只有一項：未知函數對時間微分。

∂                    ∂        
—— 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² + ...

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

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

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

symbolic solution

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

numeric solution

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

大家以動態系統的標準式求得數值解。給定當前時刻的數值解，求得下個時刻的數值解。這個手法稱作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()
{
	// initialize first frame
	f(t₀, ⋅) := initial_value;
	// draw first frame
	draw(f(t₀, ⋅));
	
	// time-stepping loop
	for (int t = t₀; t < final_time; t += Δt)
	{
		// compute next frame
		f(t+Δt, ⋅) := f(t, ⋅) + F(i, t, ⋅) * Δt;
		// draw next frame
		draw(f(t+Δt, ⋅));
	}
}

float f[100][10]; // numeric solution float Δx = 1, Δt = 1; float x0 = 0, t0 = 0; void numerical_simulation() { // initialize first frame int n = 0; for (int i=0; i<10; i++) { float t = t0; float x = x0 + Δx * j; f[n][i] = initial_value(t, x); } // draw first frame draw(f[n]); // time-stepping loop for (int n=0; n<100-1; n++) { // compute next frame for (int i=0; i<10; 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):

f(t+Δt,x) - f(t,x)
—————————————————— = (3/2) F(f(t   ,x), t   , x)
        Δt
                   - (1/2) F(f(t-Δt,x), t-Δt, x)

multistage method (RK4):

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

改寫成陣列。

Euler method:

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

multistep method (AB2):

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

multistage method (RK4):

⎧ F₁ = F(f[n][i]           , t,       x)
⎪ F₂ = F(f[n][i] + F₁ ⋅ ½Δt, t + ½Δt, x)
⎨ F₃ = F(f[n][i] + F₂ ⋅ ½Δt, t + ½Δt, x)
⎪ F₄ = F(f[n][i] + F₃ ⋅ 1Δt, t + 1Δt, x)
⎩ f[n+1][i] = f[n][i] + F₄ ⋅ Δt

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

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。

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

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

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

numerical iteration

數值遞推。套用遞迴公式，逐步求得數值解。

計算過程，建立表格，依序填值，得到數值解。

時間為主，空間為輔，使得表格對應動畫。

boundary

邊界。計算範圍。細分為時間範圍、空間範圍。

我們只關心一塊範圍，而不是從負無限大到正無限大。

邊界。輪廓形狀。例如矩形、圓形、任意封閉區域。

其中最棒的形狀是一維區間、二維矩形、三維立方體，恰好符合程式語言的陣列格式。

f(t)       5 ≤ t ≤ 10      time is an 1D interval
                           space is a 0D point
                           => boundary is two endpoints

f(t,x)     5 ≤ t ≤ 10      time is an 1D interval
           -10 ≤ x ≤ +10   space is an 1D interval
                           => boundary is a rectangle 

f(t,x,y)   5 ≤ t ≤ 10      time is an 1D interval
           x² + y² ≤ 10    space is a 2D disk
                           => boundary is a cylinder

initial value / boundary value

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

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

動態系統當中，問題被分成初始值問題、邊界值問題，邊界值被分成不是指定數值、是指定數值。

初始值問題較簡單，演算法如前所述。邊界值問題較困難，演算法更加複雜，此處省略。

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

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

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

其中最難纏的初始值是隨機函數。必須謹慎內插適當數值，以確保邊界內部數值連續。一種替代方案是【尚待確認】：函數值是隨機弦波疊加。

邊界值（不是指定數值）視實際問題而定。例如無界函數、對稱函數、週期函數。

(1) open boundary condition      開放邊界條件：無限空間
(2) symmetric boundary condition 對稱邊界條件：牆面反彈
(3) periodic boundary condition  週期邊界條件：無限循環

其中最難纏的邊界值是無界函數。必須謹慎外插適當數值，以確保邊界內緣數值不受影響。一種替代方案是sponge layer：海綿層置於邊界外緣，海綿層的函數值快速衰減至零。

邊界值（是指定數值）視實際問題而定。例如函數值、函數導數值。

(1) Dirichlet boundary condition 函數值
(2) Neumann boundary condition   函數導數值（的邊界法線分量）
(3) Cauchy boundary condition    上述兩者
(4) Robin boundary condition     上述兩者的加權總和

這些邊界值都相當難纏。必須轉換成積分方程式，而這不屬於本章範疇。一種替代方案是shooting method：亂槍打鳥、試誤法。邊界值問題簡化為初始值問題。嘗試各種初始值（發射地點），各自進行數值遞推（行進軌道），檢查邊界值是否恰好符合指定數值（恰好射中邊界值）。

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

程式碼

現在馬上來一段程式碼，讓各位實際觀察numerical 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; } }

高階時間微分方程式，化作一階時間微分方程組。

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

省略

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

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。

error metric

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

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

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

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

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

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

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

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

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

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

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
    algorithm: Cauchy condensation test
    convergence: monotone convergence theorem
                 (if monotone, then bounded = convergent)

(2) problem: equation solving
    algorithm: fixed point iteration
    convergence: Banach fixed-point theorem
                 (if contractive, then convergent)

(3) problem: differential equation solving
    algorithm: Euler method
    convergence: fundamental theorem of numerical analysis
                 (if consistent, then stable = convergent)

問題：級數收斂
演算法：柯西收縮測試
收斂性：單調收斂定理（單調時，受限即收斂。）

問題：方程式求解
演算法：不動點遞推法
收斂性：巴拿赫不動點定理（收縮則收歛。）

問題：微分方程式求解
演算法：歐拉法
收斂性：數值分析基本定理（一致時，穩定即收斂。）

fundamental theorem of numerical analysis

「數值分析基本定理」。

if consistent, then stable = convergent

拆成兩句話。第一句話是本文主角。第二句話顯而易見。

(1) consistent & stable => consistent & convergent
(2) consistent & convergent => consistent & stable   (trivial)

這三個詞彙的中文翻譯：

adjective   形容詞丨 noun        名詞　
------------一一一丨-------------一一一
consistent  一致的丨 consistency 一致性
stable      穩定的丨 stability   穩定性
convergent  收斂的丨 convergence 收斂性

這三個詞彙的原始意義：

consistency: 1. theorems have no contradiction.
             2. an equation have a unique solution.
             3. symbol and numeral are approximately equal.

stability: two sequence are bounded by a certain distance
           that is depending on initial distance
           |aₙ - bₙ| < ε when |a₀ - b₀| < δ
           where ε and δ are constants

convergence: a sequence leads to a constant.
             aₙ → c when n → ∞
             where c is a constant

一致性：定理無矛盾。等式有唯一解。符號與數值相符。表裡如一。
穩定性：兩串數列，對應項差距範圍受限。牽繩遛狗。
收斂性：一串數列，數字趨近定值。閒逛到家。

這三個詞彙在數值模擬當中的用法：

（我沒有照抄文獻裡的句子。我不確定我的理解是否正確。）

differential equation
L(f(t,x)) = 0

symbolic solution
f(t,x)

consistency: difference of the functional equation and
             its discretization approaches zero.
             L(f(t,x)) - L̂(f(t,x)) → 0 when Δt→0, Δx→0
             where L(f(t,x)) = (d/dt f(t,x)) - F(f(t,x))

             big-O notation is used somehow.
             L(f(t,x)) = L̂(f(t,x)) + O((Δt)ᵖ + (Δx)𐞥)
             where p>0, q>0

stability: distance of corresponding iteration of
           two numeric solutions with
           different initial conditions is bounded.
           ‖f̂⁽ⁿ⁾ - ĝ⁽ⁿ⁾‖ ≤ ‖f̂⁽⁰⁾ - ĝ⁽⁰⁾‖ when n≥0

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

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

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

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

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

大部分文獻又進一步修改定義，將f(t,x)改成d/dt f(t,x)再改成微分方程式的右式F(f(t,x))。我推測原因是避免在Δt→0的情況下乘以除以Δt，並且確保時間變化空間變化可以分離成左式右式。

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

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

symbolic solution
f(t,x)

consistency
F(f(t,x)) - F̂(f(t,x)) → 0 when Δt→0, Δx→0

stability with consistency
‖F(f⁽ⁿ⁾) - F(f̂⁽ⁿ⁾)‖ ≤ ‖f⁽⁰⁾ - f̂⁽⁰⁾‖ when n≥0

convergence with consistency
‖F(f⁽ⁿ⁾) - F(f̂⁽ⁿ⁾)‖ → 0 when n≥0

時間變化空間變化可以分離成左式右式，遞迴公式不受時間影響，每回合一律使用相同泛函數F̂(f(t,x))，事情比較簡單。時間變化空間變化無法分離成左式右式，遞迴公式受時間影響，每回合都要微調泛函數F̂ₜ(f(t,x))，事情更加棘手，我沒有研究。

數值分析基本定理的證明，我還沒有學會。等我學完泛函分析不知何年何月。關鍵字也許是Banach space、uniform boundedness principle、dominated convergence theorem。

https://math.stackexchange.com/questions/4460567/
https://www.colorado.edu/amath/sites/default/files/attached-files/lax_equivalence_theorem.pdf
https://users.soe.ucsc.edu/~hongwang/AM213B/Notes/Lecture02.pdf
https://users.soe.ucsc.edu/~hongwang/AM213B/Notes/Lecture03.pdf

補充說明一下。不動點定理當中，平緩函數只能是小於。數值分析基本定理當中，穩定性可以追加等於。第零回合，令數值解趨近符號解、令初始距離趨近零；往後回合，即便距離不變也無妨。一開始就夾緊，往後保持夾緊，依然收斂。

補充說明一下。數值分析基本定理當中，只有設置初始條件，沒有設置邊界條件。符號解和數值解，左右延伸至無限遠。另外也有人採用週期邊界條件，避免討論無限遠。

總之，當一致性和穩定性同時成立，就能保證數值收斂。

接下來分別介紹一致性和穩定性。

consistency與consistency error

「一致性」。符號轉數值，符號與數值相符。分為兩點討論。

一、數字大小相符。
二、運算結果相符。加減乘除微積求值代入。

「一致性誤差」。符號轉數值導致的誤差。分為兩點討論。

round-off error：四捨五入誤差。實數四則運算從符號轉數值導致的誤差。
truncation error：截斷誤差。函數微分運算從符號轉數值導致的誤差。

實數四則運算從符號轉數值，一般採用浮點數運算。四捨五入誤差，一般是指浮點數運算誤差。請見本站文件「bit」。

函數微分運算從符號轉數值，一般採用泰勒近似。截斷誤差，一般是指泰勒近似誤差。請見本站文件「approximation」。

floating-point arithmetic與round-off error

「浮點數運算」。離散化實數之四則運算。

實數資料結構採用浮點數，實施四則運算。

「四捨五入誤差」。浮點數運算總是消滅多餘的低位數。

本節不是本文重點。細節請自行上網搜尋講義。

finite difference與truncation error

「有限差分」。離散化實函數之微分運算。

實函數資料結構採用regular discretization，或者其他資料結構，實施微分運算。

「截斷誤差」。泰勒近似總是消滅多餘的高次方項。

一次微分的差分公式。細分成左邊、右邊、中央三種版本。

Taylor approximation
(1) f(x+Δx) = f(x) + f′(x)Δx + ½f″(x)(Δx)² + ⅙f‴(x)(Δx)³ + ...
(2) f(x-Δx) = f(x) - f′(x)Δx + ½f″(x)(Δx)² - ⅙f‴(x)(Δx)³ + ...

first-order backward difference
f(x-Δx) = f(x) - f′(x)Δx + O((Δx)²)
f(x) - f(x-Δx) = f′(x)Δx + O((Δx)²)
(f(x) - f(x-Δx)) / Δx = f′(x) + O(Δx)
f′(x) = (f(x) - f(x-Δx)) / Δx + O(Δx)
f′[i] = (f[i] - f[i-1]) / Δx + O(Δx)

first-order forward difference
f(x+Δx) = f(x) + f′(x)Δx + O((Δx)²)
f(x+Δx) - f(x) = f′(x)Δx + O((Δx)²)
(f(x+Δx) - f(x)) / Δx = f′(x) + O(Δx)
f′(x) = (f(x+Δx) - f(x)) / Δx + O(Δx)
f′[i] = (f[i+1] - f[i]) / Δx + O(Δx)

second-order central difference
⎰ f(x+Δx) = f(x) + f′(x)Δx + ½f″(x)(Δx)² + O((Δx)³)
⎱ f(x-Δx) = f(x) - f′(x)Δx + ½f″(x)(Δx)² + O((Δx)³)
f(x+Δx) - f(x-Δx) = 2 f′(x)Δx + O((Δx)³)
(f(x+Δx) - f(x-Δx)) / (2 Δx) = f′(x) + O((Δx)²)
f′(x) = (f(x+Δx) - f(x-Δx)) / (2 Δx) + O((Δx)²)
f′[x] = (f[i+1] - f[i-1]) / (2 Δx) + O((Δx)²)

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

Δx趨近於零、Δx遠小於一。Δx次方越高，截斷誤差越小。泰勒展開高次方項消滅越少，截斷誤差越小。

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

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

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

二次微分的差分公式。

second-order central difference of second derivative
⎰ f(x+Δx) = f(x) + f′(x)Δx + ½f″(x)(Δx)² + ⅙f‴(x)(Δx)³ + O((Δx)⁴)
⎱ f(x-Δx) = f(x) - f′(x)Δx + ½f″(x)(Δx)² - ⅙f‴(x)(Δx)³ + O((Δx)⁴)
f(x+Δx) + f(x-Δx) = 2 f(x) + f″(x)(Δx)² + O((Δx)⁴)
f(x+Δx) - 2 f(x) + f(x-Δx) = f″(x)(Δx)² + O((Δx)⁴)
(f(x+Δx) - 2 f(x) + f(x-Δx)) / (Δx)² = f″(x) + O((Δx)²)
f″(x) = (f(x+Δx) - 2 f(x) + f(x-Δx)) / (Δx)² + O((Δx)²)
f″[i] = (f[i+1] - 2 f[i] + f[i-1]) / (Δx)² + O((Δx)²)

一次微分的差分公式，增加點數以增加精度。我沒有研究。

4th-order central difference
f′[x] = (- f[i+2] + 8 f[i+1] - 8 f[i-1] + f[i-2]) / (12 Δx)
      + O((Δx)⁴)

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

想要迅速獲得中央差分公式，除了泰勒展開以外，也可以使用Lagrange interpolation然後微分。想要簡化差分公式並且維持精度，可以使用Padé approximation。我沒有研究。

consistency of numerical simulation

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

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

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

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

微分方程式重新整理成泛函數的模樣，才能相減得到差距。

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

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

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

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

truncation error:
L(f(t)) - L̂(f(t)) = O((Δt)ᵖ)

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

微分方程式的各個微分運算的截斷誤差，組成了局部誤差。截斷誤差趨近零，那麼局部誤差也趨近零。想讓截斷誤差趨近零，一種方式是符號解f(t)採用平滑函數。

順帶一提，截斷誤差同時考慮時間變數與空間變數的步進，局部誤差只考慮時間變數的步進。有些文獻誤將截斷誤差與局部誤差視為相同，誤將O((Δt)ᵖ + (Δx)𐞥)與‖f⁽¹⁾ - f̂⁽¹⁾‖視為相同。

時間變數的步進，最簡單的方式是採用一階右邊差分，p = 1。有些文獻因而省略時間變數的步進，並將截斷誤差與局部誤差視為相同，將O((Δx)𐞥)與‖f⁽¹⁾ - f̂⁽¹⁾‖視為相同。

stability

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

穩定性細分為許多類型。

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

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

stability of numerical simulation

數值模擬當中的穩定性：數值解距離受限。

穩定性細分為許多類型。

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

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

數值分析基本定理，理論上是受限穩定性，實務上是收縮穩定性。大家難以判斷受限穩定性，大家換成更強更嚴格的收縮穩定性。

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

contractive stability

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

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

當遞迴公式不因時刻而變，收縮穩定性變成：數值解每回合距離變小或不變。

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

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

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

exponential stability

「指數穩定性」。數值解距離指數衰減。

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

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

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

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

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

stability analysis with matrix
（Lax–Richtmyer stability analysis）

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

實際範例請見本站文件「linear recurrence」。

‖f̂⁽ⁿ⁺¹⁾‖ = ‖Af̂⁽ⁿ⁾‖ ≤ ‖f̂⁽ⁿ⁾‖           contractive stability

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

‖f̂⁽ⁿ⁺¹⁾‖ = ‖Af̂⁽ⁿ⁾‖ ≤ k ‖f̂⁽ⁿ⁾‖         exponential stability
k is a constant and 0 ≤ k ≤ 1

另一方面，根據範數的性質，得到不等式。

‖Af̂⁽ⁿ⁾‖ ≤ ‖A‖ ‖f̂⁽ⁿ⁾‖                  sub-multiplicative

根據不等式，收縮穩定性等價於矩陣範數小於等於1。

矩陣範數視作倍率k。收縮穩定性等價於指數穩定性！

不斷乘以一個小於等於1的倍率，顯然指數衰減。倍率從數字換成矩陣，依然指數衰減。就這樣。

‖f̂⁽ⁿ⁺¹⁾‖ = ‖Af̂⁽ⁿ⁾‖ ≤ ‖A‖ ‖f̂⁽ⁿ⁾‖       sub-multiplicative
‖f̂⁽ⁿ⁺¹⁾‖ ≤ ‖f̂⁽ⁿ⁾‖                     contractive stability
‖A‖ ≤ 1                               exponential stability

常見的範數類型是L₂-norm。矩陣範數採用譜範數：特徵值絕對值最大者。向量範數採用平方範數：元素平方和開根號。

‖f̂⁽ⁿ⁺¹⁾‖₂ = ‖Af̂⁽ⁿ⁾‖₂ ≤ ‖A‖₂ ‖f̂⁽ⁿ⁾‖₂   sub-multiplicative
‖f̂⁽ⁿ⁺¹⁾‖₂ ≤ ‖f̂⁽ⁿ⁾‖₂                   contractive stability
‖A‖₂ ≤ 1                              exponential stability

矩陣特徵值絕對值最大者小於等於1。矩陣特徵值皆介於±1（若是實數）。

‖A‖₂ = max(abs(eigval(A))) ≤ 1        spectral norm
0 ≤ abs(eigval(A)) ≤ 1                inside the unit circle
-1 ≤ eigval(A) ≤ 1                    if eigval(A) are real

想要檢查指數穩定性，只需將線性遞迴公式改寫成矩陣，然後檢查所有特徵值絕對值皆小於等於1。

0 ≤ abs(eigval(A)) ≤ 1                exponential stability

順帶一提，指數穩定性顯然滿足收斂性，畢竟指數衰減至零。大可不必煩惱數值分析基本定理。

absolute stability

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

穩定性細分許多類型。

zero-stability：零穩定性。當Δt與Δx趨近零，滿足穩定性。
absolute stability：絕對穩定性。特定大小的Δt與Δx，滿足穩定性。
stiff stability：剛硬穩定性。滿足絕對穩定性，而且，
                 當Δt與Δx趨近無限大，
                 則收斂速度趨近無限大（指數趨近負無限大）。

stiffness

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

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

interval of absolute stability

當遞迴公式恰是線性遞迴函數，特徵方程式的根皆介於±1，才會滿足指數穩定性。以此推導Δt與Δx的適當範圍。

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

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

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

(f(t + Δt) - f(t)) / Δt = f(t + Δt)   implicit 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 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 coefficient
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)

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

region of absolute stability

有些教科書喜歡畫一種圖，叫做region of absolute stability，簡稱stable region。此圖毫無實際用處，只是想讓讀者轉換心情。

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

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

套用implicit Euler method，得到左圖。

套用explicit Euler method，得到右圖。

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

【圖片尚待修正】

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

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

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

stability analysis with Fourier transform
（von Neumann stability analysis）

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

d             d
—— f(t,x) + k —— f(t,x) = 0     k is a coefficient
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])

無論如何不可能穩定。

g(w) = 1 - ½ C (exp(𝑖wΔx) - exp(-𝑖wΔx))     amplification factor
g(w) = 1 - C 𝑖 sin(wΔx)                     amplification factor
a numerical scheme is stable iff |g(w)| ≤ 1
this numerical scheme is unstable since |g(w)| > 1

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

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

g(w) = ½ (exp(𝑖wΔx) + exp(-𝑖wΔx)) - ½ C (exp(𝑖wΔx) - exp(-𝑖wΔx))
g(w) = cos(wΔx) - C 𝑖 sin(wΔx)

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

A-stability

專著《Numerical Methods for Ordinary Differential Equations》。

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

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

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

甚至深入觀察經典的演算法的穩定性。

AN-stability：滿足B穩定性的情況下，
            　multistage method的數學性質。
G-stability：滿足A穩定性的情況下，
           　multistep method的數學性質。
algebraic stability：滿足A穩定性的情況下，
                   　multistage method的數學性質。

https://encyclopediaofmath.org/wiki/Non-linear_stability_of_numerical_methods
https://www.jstor.org/stable/2156562
https://www.jstor.org/stable/2156895
http://cds.cern.ch/record/1069162

function of absolute stability

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

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

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

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

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

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

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

穩定性有三種層次。總之都是收縮。

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

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

(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

第三種穩定性分析衍生許多定理。有興趣的讀者請自行研究。

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

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

專著《Numerical Mathematics》。

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

數值分析基本定理，數值解版本，針對絕對穩定性：

(1) Courant–Friedrichs–Lewy condition：
    初始值問題IVP、數學之守恆律、平滑解
    某些線性方案、CFL condition

(2) strong stability preserving property：
    初始值問題IVP、數學之守恆律、平滑解
    Euler method家族、CFL condition

numerical scheme

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

具體來說有四個項目。

(1) property of differential equation
(2) property of dynamical system
(3) numerical scheme
(4) numerical technique

一、微分方程式的數學性質。
二、動態系統的數學性質。
三、遞迴公式設計技巧。
四、遞迴公式求解技巧。

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) smooth solution：
    平滑解。無窮次可微函數。
(b) weak solution：
    弱解。也可以是不可微函數。

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

順帶一提，數值解無法定義連續、可微。

【名稱尚待查詢】

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

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

differentiable solution

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

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

smooth solution

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

理論上來說，符號解可以是平滑函數（處處無窮次可微），也可以是分段函數（銜接處N次可微）。然而目前無人討論分段函數。

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

數學家從數值分析基本定理出發，利用微分運算的泰勒近似，求得微分方程式的數值解。數值解是近似解。全域誤差是近似誤差。想讓全域誤差趨近零，一種方式是採用平滑解。前面章節已經提及。

有限差分法可以處理平滑解。

weak solution

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

一個函數L(t,x)乘以一個平滑函數ϕ(t,x)、積分，得到新函數。一個函數經過包裝，強制變成可微函數。

⌠+∞ ⌠+∞ 
⎮   ⎮   L(t,x) ϕ(t,x) dt dx
⌡-∞ ⌡-∞ 

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

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    ⎠

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

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

⎧ 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
  ⌡-∞

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

有限體積法可以處理平滑解和弱解。

smooth approximation

為了使用有限元素法處理弱解，這裡介紹一種作弊的方法。

不可微函數強行變成平滑函數（弱解強行變成平滑解）。

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

conservative formulation【尚待確認】

為了使用有限體積法處理多值解，這裡介紹一種作弊的方法。

非函數強行變成函數（多值解強行變成單值解）。

不看單一地點的函數值，只看區間總和。

property of dynamical system

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

(1) diminishing
    遞減。符號解範數每回合變小或不變。
(2) local extremum diminishing
    局部極值遞減。符號解每回合極值數量不增加、極值大小不變廣。
(3) positivity preserving
    保正。符號解每回合保持非負數（保持正負號）。
(4) monotonicity preserving
    保單調。符號解每回合保持單調。

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

這些詞彙可以是動名詞、也可以是形容詞。形容詞可以放在be動詞之後、也可以放在名詞之前。放在名詞之前的情況下，大家習慣接上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 coefficient and k < 0
dt

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

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

diminishing if (1 + Δt k) ≤ 1

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

遞減的衍生性質。

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

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

for homogeneous scheme,
contractive => diminishing

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

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

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

local extremum diminishing（range diminishing）

專著《Computational Aerodynamics》。

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

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

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

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

1D
⎧ fᵢ⁽ⁿ⁺¹⁾ = cᵢ₋₁ fᵢ₋₁⁽ⁿ⁾ + cᵢ fᵢ⁽ⁿ⁾ + cᵢ₊₁ fᵢ₊₁⁽ⁿ⁾
⎨ cᵢ₋₁ + cᵢ + cᵢ₊₁ = 1
⎩ cᵢ₋₁ , cᵢ , cᵢ₊₁ ≥ 0

2D
⎧ fᵢ,ⱼ⁽ⁿ⁺¹⁾ = cᵢ,ⱼ fᵢ,ⱼ⁽ⁿ⁾
⎪           + cᵢ₋₁,ⱼ fᵢ₋₁,ⱼ⁽ⁿ⁾ + cᵢ₊₁,ⱼ fᵢ₊₁,ⱼ⁽ⁿ⁾
⎨           + cᵢ,ⱼ₋₁ fᵢ,ⱼ₋₁⁽ⁿ⁾ + cᵢ,ⱼ₊₁ fᵢ,ⱼ₊₁⁽ⁿ⁾
⎪ cᵢ,ⱼ + cᵢ₋₁,ⱼ + cᵢ₊₁,ⱼ + cᵢ,ⱼ₋₁ + cᵢ,ⱼ₊₁ = 1
⎩ cᵢ,ⱼ , cᵢ₋₁,ⱼ , cᵢ₊₁,ⱼ , cᵢ,ⱼ₋₁ , cᵢ,ⱼ₊₁ ≥ 0

知名範例是Laplacian運算。

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 - C) fᵢ⁽ⁿ⁾ + ½ C fᵢ₋₁⁽ⁿ⁾ + ½ C fᵢ₊₁⁽ⁿ⁾
where C = Δt/Δx

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

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

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

   ⎧ fᵢ⁽ⁿ⁺¹⁾ = sumⱼ cᵢⱼ fⱼ⁽ⁿ⁾
   ⎨ sumⱼ cᵢⱼ = 1
   ⎩ cᵢⱼ ≥ 0

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

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

=> order preserving (monotone)
   遞迴公式是單調遞增函數（任一變數皆如此）
   兩個數值解每回合保持高低關係

   ⎰ fᵢ⁽ⁿ⁺¹⁾ = Fᵢ(fⱼ₁⁽ⁿ⁾, ..., fⱼₘ⁽ⁿ⁾)
   ⎱ fⱼₖ⁽ⁿ⁾ ≤ 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)

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

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

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

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

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

   order preserving (monotone) & 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ᵢ|

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

monotonicity preserving

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

另一種解釋。數值解拆解成單調區間們與局部極值們。單調區間保持單調區間，局部極值保持局部極值。

保單調的衍生性質。針對線性遞迴公式。

    monotonicity preserving
<=> positive combination (positivity)

    monotonicity preserving & neighborhood (3-point) & diminishing
<=> local extremum diminishing (LED)

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

Harten's TVD theorem：針對線性遞迴公式，數值解邊界維持不變，數值解內部滿足total variation diminishing，那麼保單調。

進一步推理，同樣的條件之下，保單調等價於total variation diminishing。

   Dirichlet boundary condition & total variation diminishing (TVD)
   數值解兩端固定＆數值解總差值縮小或不變（事實上只能維持不變）

   ⎰ f̂ᵢ⁽ⁿ⁾ = f̂ᵢ⁽ⁿ⁺¹⁾  for all i at boundary
   ⎱ ‖f̂⁽ⁿ⁾‖ᴛᴠ ≤ ‖f̂⁽ⁿ⁺¹⁾‖ᴛᴠ

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

   (f̂⁽ⁿ⁾ is monotone => f̂⁽ⁿ⁺¹⁾ is monotone)

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

[Harten 1983]

    monotonicity preserving
<=> total variation diminishing (TVD)
under Dirichlet boundary condition 

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

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

[Harten 1983]

第一個衍生沒有設置邊界條件（數值解左右延伸至無限遠）。第二個衍生卻設置邊界條件（數值解兩端固定）。無法一路衍生。

numerical scheme

先前介紹了三個項目：

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

以下介紹更多細節，但是不細講。

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。取新值。需要移項整理，甚至需要一次方程組求解。
     d/dt f(t) = f(t) + 1
---> (f[n+1] - f[n]) / Δt = f[n+1] + 1
                            ^^^^^^

explicit。取舊值。
     d/dt f(t) = f(t) + 1
---> (f[n+1] - f[n]) / Δt = f[n] + 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。式子太多，安排順序。
     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

spatial discretization

compact finite difference：用更少的點數達到更高的精度。
degraded finite difference：達到更低的精度。自虐。

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

temporal discretization

(1) multistep method：取多處斜率。都是左方斜率。
(2) multistage method：求多次斜率。都是右方斜率。
(3) multiderivative method：取多階微分。

Adams method: Maclaurin series
https://mathworld.wolfram.com/AdamsMethod.html

Richardson extrapolation
https://en.wikipedia.org/wiki/Richardson_extrapolation

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

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

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

numerical technique

遞迴公式求解技巧。三種類型。

一、遞迴公式。

(1) finite difference method：有限差分法。求平滑解。使用原本的微分方程式。
(2) finite volume method：有限體積法。求弱解。改成通量。
(3) level set method：等高線法。求多值解。改成勢。

二、取樣擇鄰。多維空間形成網格。

(4) finite element method：有限元素法。取樣擇鄰：隨機點／網格／基底內插。
   （本站稱為函數資料結構mesh discretization。）
(5) meshfree method：無網格法。隨機點／半徑範圍／基底內插。
   （本站稱為函數資料結構meshfree discretization。）
(6) adaptive mesh refinement：自適應網格細緻化。各處根據誤差調整網格精細程度。

三、求解。多維空間形成矩陣求解。

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

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

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。時間變化不必乘上倍率。

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

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

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

【待補範例】

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

shock wave

「衝擊波」或「震波」是一種物理現象。例如空氣衝擊波，局部空氣高速移動、向前推擠、密度陡升。下圖是戰鬥機引起的空氣衝擊波，密度陡升之處形成黑線，黑線隨著戰鬥機前進。

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

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

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

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

symbolic solution of shock wave

專著《Computational Gasdynamics》。

一維空間純量守恆型的數學式子。

1D linear advection equation
du     du
—— + a —— = 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₀

scalar conservation law in one-dimensional space
du   df(u)
—— + ————— = 0
dt    dx

chain rule
du   df du            du         du
—— + —— —— = 0   =>   —— + f′(u) —— = 0
dt   du dx            dt         dx

initial condition u(0,x) = u₀(x)
solution          u(t,x) = u₀(x - f′(u(x,t)) t)
waveform          u(t,x) at certain t
wave speed        f′(u(t,x))
 rightward        f′(u(t,x)) > 0
 leftward         f′(u(t,x)) < 0
 static           f′(u) = 0
wavefront         x - f′(u)t = constant
characteristics   dx/dt = f′(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) function curve           u(t,x) for each t
(2) characteristic curves    dx/dt = f′(u) for each x

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

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

scalar conservation law in one-dimensional space
d/dt u = - d/dx f(u)   where f is known

solution
u(t,x)

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

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

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

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

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

Riemann problem:
(1) Rankine–Hugoniot jump condition
(2) Lax admissibility condition
(3) Riemann initial 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
(是咧畫啥物潲)

Riemann problem of system of conservation laws in 1D space
has new concept called discontinuity.
i.e. tangential discontinuity, contact 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》

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/

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

波形是monotonicity preserving、L_∞-diminishing、TVD。

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

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

waveforms are L∞-diminishing and TVD [Glimm-Lax 1971]
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

numeric solution of shock wave

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

conservation law：「守恆律」。微分方程式，物理量守恆。例如質量守恆、動量守恆、能量守恆。

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

conservation：「守恆」。解的數學性質。一維純量守恆律，守恆的意義是數值解每回合數值總和是定值。

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

weak solution：「弱解」。解可以是不連續函數。主要用來討論衝擊波。

admissible solution (entropy solution)：近期稱作「容許解」。早期稱作「熵解」。一種特殊的弱解，而且是唯一解。

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

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

finite difference method

有限差分法，可以處理平滑解。

兩個演算法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用來確保LED，進而確保monotonicity preserving。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條件的數學式子：波動方程式是v Δt / Δx ≤ 1。熱傳導方程式是2 v ∆t / (∆x)² ≤ 1。其中v是速度。

證明請見講義：

【待補範例】

CFL條件的幾何意義：移動距離v Δt不得超過取樣距離Δx，否則空間離散化會抓錯函數值。

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

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

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

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

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

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

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

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

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

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

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

explicit multistage method
⎧ f“⁰” = f⁽ⁿ⁾                                  skip subscript
⎨ f“ⁱ” =  sum  { aᵢⱼ f“ʲ” + Δt bᵢⱼ F(f“ʲ”) }   1 step
⎪       j=1⋯i-1                                already expanded
⎩ f⁽ⁿ⁺¹⁾ = f“ᵐ”                                m stages

convex combination property of each stages
⎧ aᵢⱼ ≥ 0       for all i and j
⎨ bᵢⱼ ≥ 0       for all i and j
⎪ sum aᵢⱼ = 1   for all i
⎩  ʲ

strong stability preserving property:

  if 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:

  ‖f“ⁱ”‖
= ‖sum { aᵢⱼ f“ʲ” + Δt bᵢⱼ F(f“ʲ”) }‖           ith stage
≤ ‖aᵢⱼ‖ ‖sum { f“ʲ” + Δt (bᵢⱼ/aᵢⱼ) F(f“ʲ”) }‖   三角不等式
=  aᵢⱼ  ‖sum { f“ʲ” + Δt (bᵢⱼ/aᵢⱼ) F(f“ʲ”) }‖   aᵢⱼ ≥ 0
≤  aᵢⱼ   sum ‖ f“ʲ” + Δt (bᵢⱼ/aᵢⱼ) F(f“ʲ”) ‖    正負抵銷
≤  aᵢⱼ   sum ‖f“ʲ”‖                             原因如後
≤        max ‖f“ʲ”‖                             凸組合

⎧ ‖f“¹”‖ ≤ max(‖f“⁰”‖)                        = ‖f“⁰”‖
⎨ ‖f“²”‖ ≤ max(‖f“⁰”‖, ‖f“¹”‖)                ≤ ‖f“⁰”‖
⎪    :           :
⎩ ‖f“ᵐ”‖ ≤ max(‖f“⁰”‖, ‖f“¹”‖, ..., ‖f“ᵐ⁻¹”‖) ≤ ‖f“⁰”‖

‖f⁽ⁿ⁺¹⁾‖ = ‖f“ᵐ”‖ ≤ ‖f“⁰”‖ = ‖f⁽ⁿ⁾‖

for each j,
explicit Euler method reaches contractive stability
‖f“ʲ” + Δt (bᵢⱼ/aᵢⱼ) F(f“ʲ”)‖ ≤ ‖f“ʲ”‖
under CFL condition Δt (bᵢⱼ/aᵢⱼ) ≤ ΔtCFL

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

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

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

finite volume method

有限差分法，無法處理弱解。函數跳躍處的差分泰勒展開，截斷誤差不受控制，收斂性不受控制。

有限體積法，可以處理弱解。物理量守恆，而且流量不受跳躍差距影響。

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

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

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

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

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

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ᵢ⁽ⁿ⁾沒有做定積分，形成數量平均數。

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ᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)

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

conservation：守恆型滿足物理量守恆，數值解每回合總和保持相同。通量的相鄰差的區間和等於頭尾相減，其餘兩兩抵銷。頭尾是零導致守恆。

assume f₋₁⸝₂⁽ⁿ⁾ = fₓ₊₁⸝₂⁽ⁿ⁾ = 0

  sum uᵢ⁽ⁿ⁺¹⁾
= sum { uᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾) }
= sum uᵢ⁽ⁿ⁾ - C sum { fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾ }
= sum uᵢ⁽ⁿ⁾ - C (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ᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)

conservation form: flux function
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (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

numerical method：有限體積法的精髓在於自訂通量。嘗試各式各樣通量，得到各式各樣演算法。本文介紹四種演算法。

(1) finite difference method & flux limiter
(2) Godunov's method
(3) Godunov's method & slope limiter
(4) essentially non-oscillatory method (ENO)

finite difference method：先前小節的演算法可以改寫成守恆型。順帶一提，upwind scheme不守恆，無法改寫成守恆型。

Lax–Friedrichs method
可以收斂至容許解。
fᵢ₊₁⸝₂⁽ⁿ⁾ = (fᵢ⁽ⁿ⁾ + fᵢ₊₁⁽ⁿ⁾) / 2
          - (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) / C / 2

Lax–Wendroff method
無法收斂至容許解。
fᵢ₊₁⸝₂⁽ⁿ⁾ = (fᵢ⁽ⁿ⁾ + fᵢ₊₁⁽ⁿ⁾) / 2
          - (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) / C / 2
          + f′ᵢ₊₁⸝₂⁽ⁿ⁾ (fᵢ₊₁⁽ⁿ⁾ - fᵢ⁽ⁿ⁾) * C / 2

flux limiter：有限差分法，引入flux averaging。微調權重以微調通量。特殊權重可以確保monotonicity preserving。不幸的是，大家沒有仔細檢查這些方法是否收斂至容許解。

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

Lax–Wendroff
Warming–Beam
Harten's incremental form
Sweby

Kurganov–Tadmor

https://en.wikipedia.org/wiki/Flux-corrected_transport
https://en.wikipedia.org/wiki/MUSCL_scheme

Godunov's method：兩兩相鄰格子皆視作Riemann problem。

此方法可以視作upwind scheme的推廣版本（此處的上風不是差分地點、而是通量地點）。當守恆律是linear advection equation，此方法退化為upwind scheme，證明在此。

講義《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ᵢ₊₁⁽ⁿ⁾

Roe's method
fᵢ₊₁⸝₂⁽ⁿ⁾ = Roe_flux(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾)

          = ⎰ min(f(uᵢ⁽ⁿ⁾), f(uᵢ₊₁⁽ⁿ⁾)) , if uᵢ⁽ⁿ⁾ < uᵢ₊₁⁽ⁿ⁾
            ⎱ max(f(uᵢ⁽ⁿ⁾), f(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, ω))

https://ocw.mit.edu/courses/16-920j-numerical-methods-for-partial-differential-equations-sma-5212-spring-2003/bafff8bb84c0ee8ac362cea89c5fea38_lec12_notes.pdf

Godunov's method：推廣版本。仔細調整水面。三個步驟。

講義《Computational Magnetohydrodynamics》。

Godunov's method
(1) cell reconstruction：每個格子，求出水面形狀。
    use uᵢ⁽ⁿ⁾ get curve from uᵢ₋₁⸝₂⁺ to uᵢ₊₁⸝₂⁻
(2) cell interface flux evaluation：每個相鄰格子介面，求出通量。
    fᵢ₊₁⸝₂⁽ⁿ⁾
(3) update：套用遞迴公式，求出新的水面平均值。
    uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)

Initialize:
 - Define the number of cells (N) and the cell width (Δx).
 - Set the initial conditions in the cells.

Time-stepping loop:
  while (time < final_time)
    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 state in cell i,
        and Fᵢ₊₁⸝₂ and Fᵢ₋₁⸝₂ are the fluxes
        at the interfaces of cell i.
    - Update the time: time = time + Δt
    - Apply boundary conditions, if necessary.

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 flux基於monotonicity preserving假設。
    根據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。
                基於monotonicity preserving。
Harten–Lax–van Leer–Einfeldt solver (HLLE)：近似解。
Harten–Lax–van Leer contact solver (HLLC)：近似解。

大家也針對特定的微分方程式，發明特定的Riemann solver。例如Euler equations of gas dynamics的Riemann solver：

Roe solver：正解。需要計算微分矩陣的反矩陣，時間複雜度較高。
Osher solver
advection upstream splitting method (AUSM)

https://en.wikipedia.org/wiki/Advection_upstream_splitting_method
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

consistent conservative sheme

本節討論得到弱解的條件。

consistent conservative scheme：一致守恆方案。能夠得到弱解的的遞迴公式。

根據下面兩個定理，一致守恆方案可以得到弱解。【尚待確認】

一、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 (any 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。若數值解滿足L_∞-diminishing與TVD，則數值解滿足L₁收斂性。一維空間的情況下，光是TVD足以滿足L₁收斂性。

做個總結。一致守恆方案收斂至弱解，當數值解滿足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

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̂⁽ⁿ⁾‖ᴛᴠ ≤ C

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

L₁-convergence：數值解與符號解每回合L₁距離趨近零。

‖f⁽ⁿ⁾ - f̂⁽ⁿ⁾‖₁ → 0 when n≥0

where ‖f‖₁ = sumᵢ |fᵢ|

注意到，此定理僅適用於一維純量守恆律。此定理無法直接用於多維空間、向量、方程組，無法直接引入散度、旋度、梯散，無法直接指定取樣間距、網格造型、邊界條件。必須另行證明或反駁。

守恆律的數值分析基本定理是數學界的大難題。每種守恆律都必須特地檢查存在性、唯一性、一致性、穩定性、收斂性。

admissible scheme（entropy satisfying scheme）

本節討論得到容許解的條件。

admissible scheme (entropy satisfying scheme)：容許方案（熵滿足方案）。能夠得到容許解（熵解）的遞迴公式。

entropy inequality：符號解滿足熵不等式則是容許解。先前提到的Oleinik admissibility condition其實是熵不等式其中一種特例。

entropy pair (η,q)
η(u) is entropy.       η″(u) > 0 (strictly convex)
q(u) is entropy flux.  q′(u) = f′(u) η′(u)

entropy inequality
d/dt η(u) + d/dx q(u) ≤ 0   for all entropy pairs
所有凸函數皆使得不等式成立。

admissible solution (entropy solution)
d/dt η(u)ϕ + d/dx q(u)ϕ ≥ 0 for all test function ϕ ≥ 0
弱解的其中一個特例：容許解。

Kruzhkov entropy inequality <=> Oleinik admissibility condition
⎧ d/dt η(u) + d/dx q(u) ≤ 0
⎨ η(u, k) = |u - k|
⎩ q(u, k) = (f(u) - f(k)) sgn(u - k)
其中一種凸函數。導致L₁-contraction。
（用來證明容許解是唯一解。）

Lₚ entropy
η(u) = |u|ᵖ / p
其中一種凸函數。導致Lₚ-contraction。

admissible solution => Lₚ-diminishing
容許解具備Lₚ-diminishing。
畢竟守恆律沒有常數項。

admissible solution => L∞-diminishing & total variation diminishing (TVD)
容許解具備L∞-diminishing和TVD。
先前已經談過。

[DiPerna 1983]
一、當任何一個凸函數使得不等式成立，而且總變差受限，就保證符號解L∞ norm受限。
二、當符號解L∞ norm受限、總變差受限，就保證符號解L₁ norm收斂。（穩定即收斂）

   single entropy inequality & bounded variation

   ⎰ d/dt η(u) + d/dx q(u) ≤ 0 for any (η,q)
   ⎱ ‖u⁽ⁿ⁾‖ᴛᴠ ≤ B when n≥0

=> bounded L∞ norm & bounded variation
   符號解每回合振幅／總差值絕對值受限

   ⎰ ‖u⁽ⁿ⁾‖∞ ≤ M when n≥0
   ⎱ ‖u⁽ⁿ⁾‖ᴛᴠ ≤ B when n≥0

=> L₁-convergence
   符號解絕對值總和收斂至零

   ‖u⁽ⁿ⁾‖₁ → 0 when n→∞

[DiPerna 1983]

https://ntrs.nasa.gov/api/citations/19840010891/downloads/19840010891.pdf
https://ocw.mit.edu/courses/16-920j-numerical-methods-for-partial-differential-equations-sma-5212-spring-2003/b83f5998bf019571a250acd9acb7aa46_lec11_notes.pdf
https://math.stackexchange.com/questions/4083618/

https://math.stackexchange.com/questions/18395/
https://math.stackexchange.com/questions/271658/
https://math.stackexchange.com/questions/782558/
https://math.stackexchange.com/questions/1529156/
L1-bounded => L2-bounded
L2-convergence => L1-convergence

discrete entropy inequality：數值解滿足離散熵不等式則是容許解。值得一提的是，目前已知有三類遞迴公式，滿足CFL condition的情況下，其數值解恰好滿足離散熵不等式。這類遞迴公式稱作entropy satisfying scheme。

conservation
sumᵢ ûᵢ⁽ⁿ⁾ = sumᵢ ûᵢ⁽ⁿ⁺¹⁾

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

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

consistent conservative scheme
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (fᵢ₊₁⸝₂⁽ⁿ⁾ - fᵢ₋₁⸝₂⁽ⁿ⁾)
where C = Δt / Δx
      f̃(u,u) = f(u)
      fᵢ₊₁⸝₂ = f̃(uᵢ,uᵢ₊₁)

discrete entropy inequality [Tadmor 1984]
η(uᵢ⁽ⁿ⁺¹⁾) ≤ η(uᵢ⁽ⁿ⁾) - C (qᵢ₊₁⸝₂⁽ⁿ⁾ - qᵢ₋₁⸝₂⁽ⁿ⁾)
where entropy pair (η,q)
      q̃(u,u) = q(u)
      qᵢ₊₁⸝₂ = q̃(uᵢ,uᵢ₊₁)

entropy satisfying scheme
目前有三種滿足discrete entropy inequality的方法
(1) artifical viscosity
    追加二次偏微分項，稱作人工黏性。
    缺點是需要手工調整係數。
(2) monotone scheme
    僅滿足其中一種凸函數，達成L₁收斂性。
    針對線性遞迴公式，至多一階精準。
(3) Osher's E scheme
    滿足所有凸函數，達成所有收斂性。
    針對線性遞迴公式，至多一階精準。

[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

artifical viscosity：此處不介紹。這裡放一篇Lax的文章。

monotone scheme：得到L₁收斂性。

   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))
二、數值解相減，結果有正有負，拆開正負，分開累計。

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

[Harten 1983]
數值解相減，令被減數是減數往右移位，硬是湊成相鄰差。

monotone自然就會有TVD了，其實不需要特別說明。
又是你，Harten......

   consistent & conservative & TVD
=> L1-difference converge to weak solution

[LeVeque 1992] theorem 15.2
基本上胡說八道

https://www.researchgate.net/publication/285599438
https://www.researchgate.net/publication/230676056

Osher's E scheme：通量不逾越Godunov flux。滿足discrete entropy inequalty所有凸函數。收斂至容許解。

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【尚待確認】

E scheme => L∞-diminishing and TVD

for consistent conservative scheme,
3-point Osher's E scheme <=> 3-point monotone scheme

(=>) [Abdi–Hansen–Schroll 2020]
An adaptive viscosity E-scheme for balance laws
http://www.math.ualberta.ca/ijnam/Volume-17-2020/No-3-20/2020-03-08.pdf

(<=)
https://math.tifrbng.res.in/~praveen/slides/gian2017_estable.pdf
since g(u,v) is non-decreasing in u and non-increasing in v
⎰ g(u,v) ≤ g(u,w) ≤ g(w,w) = f(w)   if u ≤ w ≤ v
⎱ g(u,v) ≤ g(w,v) ≤ g(w,w) = f(w)   if u ≤ w ≤ v
sign(v - u)(g(u,v) - f(w)) ≤ 0   for all w ∈ [u,v]

(<=)
https://people.math.ethz.ch/~hiptmair/tmp/NUMHYP_07.pdf
since F(v,w) is non-decreasing in v and non-increasing in w
⎰ F(v,w) - F(u,u) ≤ 0   if v < u < w
⎱ F(v,w) - F(u,u) ≥ 0   if w < u < v

finite element method

我還沒有學會。這是目前最接近此主題的台灣學校課程：

finite volume method：先前小節的演算法可以改寫成基底內插函數。

discontinuous Galerkin method (DG method)
spectral difference method / spectral volume 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

level set method

我還沒有學會。關鍵字也許是capturing multivalued solution。

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

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

一些古怪的名詞

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 method
大於一階精度的演算法們。例如MUSCL、DG method。
Harten

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

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

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ₓₓ
                 e.g. advection equation   fₜ = fₓ

微分方程式的分類：函數運算

algebraic equation：代數方程式。不做微分運算。
differential equation：微分方程式。只做微分運算。
differential-algebraic equation：微分代數方程式。兩個一起上。

differential equation：微分方程式。只做微分運算。
integral equation：積分方程式。只做積分運算。
integro-differential equation：積分微分方程式。兩個一起上。

differential equation：微分方程式。只做微分運算。
difference equation：差分方程式。只做差分運算。
differential-difference 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₀)沿著向量場跑。

微分方程式的調整

介紹幾個手法。

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

經典的微分方程式：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

Laplace equation

「拉普拉斯方程式」。處處梯散為零。

一坨東西的勢力均衡。

Poisson equation

「泊松方程式」。處處梯散已知。

常見用法是兩邊梯散相同∆f₁ = ∆f₂，f₂已知或∆f₂已知。

兩坨東西的勢力分布相等。其中一坨已知。

梯散反運算。已知梯散g，求得原函數f。

符號解：格林函數疊加而得。

數值解：擁有特殊演算法。時域一次方程組、頻域傅立葉轉換。

Helmholtz equation

「亥姆霍玆方程式」。梯散運算的特徵函數。

梯散運算的特徵函數f是各種頻率的cos和sin波。

符號解：f(x) = e^{𝑖√λ‖x‖} / 4π‖x‖。

經典的微分方程式：diffusion

引入了空間變數、時間變數，形成了守恆定律。

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

Helmholtz equation  f = ∆f     f   = fₓₓ + f₝₝ + f₞₞
Heat equation       ḟ = ∆f     fₜ  = fₓₓ + f₝₝ + f₞₞
Wave equation       f̈ = ∆f     fₜₜ = fₓₓ + f₝₝ + f₞₞

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
   勢變化

經典的微分方程式：relativity

物理元件的互動。

Maxwell's equations               電磁關係
Einstein field equation           重力時間關係
Schrödinger equation              波粒關係

Maxwell's equations

「電磁場方程組」。科學家們先後發現各式各樣的電磁關係式，由Maxwell集結成一套完整理論，由Heaviside精煉成散度與旋度。

因為是Maxwell，所以要加apostrophe s？我不知道原因。

沒有電磁感應的時候，電場靜止。靜電場散度為零。靜電場旋度為零，故靜電場是梯度場／守恆場，故靜電場可以改寫成電勢。

發生電磁感應的時候，電場運動。動電場散度依然為零。動電場旋度不再為零，故動電場不再是梯度場／守恆場。電場空間旋度恰等於磁場時間導數。電場空間變化導致磁場時間變化。電磁對調亦然。

電磁場方程組用來描述動電場。

⎧ ∇∙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

heat equation額外乘上虛數𝑖。

heat equation           ḟ = ∆f      fₜ  = fₓₓ + f₝₝ + f₞₞
Schrödinger equation    ḟ = 𝑖 ∆f     fₜ  = 𝑖(fₓₓ + f₝₝ + f₞₞)

經典的微分方程式：wave propagation

level set

eikonal equation

eikonal equation

經典的微分方程式：shock wave

符號解不連續。

Burgers' equation                 衝擊波
traffic-flow equation             衝擊波
Friedlander's 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

交通車流問題。車輛不會互相穿越，整個過程都是函數。

經典的微分方程式：solitary wave

符號解不動點。

Korteweg–de Vries equation        孤立波
shallow water equations           淺水波
Boussinesq's equations            淺水波
de Saint-Venant equations         淺水波

Korteweg–de Vries equation

符號解是solitary wave與cnoidal wave。

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

經典的微分方程式：glow

引入了合體技的概念。

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)

reaction–diffusion–advection equation

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

Turing pattern
Gray–Scott model
Ginzburg–Landau equation
Kuramoto–Sivashinsky equation  熱流交融

Turing pattern

https://www.nature.com/articles/ncomms1289

Gray–Scott model

https://pmneila.github.io/jsexp/grayscott/

Ginzburg–Landau equation

Kuramoto–Sivashinsky equation

經典的微分方程式：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

注意到此式無法用來製作動畫。

經典的微分方程式：attractor

亂繞圓圈的路線。

Lorenz equation                   大氣對流
Lorenz–Emanuel system             大氣對流
Rabinovich–Fabrikant equation     雙角