f(t) 5 ≤ t ≤ 10 time is an 1D interval
space is a 0D point
(spatial boundary is a point)
f(t,x) 5 ≤ t ≤ 10 time is an 1D interval
-10 ≤ x ≤ +10 space is an 1D interval
(spatial boundary is two endpoints)
f(t,x,y) 5 ≤ t ≤ 10 time is an 1D interval
x² + y² ≤ 10 space is a 2D disk
(spatial boundary is a circle)
隱式法。
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)
differential equation - dynamical system🚧
numerical error
數值誤差。數值遞推過程，符號解與數值解的差距。
symbolic solution / numeric solution
符號解、數值解，畫成圖片。
符號解、數值解，其數學符號的正式寫法。
symbolic solution:
f(t, x) function
f(t₀ + n Δt, x₀ + i Δx) function value at a point
numeric solution:
f̃(t, x) function
f̃(t₀ + n Δt, x₀ + i Δx) function value at a point
數學符號可以精簡成各種風格。
symbolic solution:
f(t₀ + n Δt, x₀ + i Δx) --> f(tₙ, xᵢ) --> f⁽ⁿ⁾(xᵢ)
numeric solution:
f̃(t₀ + n Δt, x₀ + i Δx) --> f̃(tₙ, xᵢ) --> f̃⁽ⁿ⁾(xᵢ) ─╮
╭──────────────────╯
╰-> f[n][i] --> fᵢ⁽ⁿ⁾
local error / global error
數值誤差，分為兩個階段：
local error: distance of first iteration of
a symbolic solution and a numeric solution.
‖f⁽¹⁾ - f̃⁽¹⁾‖
global error: distance of nth iteration of
a symbolic solution and a numeric solution.
‖f⁽ⁿ⁾ - f̃⁽ⁿ⁾‖
consistency: 1. theorems have no contradiction.
2. an equation have a unique solution.
3. symbol and numeral are approximately equal.
stability: two sequence are bounded by a certain distance
that is depending on initial distance
|aₙ - bₙ| < ε when |a₀ - b₀| < δ
where ε and δ are constants
convergence: a sequence leads to a constant.
aₙ → c when n → ∞
where c is a constant
differential equation
L(f(t,x)) = 0
symbolic solution
f(t,x)
discretization of differential equation
L̃(f(t,x)) = 0
numeric solution
f̃(t,x)
consistency: difference of the functional equation and
its discretization approaches zero.
L(f(t,x)) - L̃(f(t,x)) → 0 when Δt→0, Δx→0
where L(f(t,x)) = (d/dt f(t,x)) - F(f(t,x))
big-O notation is used somehow.
L(f(t,x)) = L̃(f(t,x)) + O((Δt)ᵖ + (Δx)𐞥)
where p>0, q>0
stability: distance of corresponding iteration of
two numeric solutions with
different initial conditions is bounded.
‖f̃⁽ⁿ⁾ - g̃⁽ⁿ⁾‖ ≤ k ‖f̃⁽⁰⁾ - g̃⁽⁰⁾‖ when n≥0
where k is a constant
with consistency,
distance of corresponding iteration of
a symbolic solution and a numeric solution with
same initial condition is bounded.
‖f⁽ⁿ⁾ - f̃⁽ⁿ⁾‖ ≤ k ‖f⁽⁰⁾ - f̃⁽⁰⁾‖ when n≥0
where k is a constant
convergence: distance of corresponding iteration of
two numeric solutions with
different initial conditions leads to zero.
‖f̃⁽ⁿ⁾ - g̃⁽ⁿ⁾‖ → 0 when n→∞
with consistency,
distance of corresponding iteration of
a symbolic solution and a numeric solution with
same initial condition leads to zero.
‖f⁽ⁿ⁾ - f̃⁽ⁿ⁾‖ → 0 when n≥0
stability analysis with Fourier transform （von Neumann stability analysis）
當遞迴公式恰是線性遞迴函數，可以逐項轉換到頻域。
舉例來說，下述微分方程式：
d d
—— f(t,x) + k —— f(t,x) = 0 k is a parameter
dt dx
遞迴公式：
f[n+1][i] = f[n][i] - ½ C (f[n][i+1] - f[n][i-1])
where C = k Δt / Δx is CFL number without absolute value
forward difference of t
central difference of x
傅立葉轉換：
DFT(f[n+1]) = DFT(f[n]) - ½ C (exp(+𝑖ωi) DFT(f[n])
- exp(-𝑖ωi) DFT(f[n]))
DFT(f[n+1]) = (1 - ½ C (exp(+𝑖ωi) + exp(-𝑖ωi)) DFT(f[n])
此遞迴公式無論如何不可能穩定。
a[i] = 1 - ½ C (exp(+𝑖ωi) - exp(-𝑖ωi)) amplification factor
a[i] = 1 - C 𝑖 sin(ωi) amplification factor
numerical scheme is stable iff |a[i]| ≤ 1
this numerical scheme is unstable since |a[i]| > 1
修改遞迴公式，採用Lax–Friedrichs method，就有機會穩定。
f[n+1][i] = ½ (f[n][i+1] + f[n][i-1])
- ½ C (f[n][i+1] - f[n][i-1])
a[i] = ½ (exp(+𝑖ωi) + exp(-𝑖ωi)) amplification factor
- ½ C (exp(+𝑖ωi) - exp(-𝑖ωi))
a[i] = cos(ωi) - C 𝑖 sin(ωi) amplification factor
let |C| ≤ 1 is CFL condition
thus |a[i]|² = cos²(ωi) + C² sin²(ωi)
≤ cos²(ωi) + sin²(ωi)
= 1
thus |a[i]| ≤ 1
a numerical scheme is stable iff |a[i]| ≤ 1
this numerical scheme is stable with CFL condition |C| ≤ 1
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
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。
f[n+1] = a f[n] where a = 1 / (1 - Δt)
再將數值a換成函數ϕ(Δt)。
f[n+1] = ϕ(Δt) f[n] where ϕ(Δt) = 1 / (1 - Δt)
數學家甚至追加參數。微分方程式／遞迴公式附帶參數，判斷各種參數數值的穩定性。
舉例來說，下述微分方程式：
d
—— f(t) = k f(t) k is a parameter
dt
套用implicit Euler method。
f[n+1] = a f[n] where a = 1 / (1 - k Δt)
再將數值a換成函數ϕ(k Δt)。
f[n+1] = ϕ(k Δt) f[n] where ϕ(k Δt) = 1 / (1 - k Δt)
再將函數輸入k Δt換成z。形成函數ϕ(z)。
f[n+1] = ϕ(z) f[n] where ϕ(z) = 1 / (1 - z)
z = k Δt
k可以是複數。z可以是複數。
數學家描述複數，習慣使用符號z。我也採用z。
收縮穩定性是|ϕ(z)| ≤ 1。可以視作一條等高線之內或之外。
數學家描述等高線，習慣使用符號ϕ。我也採用ϕ。
diagram of absolute stability
有些教科書喜歡把stable region畫出來。此圖毫無實際用處，只是想讓讀者轉換心情。
舉例來說，下述微分方程式：
d
—— f(t) = k f(t) k is a parameter
dt
套用implicit Euler method，得到左圖。
套用explicit Euler method，得到右圖。
穩定區域|ϕ(z)| ≤ 1畫在複數平面。座標軸是z的實部虛部。
implicit Euler method:
k = a + b𝑖
z = k Δt
= (a + b𝑖) Δt
= (a Δt) + (b Δt)𝑖
Re[z] = a Δt is defined as x
Im[z] = b Δt is defined as y
ϕ(z) = 1 / (1 - z)
1/ϕ(z) = 1 - z
= 1 - k Δt
= 1 - (a + b𝑖) Δt
= (1 - a Δt) - (b Δt)𝑖
Re[1/ϕ(z)] = 1 - a Δt
Im[1/ϕ(z)] = -b Δt
|ϕ(z)| ≤ 1 contractive stability
=> |1/ϕ(z)| ≥ 1
=> |1/ϕ(z)|² ≥ 1
=> Re[1/ϕ(z)]² + Im[1/ϕ(z)]² ≥ 1
=> (1 - a Δt)² + (-b Δt)² ≥ 1
=> (1 + Re[z])² + (-Im[z])² ≥ 1
=> (1 + Re[z])² + Im[z]² ≥ 1
=> (1 - x)² + y² ≥ r²
=> outside a circle with center (1,0) and radius 1
explicit Euler method:
k = a + b𝑖
z = k Δt
= (a + b𝑖) Δt
= (a Δt) + (b Δt)𝑖
Re[z] = a Δt is defined as x
Im[z] = b Δt is defined as y
ϕ(z) = 1 + k Δt
= 1 + (a + b𝑖) Δt
= (1 + a Δt) + (b Δt)𝑖
Re[ϕ(z)] = 1 + a Δt
Im[ϕ(z)] = b Δt
|ϕ(z)| ≤ 1 contractive stability
=> |ϕ(z)|² ≤ 1
=> Re[ϕ(z)]² + Im[ϕ(z)]² ≤ 1
=> (1 + a Δt)² + (b Δt)² ≤ 1
=> (1 + Re[z])² + Im[z]² ≤ 1
=> (1 + x)² + y² ≤ r²
=> inside a circle with center (-1,0) and radius 1
此例當中，z = k Δt。k可能是複數，於是畫在複數平面、視作二維向量。Δt只能是實數，於是視作向量縮放倍率。
zero stability: |ϕ(z)| ≤ 1 when z → 0
absolute stability: |ϕ(z)| ≤ 1 for given z
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 → -∞)
針對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)定義穩定性。此例當中：
zero stability: |r(z)| ≤ 1 when z → 0
absolute stability: |r(z)| ≤ 1 for given z
A-stability: |r(z)| ≤ 1 for all Re[z] ≤ 0
L-stability: A-stability and r(z) → 0 as z → -∞
(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 with certain condition.
explicit multistep method is never A-stable.
[Dahlquist 1963]
(3) Dahlquist's order barrier theorem
implicit multistep method with m steps is at most
m+1 (m is odd) / m+2 (m is even) order accurate.
explicit multistep method with m steps is at most
m order accurate.
[Dahlquist 1978]
(4) Dahlquist's order barrier theorem
A-stable implicit multistep method is at most
second order accurate.
[Dahlquist 1963]
differential equation - discretization🚧
discretization
離散化。未知函數實施離散化，符號轉數值，以便計算數值解。
discretization of real number:
(0) number discretization 數字離散化
(a) floating-point arithmetic 四則運算離散化（浮點數運算）
discretization of real function:
(1) domain discretization 定義域離散化
(a) grid (mesh-based method) 方格（網格法）
(b) mesh (mesh-based method) 網格（網格法）
(c) particles (meshfree method) 粒子（無網格法）
(2) codomain discretization 對應域離散化
(a) function value (FDM) 單點函數值（有限差分法）
(b) average function value (FVM) 區間平均函數值（有限體積法）
(c) shape function (FEM) 形狀函數（有限元素法）
(d) sinusoid (spectral method) 弦波（頻譜法）
(3) operator discretization 運算離散化
(a) numerical arithmetic 四則運算離散化
(b) numerical differentiation 微分運算離散化
(c) numerical integration 積分運算離散化
(2) codomain discretization 對應域離散化
(a) function value 單點函數值
at a point
(b) average function value 區間平均函數值
over an interval
(c) shape function 形狀函數
(d) sinusoid 弦波
(1) shape function at samples forms Kronecker delta function.
⎰ Nᵢ(xⱼ) = 1 if i = j
⎱ Nᵢ(xⱼ) = 0 if i ≠ j
(2) sum of all shape functions equals 1.
sum Nᵢ(x) = 1 for all x
ⁱ
Crank–Nicolson method：Euler method with
fᵢ⁽ⁿ⁺¹⸍²⁾ = (fᵢ⁽ⁿ⁺¹⁾ + fᵢ⁽ⁿ⁾) / 2
for all i
Du Fort-Frankel method：leapfrog method with
fᵢ⁽ⁿ⁾ = (fᵢ⁽ⁿ⁺¹⁾ + fᵢ⁽ⁿ⁻¹⁾) / 2
(1) height avoids negative value. (cutoff threshold)
if h < ε, then h = ε.
h = max(ε, h)
(2) velocity avoids division by zero.
if h < ε, then u = 0.
v = ⎰ (hv)/h , if h > ε
⎱ 0 , otherwise
where ε = 10⁻⁶
constant-coefficient advection equation:
d d
—— f(t,x) + k —— f(t,x) = 0 k is a parameter
dt dx
forward-time central-space method (FTCS):
f[n+1][i] = f[n][i] - ½ C (f[n][i+1] - f[n][i-1])
where C = k Δt / Δx is CFL number without absolute value
Lax–Friedrichs method:
1. f[n+1][i] = f̄[n][i] - ½ C (f̄[n][i+1] - f̄[n][i-1])
2. f̄[n+1][i] = f[n+1][i] + ½ (f̄[n][i+1] - 2 f̄[n][i] + f̄[n][i-1])
Lax–Friedrichs method:
f̄[n+1][i] = f̄[n][i] - ½ C (f̄[n][i+1] - f̄[n][i-1])
+ ½ (f̄[n][i+1] - 2 f̄[n][i] + f̄[n][i-1])
Lax–Friedrichs method:
f̄[n+1][i] = ½ (f̄[n][i+1] + f̄[n][i-1])
- ½ C (f̄[n][i+1] - f̄[n][i-1])
shallow water equations:
⎧ ∂/∂t h + ∂/∂x (hu) + ∂/∂y (hv) = 0
⎨ ∂/∂t u + u ∂/∂x u + v ∂/∂y u + g ∂/∂x h = 0
⎩ ∂/∂t v + u ∂/∂x v + v ∂/∂y v + g ∂/∂y h = 0
h: height
u: x velocity
v: y velocity
g: gravity (constant)
differential equation:
df df
a —— + b —— = c
dx dy
chain rule:
dx dy
let a = —— and b = ——
dt dt
dx df dy df df
LHS = —— —— + —— —— = ——
dt dx dt dy dt
characteristics:
df
—— = c
dt
characteristic curve:
f = f₀ + ∫ c dt
differential equation:
d d
a(x,y) —— f(x,y) + b(x,y) —— f(x,y) = c(x,y)
dx dy
differential equation with additional variable t:
d d
a(x(t),y(t)) —— f(x(t),y(t)) + b(x(t),y(t)) —— f(x(t),y(t))
dx dy
= c(x(t),y(t))
chain rule:
d
let a(x(t),y(t)) = —— x(t)
dt
d
b(x(t),y(t)) = —— y(t)
dt
d d d d
LHS = —— x(t) —— f(x(t),y(t)) + —— y(t) —— f(x(t),y(t))
dt dx dt dy
d
= —— f(x(t),y(t))
dt
characteristics:
d
—— f(x(t),y(t)) = c(x(t),y(t))
dt
characteristic curve:
⌠t
f(x(t),y(t)) = f(x(0),y(0)) + ⎮ c(x(t),y(t)) dt
⌡0
https://math.stackexchange.com/questions/86570/
特徵曲線。選定初始值，根據速度，畫出曲線。
特徵法。微分方程式，利用特徵，求得符號解。解微分方程組。
method of characteristics:
⎧ d
⎪ —— x(t) = a(x(t),y(t))
⎪ dt
⎪
⎪ d
⎨ —— y(t) = b(x(t),y(t))
⎪ dt
⎪
⎪ d
⎪ —— f(t) = c(x(t),y(t))
⎩ dt
advection equation:
d d
—— f(t,x) + v(t,x) —— f(t,x) = 0
dt dx
explicit Euler method:
v[n][i] Δt
f[n+1][i] = f[n][i] + —————————— (f[n][i] - f[n][i-1])
Δx
CFL condition for advection equation:
|vₘₐₓ| Δt
————————— ≤ 1 where vₘₐₓ is the extremum value of v
Δx in numerical simulation
exponential decay equation:
d
—— f(t) = k f(t) k is a parameter and k < 0
dt
explicit Euler method:
fᵢ⁽ⁿ⁺¹⁾ - fᵢ⁽ⁿ⁾
——————————————— = k fᵢ⁽ⁿ⁾
Δt
fᵢ⁽ⁿ⁺¹⁾ = (1 + Δt k) fᵢ⁽ⁿ⁾
diminishing if (1 + Δt k) ≤ 1
L∞-diminishing ⇏ total variation diminishing (TVD)
L∞-diminishing ⇍ total variation diminishing (TVD)
四、當遞迴公式不因地點而變，且數值解兩端固定（狄利克雷邊界條件），LED即TVD。
兩者差別在於加上常數函數（整體垂直平移）。利用邊界條件排除此例外。【尚待確認】
under Dirichlet boundary condition,
local extremum diminishing (LED) for every closed intervals
<=> total variation diminishing (TVD) for every closed intervals
constant-coefficient advection equation:
d d
—— u(t,x) + c —— u(t,x) = 0
dt dx
semi-discretization (explicit Euler method):
u(t+Δt,x) - u(t,x) d
—————————————————— + c —— u(t,x) = 0
Δt dx
spectral method:
⎛ u(t+Δt,x) - u(t,x) d ⎞
DFT⎜ —————————————————— + c —— u(t,x) ⎟ = 0
⎝ Δt dx ⎠
改寫成陣列：
(DFT(u[n+1]) - DFT(u[n])) / Δt + c (𝑖ωi/Δx) DFT(u[n]) = 0
DFT(u[n+1]) = (1 + C 𝑖ωi) DFT(u[n])
u[n+1] = IDFT((1 + C 𝑖ωi) DFT(u[n]))
where C = c Δt/Δx is CFL number without absolute value
𝑖 is imaginary unit
ω = 2π/N is fundamental angular frequency of DFT
i = (0, 1, ⋯, N-1) is a sequence from 0 to N-1
1 = (1, 1, ⋯, 1) is a sequence of 1s.
N is sampling number
也可以全部離散化。優點是adaptive mesh refinement容易收斂。【尚待確認】
constant-coefficient advection equation:
d d
—— u(t,x) + c —— u(t,x) = 0
dt dx
full-discretization:
u(t+Δt,x) - u(t,x) u(t,x+Δx) - u(t,x-Δx)
—————————————————— + c ————————————————————— = 0
Δt 2 Δx
spectral method:
⎛ u(t+Δt,x) - u(t,x) u(t,x+Δx) - u(t,x-Δx) ⎞
DFT⎜ —————————————————— + c ————————————————————— ⎟ = 0
⎝ Δt 2 Δx ⎠
改寫成陣列：
(DFT(u[n+1]) - DFT(u[n])) / Δt
+ c (exp(+𝑖ωi) DFT(u[n]) - exp(-𝑖ωi) DFT(u[n])) / (2 Δx) = 0
DFT(u[n+1]) = (1 - ½ C (exp(+𝑖ωi) + exp(-𝑖ωi)) DFT(u[n])
u[n+1] = IDFT((1 - ½ C (exp(+𝑖ωi) + exp(-𝑖ωi)) DFT(u[n]))
[Rossby 1942]
Forecasting of flow patterns in the free atmosphere by a trajectory method.
Basic Principles of Weather Forecasting by V. P. Starr (appendix), pp.268-284
https://photino.cwa.gov.tw/rdcweb/lib/cd/cd07mb/MB/PDF/04/No.1/03.pdf
https://en.wikipedia.org/wiki/Carl-Gustaf_Rossby
Hamilton–Jacobi equation d/dt f + H(x, df/dx, t) = 0
where f has the form f(x(t), t)
eikonal equation |∇f| = n
level set equation d/dt φ + v|∇φ| = 0
https://math.mit.edu/classes/18.086/2007/levelsetpres.pdf
https://math.mit.edu/classes/18.086/2007/levelsetnotes.pdf
https://math.mit.edu/classes/18.376/Notes/LectureTopics18376.pdf
https://en.wikipedia.org/wiki/Viscosity_solution
https://benjaminmoll.com/wp-content/uploads/2019/07/viscosity_slides.pdf
In the case of hyperbolic systems, the notion of weak
solution based on distributions does not guarantee
uniqueness, and it is necessary to supplement it with
entropy conditions or some other selection criterion.
In fully nonlinear PDE such as the Hamilton–Jacobi
equation, there is a very different definition of weak
solution called viscosity solution.
level set ghost point
https://www.researchgate.net/publication/334643942
Stanley Osher
https://link.springer.com/book/10.1007/b98879
James Sethian
https://math.berkeley.edu/~sethian/2006/Publications/Book/book.html
Hailiang Liu
https://faculty.sites.iastate.edu/hliu/
Ken Museth
https://ken.museth.org/Publications.html
Fronts propagating with curvature-dependent speed
https://math.berkeley.edu/~sethian/Papers/sethian.osher.88.pdf
capturing multivalued solution
Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations
https://math.berkeley.edu/~sethian/Papers/sethian.osher.88.pdf
Multi-Valued Solution and Level Set Methods in Computational High Frequency Wave Propagation
https://faculty.sites.iastate.edu/hliu/files/inline-files/Liu-osher-tsai-CICP06_1.pdf
Front Tracking for Hyperbolic Conservation Laws
https://www.uio.no/studier/emner/matnat/math/MAT4380/v15/pensumliste/chapter3.pdf
A Level Set Method for the Computation of Multivalued Solutions to Quasi-Linear Hyperbolic PDEs and Hamilton-Jacobi Equations
https://www.researchgate.net/publication/2855502
需要判斷乾溼。水深歸零稱作乾，反之稱作濕。
(1) wet–dry front tracking：找到乾濕交界之處。
(2) wet–dry transition：乾濕互動方式。
(1) particle is vorticity.
(2) calculate velocity field by Biot—Savart law.
summing the contributions of all particles.
(3) particle advection by velocity field.
updating the positions of particles.
vorticity is constant for non-viscous fluid.
(4) apply boundary condition.
symplectic method
扭對稱法。微分方程式改寫成扭對稱動態系統。
詳情請見後面章節Hamilton's equation。
differential equation - conservation law🚧
本文介紹動態系統的其中一種特例：守恆律
conservation law
物理學之守恆定律：物理量總和固定。
例如質量守恆、動量守恆、能量守恆。不生不滅，相加為常數。
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
linear
conservation law d/dt u(t,x) + a d/dx u(t,x) = 0
initial condition u(0,x) = u₀(x)
solution u(t,x) = u₀(x - at)
waveform u(t,x) at certain t
wave speed a = constant
rightward a > 0
leftward a < 0
static a = 0
wavefront x - at = constant
characteristics dx/dt = a
expansion a u(t,x) ≤ a u(t,y) for x₀ ≤ x ≤ y ≤ y₀
compression a u(t,x) ≥ a u(t,y) for x₀ ≤ x ≤ y ≤ y₀
非線性守恆律、或者說一般的守恆律。
conservation law:
d d
—— u(t,x) + a(u(t,x)) —— u(t,x) = 0
dt dx
d d
—— u(t,x) + f′(u(t,x)) —— u(t,x) = 0 define f′(u) = a(u)
dt dx
d d
—— u(t,x) + —— f(u(t,x)) = 0 chain rule
dt dx
conservation law d/dt u(t,x) + a(u(t,x)) d/dx u(t,x) = 0
conservation law d/dt u(t,x) + d/dx f(u(t,x)) = 0
where f′(u) = a(u)
initial condition u(0,x) = u₀(x)
solution u(t,x) = u₀(x - a(u(t,x))⋅t)
waveform u(t,x) at certain t
wave speed a(u(t,x))
rightward a(u(t,x)) > 0
leftward a(u(t,x)) < 0
static a(u(t,x)) = 0
wavefront x - a(u)t = constant
characteristics dx/dt = a(u)
expansion f(u(t,x)) ≤ f(u(t,y)) for x₀ ≤ x ≤ y ≤ y₀
compression f(u(t,x)) ≥ f(u(t,y)) for x₀ ≤ x ≤ y ≤ y₀
符號解可以畫成函數曲線圖（橫軸位置x、縱軸函數值u）、特徵曲線圖（橫軸位置x、縱軸時間t）。
(1) waveforms (function curves) u(t,x) for each t
(2) characteristic curves dx/dt = a(u) for each x
符號解是LED。LED導致保單調、L∞遞減、TV遞減。
monotonicity preserving
https://books.google.com.tw/books?id=pcwLAQAAQBAJ&pg=PA67
L∞-diminishing and total variation diminishing (TVD)
https://books.google.com.tw/books?id=-gHkBwAAQBAJ&pg=PA291
tame oscillation condition
https://bpb-us-e1.wpmucdn.com/sites.psu.edu/dist/d/80666/files/2019/09/clawtut4-Oxford.pdf
[Glimm–Lax 1970]
Riemann problem = conservation law & Riemann initial condition
conservation law:
d d
—— u(t,x) + —— f(u(t,x)) = 0
dt dx
solution:
u(t,x)
Riemann initial condition:
方便教學的簡易範例。不連續處只有一個，函數值是兩個常數。
u(0,x) = ⎰ uʟ(0) if x ≤ p(0)
⎱ uʀ(0) if x > p(0)
discontinuous position of solution:
p(t)
moving speed of discontinuous position:
s(t) = p′(t)
function value of discontinuous position of solution:
uʟ(t) = lim u(t,x)
x→p(t)⁻
uʀ(t) = lim u(t,x)
x→p(t)⁺
compression:
壓縮。往右推擠。
f(uʟ) ≥ f(uʀ) at certain t
Rankine–Hugoniot jump condition:
左側溢出等於右側溢出。物理量總和守恆。
f(uʟ) - s uʟ = f(uʀ) - s uʀ for all t
Rankine–Hugoniot jump condition:
改寫數學式子，宛如截線斜率，宛如等比例，宛如物理學的平均速度公式。
f(uʀ) - f(uʟ)
————————————— = s for all t
uʀ - uʟ
Oleĭnik entropy condition:
內插新增一處不連續，仍往右推擠。（內插新增多處不連續，則速度s(t)遞減。）
f(u*) - f(uʟ) f(uʀ) - f(u*)
————————————— ≥ ————————————— for all t
u* - uʟ uʀ - u* for all u* that uʟ > u* > uʀ
Lax entropy condition:
a(u)恰是嚴格遞增函數、f(u)恰是嚴格凸函數。得以簡化式子、追加中間項s(t)。
f′(uʟ) > s > f′(uʀ) when f″(u) = a′(u) > 0 for all t
Rankine–Hugoniot jump condition:
f(uʀ) - f(uʟ)
————————————— = s for all t
uʀ - uʟ
Oleĭnik entropy condition:
f(u*) - f(uʟ) f(uʀ) - f(u*)
————————————— ≥ ————————————— for all t
u* - uʟ uʀ - u* for all u* that uʟ > u* > uʀ
新式也等價於下述不知名entropy condition。證明省略。
weak formulation of entropy inequality
<=> an entropy condition
an entropy condition:
q(uʀ) - q(uʟ)
s ≥ ————————————— for all t
η(uʀ) - η(uʟ)
Chapter 8
守恆律的初始值的L²和L∞受限，
其中f二次微分處處不是零（f是凹函數或凸函數）。
（證明是用Compensated Compactness / Young Measure）
Corollary 5.2
一個函數序列uₙ(t,x)，屬於Sobolev space W¹¹：
每個函數，每個時刻，本身、對空間微分、對時間微分，三者的L¹同時受限。
那麼可以找到一個函數子序列：
TV受限、L¹受限、L¹平緩、Lᵖloc逐點收斂至特定函數
Theorem 5.6 BV compactness theorem [Helly]
A sequence {uₙ}⊂BV(Ω) with uniformly bounded TV and L∞
has a subsequence that converges to u in Lᵖloc.
Arzelà–Ascoli compactness theorem:
a sequence {uₙ}
1. uniformly bounded
‖uₙ(x)‖ < M for all n and x
2. uniformly equicontinuous
‖uₙ(x) - uₙ(y)‖ < ε when ‖x - y‖ < δ for all n and x
then exist a subsequence of {uₙ} that uniformly converges to u
‖uₙ(x) - u(x)‖ → 0 for all x when n→∞
dominated convergence theorem:
a sequence {uₙ}⊂L¹(Ω)
1. uniformly bounded by v in L¹
‖uₙ(x)‖ < v(x) for all n and x where v⊂L¹(Ω)
2. pointwise converges to u in Lᵖ
‖uₙ(x) - u(x)‖ → 0 for all x when n→∞
then sequence {uₙ} converges to u in L¹
‖uₙ - u‖₁ → 0 when n→∞
uniformly有時候是指整個數列（全部的n），有時候是指穩定性。
數學元件。
Lebesgue space p次方和仍受限的函數
Lᵖ(Ω) = {f:Ω→ℝ | ∫Ω|f(x)|ᵖdx < ∞}
norm
‖f‖Lᵖ(Ω) = (∫Ω|f(x)|ᵖdx)¹⸍ᵖ
space of locally integrable functions 積分範圍只考慮連通受限區域
Lᵖloc(Ω) = {f∈Lᵖ(K) ∀K⊂⊂Ω}
relatively compact (compact and contained by)
K⊂⊂Ω
Sobolev space 微分1次到微分α次，p次方和仍受限的函數。
Wᵏᵖ(Ω) = {u∈Lᵖ(Ω) | Dαu∈Lᵖ(Ω), ∀|α|≤k}
norm
‖u‖Wᵏᵖ(Ω) = (sum|α|≤k ‖Dαu‖Lᵖ(Ω)ᵖ)¹⸍ᵖ
space of functions of bounded variation 升降高度總和受限
BV(Ω) = {f∈L¹(Ω) | ‖f‖ᴛᴠ < ∞}
total variation
‖f‖ᴛᴠ = sup sum |f(xᵢ) - f(xᵢ₋₁)|
x₀<⋯<xₙ 1⋯n
Banach space complete normed space, e.g. Lᵖ(Ω)
Fréchet space with semi-norm and metric
Hibert space L²(Ω) with inner product
⟨f,g⟩ = ∫Ωf(x)g(x)dx
Hölder's inequality ∫Ω|f(x)g(x)|dx ≤ ‖f‖Lᵖ(Ω)‖f‖L𐞥(Ω)
where 1/p + 1/q = 1
Lᵖ(Ω) ⊂ L𐞥(Ω) if p > q
Lᵖ-stability => L𐞥-stability if p > q
https://math.stackexchange.com/questions/18395/
https://math.stackexchange.com/questions/271658/
https://math.stackexchange.com/questions/782558/
https://math.stackexchange.com/questions/1529156/
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
Finite volume schemes and Lax–Wendroff consistency
https://hal.science/hal-03257774v2/document
躼躼長講了一大堆，上尾仔干焦保證物理量總和守恆。無彩工。
consistent conservative scheme
一致守恆方案，又滿足任意一種收縮穩定性，可以得到其中一個弱解。然並卵。
目前有兩類線性遞迴公式可以得到弱解。
(1) total variation diminishing scheme (TVD scheme)
總變差遞減方案。
數值解每回合總變差變小或不變。
(2) local extremum diminishing scheme (LED scheme)
局部極值遞減方案。
數值解每回合局部極值數值範圍變小或不變。
遞迴公式係數是凸組合。
The set of TVD schemes contains the set of monotone schemes.
However, unlike the case of monotone schemes,
the fact that a scheme is TVD does not automatically imply that
it is consistent with the entropy inequality.
[Harten 1984]
for consistent conservative scheme,
3-point monotone scheme <=> its flux function f̌(u,v)
is non-decreasing in u
and non-increasing in v
[Harten–Hyman–Lax 1976]
3-point consistent conservative monotone scheme:
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - λ (f̌(uᵢ⁽ⁿ⁾, uᵢ₊₁⁽ⁿ⁾) - f̌(uᵢ₋₁⁽ⁿ⁾, uᵢ⁽ⁿ⁾))
where ⎰ d/du f̌(u,v) ≥ 0 (f̌(u,v) is non-decreasing in u)
⎱ d/dv f̌(u,v) ≤ 0 (f̌(u,v) is non-increasing in v)
此時，遞減保序方案＝遞減Godunov方案。證明省略。
for consistent conservative scheme with CFL condition,
3-point monotone scheme <=> 3-point E scheme
(⟹)
https://math.tifrbng.res.in/~praveen/slides/gian2017_estable.pdf
since f̌(u,v) is non-decreasing in u and non-increasing in v
⎰ f̌(u,v) ≤ f̌(u,w) ≤ f̌(w,w) = f(w) if u ≤ w ≤ v
⎱ f̌(u,v) ≤ f̌(w,v) ≤ f̌(w,w) = f(w) if u ≤ w ≤ v
sign(v - u)(f̌(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
(⟸) [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://link.springer.com/chapter/10.1007/978-3-319-96415-7_33
稍後介紹的數值方案，其中一些是熵滿足方案。
Lax–Friedrichs method is monotone scheme when CFL ≤ 1.
Lax-Wendroff scheme is NOT monotone scheme.
[Tadmor 1984]
Godunov's method is E scheme when CFL ≤ 1.
MUSCL is E scheme when CFL ≤ 1.
[Osher 1985]
linear conservation law:
d d
—— u(t,x) + a —— u(t,x) = 0
dt dx
forward-time backward-space method:
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (uᵢ⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾)
forward-time forward-space method:
uᵢ⁽ⁿ⁺¹⁾ = uᵢ⁽ⁿ⁾ - C (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾)
upwind difference method:
uᵢ⁽ⁿ⁺¹⁾ = ⎰ uᵢ⁽ⁿ⁾ - C (uᵢ⁽ⁿ⁾ - uᵢ₋₁⁽ⁿ⁾) , if a > 0
⎱ uᵢ⁽ⁿ⁾ - C (uᵢ₊₁⁽ⁿ⁾ - uᵢ⁽ⁿ⁾) , if a < 0
where C = a Δt / Δx is CFL number without absolute value
conservation law:
d d
—— u(t,x) + —— f(u(t,x)) = 0
dt dx
upwind difference method:
uᵢ⁽ⁿ⁺¹⁾ = ⎰ uᵢ⁽ⁿ⁾ - λ (f(uᵢ⁽ⁿ⁾) - f(uᵢ₋₁⁽ⁿ⁾)) , if f(uᵢ⁽ⁿ⁾) > 0
⎱ uᵢ⁽ⁿ⁾ - λ (f(uᵢ₊₁⁽ⁿ⁾) - f(uᵢ⁽ⁿ⁾)) , if f(uᵢ⁽ⁿ⁾) < 0
where λ = Δt / Δx is CFL number without speed
reconstruction:
(1) piecewise constant function
分段常數函數：單一格子的水面形狀是水平線。形成LED。
形成線性遞迴公式。根據Godunov's order barrier theorem，一階準度。
(2) piecewise linear function with slope limiter
分段一次函數：單一格子的水面形狀是斜線。（不是一次內插）
斜率限制器：所有格子的水面形狀以slope limiter將強行調整成LED。
形成非線性遞迴公式。二階準度。
然而捨棄一致性也無所謂準度了。二階準度只是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)
van Leer's method
of linear conservation law:
fᵢ₊₁⸝₂⁽ⁿ⁾ = ⎰ a (uᵢ⁽ⁿ⁾ + ½ σᵢ⁽ⁿ⁾ (Δx - a Δt)) , if a > 0
⎱ a (uᵢ₊₁⁽ⁿ⁾ - ½ σᵢ₊₁⁽ⁿ⁾ (Δx + a Δt)) , if a < 0
Godunov-type method:
for each volume:
set initial value
by initial condition.
for (t = beginning_time; t <= ending_time; t += Δt)
for each volume:
compute quantities at interfaces
by piecewise polynomial function.
for each interface:
compute flux at interface
by Riemann solver.
for each volume:
update the solution
by recurrence.
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
一些古怪的名詞
數值方法系列。第二行是作者原創名稱。
Godunov's method with slope limiter
reconstruct–evolve–average method (REA method)
LeVeque
Godunov's method with slope limiter
monotonic upstream-centered scheme for conservation laws (MUSCL)
van Leer
Godunov's method with shape function
discontinuous Galerkin method (DG method)
Cockburn
Godunov's method with higher-order shape function
arbitrary high order using derivatives (ADER-DG method)
Toro
higher-order accuracy numerical scheme
high-resolution scheme
大於一階準度的數值方案們。例如MUSCL、DG method。
Harten
[linear algebra perspective]
(1) input x is a vector.
output too.
x = ⎡ p ⎤ y = ⎡ p' ⎤
⎣ q ⎦ ⎣ q' ⎦
(2) transformation A is a matrix.
A = (∇F)ᵀ
y = A x
(3) area preservation
(a) determinant
⎡ | | ⎤ ⎡ | | ⎤
det⎢ x₁ x₂ ⎥ = det⎢ Ax₁ Ax₂ ⎥
⎣ | | ⎦ ⎣ | | ⎦
(b) cross product
x₁ × x₂ = Ax₁ × Ax₂ (add a padding zero)
(c) bilinear form
x₁ᵀ J x₂ = (Ax₁)ᵀ J (Ax₂)
where J = ⎡ 0 1 ⎤
⎣ -1 0 ⎦
[exterior algebra perspective]
(1) input (p, q) is a 2-form.
output too.
x = (p, q)
y = (p', q')
(2) transformation F
y = F(x)
(3) area preservation
(a) wedge product
p Λ q = p' Λ q'
[Cheng–Knorr 1976]
The integration of the vlasov equation in configuration space
https://www.osti.gov/servlets/purl/4200114
(6) and (7)
dv/dt = -E(x)
dx/dt = v
f⁽⃰⁾(x,v) = f⁽ⁿ⁾(x-vΔt/2, v)
f⁽⃰⃰⁾(x,v) = f⁽⃰⁾(x, v+E(x)Δt)
f⁽ⁿ⁺¹⁾(x,v) = f⁽⃰⃰⁾(x-vΔt/2, v)
[Charney–Eliassen 1949]
A Numerical Method for Predicting the Perturbations of
the Middle Latitude Westerlies
(7)(8)相差一個負號，與眾不同。導致旋度公式不成立。
[Charney–Fjørtoft–von Neumann 1950]
Numerical Integration of the Barotropic Vorticity Equation
(I3)上方提到Laplacian和Jacobian。
Jacobian定義成交叉微分項，與眾不同。其實是determinant of Jacobian。
[Arakawa 1966]
Computational design for long-term numerical integration of
the equations of fluid motion: Two-dimensional incompressible
flow. Part I
(2)(3)(13)沿襲前人的定義。
[Arakawa–Lamb 1980]
A Potential Enstrophy and Energy Conserving Scheme for the
Shallow Water Equations
微分方程式取名為potential vorticity advection equation。
shallow water equations:
⎧ ∂h/∂t + ∂(hu)/∂x + ∂(hv)/∂y = 0 conservation of mass
⎨ ∂u/∂t + u ∂u/∂x + v ∂u/∂y + g ∂h/∂x = 0 advection of velocity
⎩ ∂v/∂t + u ∂v/∂x + v ∂v/∂y + g ∂h/∂y = 0 advection of velocity
shallow water equations (Poisson bracket formulation):
∂f/∂t + {f,H}(h,u) = 0
H ≔ ½hu² + ½hv² + ½gh² energy
f ≔ F(h,u,v) height and velocity
p ≔ h height
q ≔ u x velocity
懸案：守恆律可以強行轉換成Hamilton's equations，然後套用扭對稱法。【尚待確認】
Euler–Lagrange equation
我不知道這些主題對於設計數值方案有什麼幫助。有些文獻提到這些主題，於是順手紀錄起來。
Eulerian：用來描述格子流動變化。
Lagrangian：用來描述粒子移動軌跡。
Lagrangian L(q,q̇,t)
generalized coordinates q = x
generalized velocity q̇ = dx/dt
time t
(1) first-order linear differential equation
a₀ f + a₁ f′ = b
approach: integrator
solution: exponentiation
https://tutorial.math.lamar.edu/classes/de/linear.aspx
(2) second-order linear differential equation
a₀ f + a₁ f′ + a₂ f″ = b
approach: characteristic equation
solution: Euler's identity
https://www.ryanjuckett.com/damped-springs/
(3) eigenvalue problem
f′ = k f where k is constant function
approach: characteristic equation
solution: eigenfunction (eigenmode)
https://tutorial.math.lamar.edu/classes/de/bvpevals.aspx
(4) linear differential equation
a₀ f + a₁ f′ + a₂ f″ + ... = b
approach: Laplace transform (spectral method)
solution: fundamental solution
https://tutorial.math.lamar.edu/Classes/DE/IVPWithLaplace.aspx
https://tutorial.math.lamar.edu/Classes/DE/Laplace_Table.aspx
https://tutorial.math.lamar.edu/classes/de/FundamentalSetsofSolutions.aspx
(5) system of linear differential equation
f⃗′ = A f⃗ + b⃗
approach: eigendecomposition --> Laplace transform
solution: fundamental solution
https://tutorial.math.lamar.edu/classes/de/SystemsIntro.aspx
微分方程式的化簡
(1) principle of superposition
疊加原理。強行分開邊界條件，成為多個問題，分別求解，疊加答案。
(2) separation of variables
分離變數法。強行分開變數，成為多個問題，分別求解，融合答案。
(3) method of lines
線法。強行將變數離散化，化作常微分方程式。僅適用IVP。
(4) method of characteristics
特徵線法。強行將變數對時間偏微分，化作常微分方程組。僅適用IVP。
各種生醫理工科系（除了資工系）都會開設專業課程，講解自身領域經常使用的微分方程式。舉例來說，淺水方程組shallow water equations在大氣系、海洋系。電磁方程組Maxwell's equations在電機系。柔體方程組Navier–Cauchy equations在機械系、土木系。
經典的微分方程式：exponentiation
基礎的微分方程式。
exponential equation ḟ = f fₜ = f
phase rotation equation ḟ = 𝑖f fₜ = 𝑖f
0D wave equation f̈ = f fₜₜ = f
oscillator equation f̈ = -f fₜₜ = -f
exponential equation
「指數方程式」。
ḟ = f，微分之後還是一樣。微分不動點。
ḟ = cf ，添上倍率c。c>1是成長率，c<1是衰減率。
符號解是指數函數。
f(t) = f(0)exp(ct)
phase rotation equation
「相位旋轉方程式」。
ḟ = 𝑖f，添上虛數𝑖。
ḟ = 𝑖ωf ，添上角速度ω。
符號解是歐拉公式。
f(t) = f(0)exp(𝑖ωt) = f(0)(cos(ωt) + 𝑖sin(ωt))
0D wave equation【尚無正式名稱】
「零維波動方程式」。
⎰ ṗ = cq
⎱ q̇ = cp
一式對時間微分、替換二式，簡化為一道式子。
f̈ - c²f = 0
符號解是雙曲函數。A和B取決於初始條件。
f(t) = A exp(ct) + B exp(-ct) = C cosh(ct) + D sinh(ct)
A + B = f(0) initial position
A - B = ḟ(0)/c initial velocity
A = (f(0) + ḟ(0)/c) / 2
B = (f(0) - ḟ(0)/c) / 2
C = f(0) initial position
D = ḟ(0)/c initial velocity
oscillator equation
「振子方程式」。簡諧運動。扭對稱動態系統經典範例。
⎰ ṗ = -ωq
⎱ q̇ = +ωp
一式對時間微分、替換二式，簡化為一道式子。形成守恆律。
f̈ + ω²f = 0
符號解是弦波。A和B取決於初始條件。
f(t) = A cos(ωt) + B sin(ωt) = C cos(ωt - ϕ)
A = f(0) initial position
B = ḟ(0)/ω initial velocity
C = magnitude(A+B𝑖) = sqrt(A²+B²) initial magnitude
ϕ = phase(A+B𝑖) = tan⁻¹(B/A) initial phase
mg: gravity (force) vector
p: pressure (force per area) scalar
σ: stress (force per area) matrix
ρg: force density from gravity vector
∇p: force density from pressure vector
∇∙σ: force density from stress vector
∫ { ∂/∂t ρ + ∇∙(ρv) } dz = 0 integral
∂/∂t (ρh) + ∇∙(ρhv) = 0 ∫ρdz = ρh
∂/∂t h + ∇∙(hv) = 0 ρ is constant
第二式與第三式可以推導出知名公式p = ρgh。
ρ D/Dt v + ∂/∂z p + ρg = 0 z-axis of advection of velocity
D/Dt v = 0 hydrostatic balance
∂/∂z p + ρg = 0 subtraction
∂/∂z p = -ρg hydrostatic pressure distribution
p = -ρgz ρ and g are constants
p = ρgh unit_vector(h) = - unit_vector(z)
利用p = ρgh改寫第二式。
ρ D/Dt vₓ + ∂/∂x p + ρgₓ = 0 x-axis of advection of velocity
ρ D/Dt vₓ + ∂/∂x p = 0 gₓ = 0
ρ D/Dt vₓ + ∂/∂x (ρgh) = 0 p = ρgh
ρ D/Dt vₓ + ρg ∂/∂x h = 0 ρ and g are constants
D/Dt vₓ + g ∂/∂x h = 0 divides ρ
ρ D/Dt v + ∂/∂y p + ρg = 0 y-axis of advection of velocity
D/Dt v + g ∂/∂y h = 0
z軸0 + 0 = 0沒有意義。精簡成二維空間。
ρ D/Dt v + ∂/∂z p + ρg = 0 z-axis of advection of velocity
0 + 0 = 0
shallow water equations
淺水方程組有兩種寫法。接著介紹動量密度連續方程式。
⎰ ∂/∂t h + ∇∙(hv) = 0 conservation of mass
⎱ ∂/∂t (hv) + ∇∙(hv⊗v) + hg ∇h = 0 conservation of momentum
(1) isobaric process: p = constant
(2) isochoric process: V = constant
(3) isothermal process: T = constant => pρ = constant
(4) isentropic process: S = constant => adiabatic and reversible
(5) adiabatic process: Q = 0 => pρ⁻ᵞ = constant
多變過程是這些過程的通用公式：等壓過程n=0，等容過程n=∞，等溫過程n=1，等熵過程n=γ。
polytropic process: pVⁿ = constant
heat transfer formula
順便介紹一些物理量。
專著《Thermodynamics: An Engineering Approach》。
總能＝內能＋機械能 total energy E = U + M
內能＝顯熱＋潛熱＋鍵能＋核能 internal energy U = ......
機械能＝動能＋位能＋壓力能 mechanical energy M = ½mv² + mgh + pV
內能 internal energy U
外功 external work W
焓＝熱能總和＝內能＋壓力能 enthalpy H = U + pV
熱＝熱能差異＝內能差異＋外功 heat Q = ΔU + W
機械能觀點：能量＝壓力Ｘ體積 p-V diagram energy = pV
熱能觀點：能量＝溫度Ｘ熵 T-S diagram energy = TS
r: specific gas constant
γ: heat capacity ratio
cₚ: specific heat capacity at constant pressure
cᵥ: specific heat capacity at constant volume
formulas of ideal gas：等容過程的一些關係式。用於Euler's equation，將壓力p改寫成密度ρ與比內能eᵢ，減少未知函數。
⎰ pV = mrT ideal gas law
⎱ Qᵥ = cᵥmΔT heat transfer at constant volume
=> ΔU = Qᵥ = cᵥmΔT = (r/(γ-1))mΔT = ΔpV/(γ-1)
=> Δeᵢ = ΔU/m = Δp/ρ(γ-1)
=> eᵢ = p/ρ(γ-1)
=> ρeᵢ = p/(γ-1)
=> p = (γ−1)ρeᵢ
eᵢ: specific internal energy
經典的微分方程式：electrodynamics
引入了散度、旋度。
electric force:
Coulomb's force
Lorentz's force
Abraham–Lorentz force
electrostatics:
Gauss's law
Coulomb's law
magnetostatics:
Ampère's law
Biot–Savart law
Ampère's force law
electromagnetic induction:
Ampère–Maxwell law
Faraday's law
electric energy:
Faraday's law
electrical resistivity and conductivity:
Ohm's law
Gauss's law:
∇∙E = ρ/ε₀ define ρ = q/V = ∂q/∂V (V is constant)
l length 長度 m
A area 面積 m²
V volume 體積 m³
q charge 電荷 C
ρ charge density 電荷密度 C/m³
F electric force 電力＝力 N
E electric field 電場＝電力／電荷 N/C
F∙A 電力通量 Nm²
E∙A electric flux 電通量 Nm²/C
∇∙E 電場散度 N/Cm²
ε₀ permittivity 真空電容率 Nm²/C²
of free space
Coulomb's law (force between two point charges):
F = kq₁q₂r̂/‖r‖² where k = 1/(4πε₀)
F: electric force 電力＝力
q: electric charge 電荷
x: position 位置
r: displacement 移位
π: 3.14159... 圓周率
ε₀: permittivity of free space 真空電容率
或者簡單來說：球面的電力通量，在球心設置一個假想電荷。
FA = (1/ε₀)q² Gauss's law
∮FdA = F(4πr²) = (1/ε₀)q² surface integral over sphere
F = (1/ε₀)q²/(4πr²) = kq²/r² define k = 1/(4πε₀)
F = kq²r̂/‖r‖² scalar r → vector r
F = kq₁q₂r̂/‖r‖² test charge at sphere center
l length 長度 m
A area 面積 m²
V volume 體積 m³
q charge 電荷 C
ρ charge density 電荷密度 C/m³
j current density 電流密度 (C/m³)(m/s)
qB ／人◕ ‿‿ ◕人＼ N/(m/s)
B magnetic field 磁場 (N/C)/(m/s)
A×B 面積與磁場的叉積 (Nm²/C)/(m/s)
∇×B 磁場旋度 (N/Cm²)/(m/s)
μ₀ permeability 真空磁導率 Ns²/C²
of free space
Ampère's force law (force between two parallel currents):
F/L = μ₀I₁I₂/2πd
ρ: charge density 電荷密度
j: current density 電流密度
x: position 位置
r: displacement 移位＝位置－位置
I: electric current 電流＝電流密度Ｘ面積
π: 3.14159... 圓周率
μ₀: permeability of free space 真空磁導率
μ: electron mobility 電子遷移率
v: drift velocity 漂移速度＝速度
j: current density 電流密度＝每單位體積的電荷速度
ρ: charge density 電荷密度＝每單位體積的電荷
n: charge carrier density 載子密度＝每單位體積的載子
q: charge per carrier 每單位載子的電荷（不是之前的q）
經典的微分方程式：magnetohydrodynamics
集大成。
Maxwell's equations 電磁場
electromagnetic wave equation 電磁波
equation of motion of charged particle 電磁砲
induction equation 電磁牆
magnetohydrodynamic equations 電磁流
F = Fele + Fmag
F = qEele + qEmag define E = F/q
F = qE + q(v × B) Lorentz's law: Emag = v × B
F = q (E + v × B)
f = ρ (E + v × B) define f = F/V and ρ = q/V
f = ρE + (j × B) define j = ρv
∆B: magnetic diffusion 磁擴散
η: magnetic diffusivity 磁擴散率
原理：帶電粒子運動方程式、歐姆定律、安培定律、法拉第定律，四式聯立。
方便起見，捨棄安培－馬克士威定律，採用安培定律。省略的那一項，數值相對較小，無傷大雅。
⎧ F = q (E + v × B) equation of motion of charged particle
⎨ j = σ E Ohm's law
⎪ ∇×B = μ₀j Ampère's law
⎩ ∇×E = - ∂/∂t B Faraday's law
F = q (E + v × B) equation of motion of charged particle
Eeff = E + v × B effective electric field
j = σ (E + v × B) apply Ohm's law
∇×B = μ₀σ (E + v × B) apply Ampère's law
E = (1/μ₀σ) ∇×B - v × B
∇×((1/μ₀σ) ∇×B - v × B) = - ∂/∂t B apply Faraday's law
∂/∂t B + (1/μ₀σ) ∇×∇×B - ∇×(v×B) = 0
∂/∂t B - (1/μ₀σ) ∆B - ∇×(v×B) = 0 since ∇×∇×B = -∆B
F = q (E + v × B) equation of motion of charged particle
Eeff = E + v × B effective electric field
j = σ (E + v × B) apply Ohm's law
j/σ = E + v × B
0 = E + v × B perfect conductor: σ→∞
Lorentz force density:
j×B
force density from pressure:
∇(½(B∙B)/μ₀)
force density from stress:
(B∙∇B)/μ₀
force density from stress (under magnetostatics):
(B∙∇B)/μ₀ = ∇∙(B⊗B)/μ₀ when ∇∙B = 0
magnetic pressure:
½(B∙B)/μ₀
magnetic stress:
(B⊗B)/μ₀
p: total pressure
p₀: static pressure (equilibrium pressure)
p̃: acoustic pressure (pressure perturbation)
ρ: total density
ρ₀: static density (equilibrium density)
ρ̃: density perturbation
c: propagation speed in medium
γ: ratio of specific heat capacity
ₘ: membrane 細胞膜
ᴋ: potassium ion channel 鉀離子通道
ɴₐ: sodium ion channel 鈉離子通道
ₗ: leakage ion channel 洩漏離子通道
I: electric current 電流
V: electric potential 電位
C: capacitance per unit area 電容率
g: conductance per unit area 電導率
n: potassium channel subunit activation 鉀離子通道活化量
m: sodium channel subunit activation 鈉離子通道活化量
h: sodium channel subunit inactivation 鈉離子通道不活化量
α: forward transport rate 正向傳輸率
β: backward transport rate 逆向傳輸率
.jpg)
