function

陣列：各個整數位置，各有一個數值。

函數：可以想成是陣列的連續版本。

Plot3D[Sin[x * y] * Cos[x + y], {x, -Pi, Pi}, {y, -Pi, Pi}, PlotRange -> {-2, 2}, Boxed -> False, Axes -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &) ]

function運算

函數就和實數、複數一樣，有著各種運算。以下將介紹函數的運算：求值、加、減、乘、除、微、積。

       operation (noun)      operation (verb)    result (noun)
-----  --------------------  ------------------  --------------
( )    evaluation      求值  evaluate      求值  value       值
+      addition        加法  add           加    sum         和
-      subtraction     減法  subtract      減    difference  差
       multiplication  乘法  multiply      乘    product     積
/      division        除法  divide        除    quotient    商
d/dx   differentiation 微分  differentiate 微    derivative  導
∫ dx   integration     積分  integrate     積    integral    積分

求值

給定輸入數值，計算輸出數值。

f(x)      輸入變數剛好一個

f(x, y)   輸入變數剛好兩個

// sin x cos y // f(x, y) = ----- + ----- // y x float f(float x, float y) { // x和y不能是0 return sin(x) / y + cos(y) / x; }

加減乘除

相同的輸入數值，其輸出數值相加，得到新函數。

各處各自疊加，各處互不干涉。

f + g

float f(float x) {return sin(x);} float g(float x) {return exp(x);} float add(float x) { return f(x) + g(x); }

微分

相鄰數字差，通通除以dx，得到新函數。

本站文件「sequence」介紹的是離散版本。此處介紹的是連續版本，只多了個dx：一個無限微小、略大於零的數值。

 d
—— f   輸入變數剛好一個
dx

 ∂
—— f   輸入變數多於一個
∂x

當輸入變數只有一個，導數是座標(x,f(x))的「斜率slope」。當輸入變數有許多個，各個輸入變數分別求得斜率，合稱「梯度gradient」。

float f(float x) {return sin(x) / x;} // x不能是0 float differentiate(float x) { // 令dx是0.001 float dx = 1e-3; // 理論上dx越小則導數越準確， // 實際上則必須考慮浮點數精確度。 return (f(x) - f(x - dx)) / dx; // slope } float f(float x, float y) { return sin(x) / y + cos(y) / x; } struct Vector2D {float x, y;}; Vector2D differentiate(float x, float y) { static float dx = 1e-3, dy = 1e-3; return (Vector2D){ // gradient (f(x, y) - f(x - dx, y)) / dx, (f(x, y) - f(x, y - dy)) / dy }; } Vector2D differentiate(float x, float y) { static float d = 1e-3; float fxy = f(x, y); return (Vector2D){ (fxy - f(x - d, y)) / d, (fxy - f(x, y - d)) / d }; }

積分

從負無限大開始的連續數字和，通通乘以dx，得到新函數。

∫ f dx

當輸入變數只有一個，積分是左至-∞、右至x、下至0、上至f(x)，四個邊界所包圍的「面積area」，面積可正可負。當輸入變數有許多個，各個輸入變數一齊累計，得到「容積volume」。

函數積分最簡單的演算法是rectangle rule：按照定義來，將面積切割成數條寬度為dx的矩形。

然而，左至負無限大，演算法永不結束，怎麼辦？解決方式是增設左邊界，想訂多少就多少。數學家把「自訂左右邊界的積分運算」的結果叫做「定積分definite integral」。

對計算學家來說，積分就是累積和，定積分就是區間和啦，就這樣而已。

float f(float x) {return log10(x);} float integrate(float a, float b) { float dx = 1e-3; int n = floor((b - a) / dx); float integral = 0; for (int i=0; i<n; ++i) // 矩形面積等於長乘寬 integral += f(a + dx*i) * dx; // 剩下的一小段 integral += f(a + dx*n) * (b - (a + dx*n)); return integral; }

矩形畢竟不是無限薄。當函數是斜線，仍有縫隙。只好改用trapezoidal rule：將面積切割成數條寬度為dx的梯形。

float trapezoid(float a, float b) { // 梯形面積等於上底加下底括號乘以高除以二 return (f(a) + f(b)) * (b-a) / 2; } float integrate(float a, float b) { float dx = 1e-3; int n = floor((b - a) / dx); float integral = 0; for (int i=0; i<n; ++i) integral += trapezoid(a + dx*i, a + dx*(i+1)); integral += trapezoid(a + dx*n, b); return integral; }

當函數是曲線，仍有縫隙。只好改用parabolic rule又稱Simpson's rule：梯形的斜邊改成拋物線。

float simpson(float a, float b) { // 此面積公式的由來，請自行搜尋資料，此處不解釋。 float c = (a + b) / 2; return (f(a) + 4*f(c) + f(b)) * (b-a) / 6; } float integrate(float a, float b) { float dx = 1e-3; int n = floor((b - a) / dx); float integral = 0; for (int i=0; i<n; ++i) integral += simpson(a + dx*i, a + dx*(i+1)); integral += simpson(a + dx*n, b); return integral; }

長方形的邊是零次、梯形的斜邊是一次、拋物線是二次。當函數是三次曲線，仍有縫隙。無論如何修正，只要函數是更高次曲線就仍有縫隙，況且運算量也更大了。只好改用adaptive Simpson's rule：dx最初是b-a；當縫隙太大，就讓dx變成一半。

有些區間已經幾乎沒有縫隙，大可不必再切半。

函數起起伏伏，難以估量縫隙大小。只好比較前後算得的面積，當面積差異足夠小，就推定面積已經足夠準確、足夠穩定了。

float simpson(float a, float b) { float c = (a + b) / 2; return (f(a) + 4*f(c) + f(b)) * (b-a) / 6; } float asr(float a, float b, float A, float eps) { float c = (a + b) / 2; float L = simpson(a, c); // dx變一半了 float R = simpson(c, b); // dx變一半了 // 此誤差公式的由來，請自行搜尋資料，此處不解釋。 if (fabs(L+R-A) <= 15*eps) // 前後面積差異足夠小 return L+R+(L+R-A)/15; return asr(a, c, L, eps/2) // 前後面積差異仍太大 + asr(c, b, R, eps/2); } float integrate(float a, float b) { // dx最初是b-a return asr(a, b, simpson(a, b), 1e-3); } float simpson(float a, float b, float fa, float fb, float fc) { return (fa + fb + 4*fc) * (b-a) / 6; } // 減少f(x)的呼叫次數 float asr(float a, float b, float c, float fa, float fb, float fc, float A, float eps) { float d = (a+c) / 2, e = (c+b) / 2; float fd = f(d), fe = f(e); float L = simpson(a, c, fa, fc, fd); float R = simpson(c, b, fc, fb, fe); if (fabs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15; return asr(a, c, d, fa, fc, fd, L, eps/2) + asr(c, b, e, fc, fb, fe, R, eps/2); } float integrate(float a, float b) { float c = (a + b) / 2; float fa = f(a), fb = f(b), fc = f(c); return asr(a, b, c, fa, fb, fc, simpson(a, b, fa, fb, fc), 1e-3); }

除了即時調整dx，亦得事前確立dx。Gauss quadrature：當函數內急外緩，則讓dx內密外疏。使用特殊函數自訂dx。

最後介紹MCMC integration：使用「rejection sampling」，以函數曲線下方的取樣次數比例，推定面積。取樣時，自訂x的分布，可達到自訂dx的效果。

UVa 1356 1280 12528 ICPC 3001

multivariate function

多個函數。

⎧ f₁(x, y) = x + y
⎨ f₂(x, y) = 2x² + 1
⎩ f₃(x, y) = 1 / y

數學家重新解讀為：一個函數，輸入與輸出推廣成多個數值。微觀是多個函數，宏觀是一個函數。

F(x, y) = (x + y, 2x² + 1, 1 / y)

多個數值甚至組成向量、矩陣。多個數值視作一個元件，佯裝成普通的函數。寫來簡便，讀來卻傷腦筋。

             ⎡  x + y  ⎤
F( ⎡ x ⎤ ) = ⎢ 2x² + 1 ⎥
   ⎣ y ⎦     ⎣  1 / y  ⎦

   F(x⃗)    =      y⃗

原本是「多個函數」，卻被重新解讀，稱作「多變量函數」。

multivariate function可以解讀成「空間處處有數據」

空間處處有數據。（圖例只畫其中幾處。）

空間座標是函數輸入，數據是函數輸出。

空間可以是一維直線、二維平面、三維空間、……。

數據可以是一個純量（數值）、一個向量、一個矩陣、……。

空間處處有數據，物理學家稱為「場field」。

例如磁場：空間處處有磁力數據。空間是三維空間，數據是一個三維向量。磁力大小是向量長度、磁力方向是向量角度。

multivariate function可以畫成圖形

純量場：畫成透視圖、濃度圖、等高線圖。等距取樣之後：畫成長條圖、折線圖、點陣圖、數值陣列。

（圖例是二維空間，省略座標軸。）

Plot3D[Cos[x]+Cos[y], {x, -2Pi, 2Pi}, {y, -2Pi, 2Pi}, Axes -> False, Boxed -> False, Mesh -> None, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &) ] DensityPlot[Cos[x]+Cos[y], {x, -2Pi, 2Pi}, {y, -2Pi, 2Pi}, Axes -> False, Frame -> False, MaxRecursion -> 1, PlotPoints -> 5, ColorFunction -> "SolarColors"] ContourPlot[Cos[x]+Cos[y], {x, -2Pi, 2Pi}, {y, -2Pi, 2Pi}, Axes -> False, Frame -> False]

向量場：畫成軌跡圖。等距取樣之後：畫成箭矢圖。

VectorPlot3D[{-y, x, 0.5 * z}, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Boxed -> False, Axes -> False, VectorPoints -> 8, VectorScale -> {Small, Scaled[0.5]}, VectorStyle -> {Black}]

矩陣場：我畫不出來。

分開觀察每個維度，就能勉強畫出來。向量場可以拆成多個純量場。矩陣場可以拆成多個向量場、甚至多個純量場。

multivariate function可以描述世界

物理學家習慣以多變量函數描述物理現象，諸如水流、氣流、熱傳導、電磁場、重力場、聲波、震波、……。精彩的思想突破有：

Faraday以三維場解釋電磁現象，衍生古典場論。

Einstein以四維場解釋重力時間現象，衍生相對論。

Schrödinger以複數四維場解釋波粒現象，衍生量子場論。大家仍在驗證。質量擁有實部虛部，若隱若現，非常奇葩。

想要用計算機模擬物理現象，首先必須瞭解多變量函數！

multivariate function運算

除了原本的加減乘除微積，又多了幾種運算：散、旋、梯、梯散、梯旋。數學家並未命名運算名稱，只命名運算結果。

    	operator  result (noun)
----	--------  --------------------------------------------
∇∙  	div       divergence  散度
∇×  	curl      curl        旋度
∇   	grad      gradient    梯度
∇∙∇ 	divgrad   Laplacian   梯度的散度（運算符號從右往左疊）
∇×∇ 	curlgrad  ---         梯度的旋度（總是0，缺乏討論意義）

註：∇∙∇經常簡寫成∇²或∆，稱作Laplace運算子。
　　Laplace運算子(U+2206)、大寫希臘字母delta(U+0394)，兩者是相異字元。
　　然而當今所有字體都沒有特地設計Laplace運算子的造型，
　　而是直接複製大寫希臘字母delta的造型，導致兩者外觀一模一樣。
　　https://tex.stackexchange.com/questions/76553/

各種運算如圖所示：

散度

向量場，求散度，得純量場。模仿了點積。

                dFx
div(F) = ∇∙F = ———               1D
                dx 

                ∂Fx   ∂Fy
div(F) = ∇∙F = ——— + ———         2D
                ∂x    ∂y 

                ∂Fx   ∂Fy   ∂Fz
div(F) = ∇∙F = ——— + ——— + ———   3D
                ∂x    ∂y    ∂z

向量場拆成XY兩個分量。觀察座標(x,y)。（三維很難畫，因此圖例為二維。）

X分量的左右差異、Y分量的上下差異，引發了伸縮。

分別除以dx與dy，加總兩個差異。以朝外為正值。

物理意義是每一處的縮放多寡。

實際應用，例如推擠、流動。

struct Vector2D {float x, y;}; // vector field Vector2D F(float x, float y) { return (Vector2D){x * x, x * y}; } float divergence(float x, float y) { static float dx = 1e-3, dy = 1e-3; return (F(x, y).x - F(x - dx, y).x) / dx + (F(x, y).y - F(x, y - dy).y) / dy; }

旋度

向量場，求旋度，得向量場。模仿了叉積。

curl(F) = ∇×F = 0                                       1D

                           ∂Fy   ∂Fx         ∂Fy   ∂Fx
curl(F) = ∇×F = ( 0 , 0 , ——— - ——— )  ==>  ——— - ———   2D
                           ∂x    ∂y          ∂x    ∂y 

                   ∂Fz   ∂Fy   ∂Fx   ∂Fz   ∂Fy   ∂Fx  
curl(F) = ∇×F = ( ——— - ——— , ——— - ——— , ——— - ——— )   3D
                   ∂y    ∂z    ∂z    ∂x    ∂x    ∂y   

Y分量的左右差異、X分量的上下差異，引發了轉動。

分別除以dx與dy，加總兩個差異。以逆時針為正值。

物理意義是每一處的旋轉多寡。

二維旋度是XY平面上的旋轉。三維旋度則是分別於YZ、ZX、XY三個平面上的旋轉。

實際應用，例如渦度。

struct Vector2D {float x, y;}; // vector field Vector2D F(float x, float y) { return (Vector2D){x * x, x * y}; } float curl(float x, float y) { static float dx = 1e-3, dy = 1e-3; return (F(x, y).y - F(x - dx, y).y) / dx - (F(x, y).x - F(x, y - dy).x) / dy; }

散度定理、旋度定理

任意向量場、任意封閉區域。觀察邊界內、邊界上。

散度定理：邊界內的向量的散度的總和、邊界上的向量的法線分量的總和（通量），兩者相等。

旋度定理：邊界內的向量的旋度的總和、邊界上的向量的切線分量的總和（環量），兩者相等。

原理是(b-a)+(c-b)+(d-c)+...+(z-y) = z-a。導數的區間和＝尾項減去頭項。

梯度

純量場，求梯度，得向量場。

               df
grad(f) = ∇f = ——                 1D
               dx

                 ∂f   ∂f  
grad(f) = ∇f = ( —— , —— )        2D
                 ∂x   ∂y  

                 ∂f   ∂f   ∂f  
grad(f) = ∇f = ( —— , —— , —— )   3D
                 ∂x   ∂y   ∂z  

物理意義是相鄰位勢差距。

純量場與梯度向量場，物理學家稱作「勢potential」與「梯度場gradient field」。

例如電勢與電場。電勢的梯度是電場。電勢是一個數值，代表正電荷多寡。正電荷四散，從電勢高往電勢低，形成了電流。空間處處有電流，形成了電場。

struct Vector2D {float x, y;}; // scalar field float f(float x, float y) { return x * y; } Vector2D gradient(float x, float y) { static float dx = 1e-3, dy = 1e-3; return (Vector2D){ (f(x, y) - f(x - dx, y)) / dx, (f(x, y) - f(x, y - dy)) / dy }; }

梯度的散度【目前稱作Laplace operator】

純量場，先求梯度，再求散度，得純量場。

                  d²Fx
divgrad(F) = ∆F = ————                 1D
                  dx² 

                  ∂²Fx   ∂²Fy
divgrad(F) = ∆F = ———— + ————          2D
                  ∂x²    ∂y² 

                  ∂²Fx   ∂²Fy   ∂²Fz
divgrad(F) = ∆F = ———— + ———— + ————   3D
                  ∂x²    ∂y²    ∂z²

物理意義有兩種等價的解釋方式：

一、考慮相鄰差距，差距大小產生對應流量，每一處的瞬間升降多寡，discharge。

二、無視相鄰差距，各處同時均勻四散、離開原處、前往隔壁，每一處的瞬間升降多寡，diffuse。

const float dx = 1e-3, dy = 1e-3; // scalar field float f(float x, float y) { return x * y; } float diffx(float x, float y) { return (f(x, y) - f(x - dx, y)) / dx; } float diffxx(float x, float y) { return (diffx(x, y) - diffx(x - dx, y)) / dx; } float diffy(float x, float y) { return (f(x, y) - f(x, y - dy)) / dy; } float diffyy(float x, float y) { return (diffy(x, y) - diffy(x, y - dy)) / dy; } float divgrad(float x, float y) { return diffxx(x, y) + diffyy(x, y); } float divgrad(float x, float y) { static float dx = 1e-3, dy = 1e-3; return (f(x, y) - f(x-dx, y) * 2 + f(x-dx-dx, y)) / dx / dx + (f(x, y) - f(x, y-dy) * 2 + f(x, y-dy-dy)) / dy / dy; }

梯度的旋度

一定是0。不存在樓梯幻覺。

// scalar field float f(float x, float y) { return x * y; } float curlgrad(float x, float y) { return 0; }

常數函數、和諧函數

常數函數：相鄰高低差距，處處等於零。數學式子是∇f = 0，梯度處處為零。

換句話說，處處相等。

不生不滅。高者抑之，下者舉之。

和諧函數：相鄰出入差距，處處等於零。數學式子是∆f = 0，梯度的散度處處為零；梯度的旋度本來就是零。

換句話說，處處等於四周平均。

因果輪迴。萬物並作，吾以觀復。

物理意義是靜態平衡、動態平衡。引入時間變數之後，則是慣性運動、簡諧運動。

反梯度

反梯度有無限多個。常數函數的梯度是零；正解加上任意常數函數，仍是正解；如同「不定積分的常數項」。

反梯度可能不存在。向量場存在樓梯幻覺，就沒有反梯度。但是還是可以找到平方誤差最小的場，讓其擁有反梯度。

反梯散【目前稱作Poisson equation】

反梯散有無限多個，而且一定存在。和諧函數的梯度的散度是零；正解加上任意和諧函數，仍是正解。

反梯旋【有點類似Clebsch decomposition】

梯度的旋度一定是零。反梯旋似乎缺乏討論意義，查無資料。

Clebsch decomposition有點類似反梯旋，不過我還沒學會。一個向量場等於兩個梯度場的加權總和F=aX+bY，令其中一個權重是常數函數a=1。分解方式無限多種。

F = grad(x) + b grad(y) = X + b Y

整個式子求旋度。原本向量場的旋度，垂直於b與y的梯度。

curl(F) = grad(b) × grad(y)

找到適當的b與y，旋度可以畫成圖形。請參考《Inside Fluids: Clebsch Maps for Visualization and Processing》。

反散度

反散度有無限多個。散度處處為零的場，稱作「無散場solenoidal field」；正解加上任意無散場，仍是正解。

反旋度【目前稱作Biot–Savart operator】

已知旋度，求得原本的向量場。

反旋度有無限多個。旋度處處為零的場，稱作「無旋場irrotational field」；正解加上任意無旋場，仍是正解。

邊界條件：已知函數值
【目前稱作Dirichlet boundary condition】

上述的反運算系列，正解都有無限多個。為了有戲可唱，於是數學家追加限制條件，以得到唯一解。

經典的限制條件是邊界條件：指定邊界的輪廓形狀（可以無限寬闊），以及指定邊界的每一個函數值多寡。

電腦做運算，數值要限量。邊界條件派上用場。

反梯度：任取兩解，相減必為零。沒有相異解，只有唯一解。

已知梯度G，已知反梯度的邊界函數值b，欲求反梯度f。
⎰ grad(f) = G
⎱ bound(f) = b

任取兩解f₁與f₂。
⎰ grad(f₁) = G    ⎰ grad(f₂) = G
⎱ bound(f₁) = b   ⎱ bound(f₂) = b

兩式相減，抵銷梯度，抵銷邊界。
⎰ grad(f₁)  - grad(f₂)  = grad(f₁ - f₂)  = grad(fdiff)  = 0
⎱ bound(f₁) - bound(f₂) = bound(f₁ - f₂) = bound(fdiff) = 0

已知梯度0，已知反梯度的邊界函數值0，欲求反梯度fdiff。
⎰ grad(fdiff) = 0      梯度是零，即是常數函數。
⎱ bound(fdiff) = 0     邊界是零，又是常數函數，必然處處是零。

得到fdiff = 0。因此不存在相異解、只存在唯一解。

反梯散：最後一步是和諧函數，處處等於四周平均。若邊界是零，則內部不高於零、不低於零、處處是零。

⎰ graddiv(fdiff) = 0   梯度的散度是零，即是和諧函數。
⎱ bound(fdiff) = 0     邊界是零，又是和諧函數，必然處處是零。

反散度：最後一步是旋度定理，邊界上的切線分量總和、邊界內的旋度總和，兩者相等。若邊界是零，則內部旋度總和是零。【尚待確認】

⎰ div(Fdiff) = 0       散度是零，散度處處是零。即是無散場。
⎱ bound(Fdiff) = 0     邊界是零，旋度總和是零。必然有解？

散旋諧分解【目前稱作Helmholtz–Hodge decomposition】

v = (sin(x+y) + 0.5, sin(x-y) - 1) curl = (sin(x) cos(y), -sin(y) cos(x)) div = (sin(y) cos(x), sin(x) cos(y)) har = (0.5, -1) div potential = curl potential = sin(x) sin(y)

一個向量場等於三個向量場相加：散場、旋場、諧場。只有唯一一種分解方式。

F = D + C + H
such that div(D) ≠ 0, curl(D) = 0,
          div(C) = 0, curl(C) ≠ 0,
          div(H) = 0, curl(H) = 0

散場：旋度處處是零
旋場：散度處處是零
諧場：兩者處處是零

散場、旋場，處處互相垂直。

諧場，是和諧函數的梯度。

散度、旋度的重要公式：

1. div(grad(f)) = divgrad(f)
2. curl(grad(f)) = 0
3. curl(curl(F)) = grad(div(F)) - div(grad(F))
4. div(curl(F)) = 0

註：向量場F的梯度的散度：Fx Fy Fz三個純量場各自求梯度的散度。

散場、旋場代入上述公式：

  div(grad(f)) = divgrad(f)  	 curl(curl(F)) = -divgrad(F)

          divgrad            	          -divgrad
    ———————————————————      	     ———————————————————
   |                   |     	    |                   |
   |   grad      div   ↓     	    |   curl      curl  ↓
  散勢 ———> 散場 ———> 散度   	   旋勢 ———> 旋場 ———> 旋度
(純量場)  (向量場)  (純量場) 	 (向量場)  (向量場)  (向量場)

註：旋場使得grad(div(F)) = 0。

散旋諧分解的計算公式：

          divgrad            	          -divgrad
    ———————————————————      	     ———————————————————
   |                   |     	    |                   |
   |   grad      div   ↓     	    |   curl      curl  ↓
  散勢 ———> 散場 ———> 散度   	   旋勢 ———> 旋場 ———> 旋度
(純量場)  (向量場)  (純量場) 	 (向量場)  (向量場)  (向量場)
             ↑         ↑     	              ↑         ↑
          DCH|         |equal	           DCH|         |equal
             |   div   |     	              |   curl  |
             場  ———> 散度   	              場  ———> 旋度

⎧ D = grad(divgrad⁻¹(div(F)))
⎨ C = curl(-divgrad⁻¹(curl(F)))
⎩ H = F - D - C

註：令邊界函數值為零，確保divgrad⁻¹有唯一解、不含和諧函數。

散場：求散度、求反梯散、求梯度
旋場：求旋度、三個值分別求反梯散、三個值變號、求旋度
旋場（二維）：求旋度（只得z值）、求反梯散、變號、求梯度、轉90度
諧場：減散場、減旋場

散旋諧分解，也有人寫成勢的風格，而非場的風格。

F = D + C + H = grad(d) + curl(c) + grad(h)

⎧ d = grad⁻¹(D) =  divgrad⁻¹(div(F))
⎨ c = curl⁻¹(C) = -divgrad⁻¹(curl(F))
⎩ h = grad⁻¹(H)

散勢：散場的反梯度
旋勢：旋場的反旋度
旋勢（二維）：旋場轉-90度的反梯度
諧勢：諧場的反梯度

諧場可以隨意分攤給散場、旋場。散場摻加諧場，形成無旋場。旋場摻加諧場，形成無散場。

F = Dirr + Csole
F = grad(dirr) + grad(csole)

散旋諧分解的傅立葉轉換，我還沒有學會。

縱場：散場的傅立葉轉換
橫場：旋場的傅立葉轉換
查無資料：諧場的傅立葉轉換
http://physics.stackexchange.com/questions/1115/

原本的場，平方誤差最小的梯度場，兩者散度一樣。

散場加諧場，正是平方誤差最小的梯度場，擁有反梯度。

          divgrad            
    ———————————————————      
   |                   |     
   |   grad      div   ↓     
  散勢 ———> 散場 ———> 散度   
 +諧勢     +諧場      +0     
  (勢)    (梯度場)           
             ↑         ↑     
         proj|         |equal
             |   div   |     
             場  ———> 散度   

符號、數值

數有兩種記載方式：符號、數值。

               1        |                  
      1 ÷ 3 = ———       | 1 ÷ 3 = 0.333333 
               3        |                  
        symbol          |      numeral     
      橫槓代表分數       |    無法精準記錄 

 f(x) = ⎰ x   if x > 0  | f(0.01) = 0.01  
        ⎱ 0   otherwise | f(0.02) = 0.02  
                        |     :       :   
        symbol          |     numeral     
      以文字描述         |   無法盡數記錄 

採用符號，函數可以寫成程式碼、寫成函式。

// f(x) = x² + 2x + 1 float f(float x) {return x*x + x*2 + 1;} // f(x) = ⎰ x if x > 0 // ⎱ 0 otherwise float f(float x) {return x > 0 ? x : 0;}

採用數值，函數可以儲存於陣列。

// f(x) = x² + 2x + 1 float f[10] = {16, 9, 4, 1, 0, 1, 4, 9, 16, 25}; // f(x) = ⎰ x if x > 0 // ⎱ 0 otherwise float f[10] = {0, 0, 0, 0, 0, 1, 2, 3, 4, 5};

符號轉數值，就是函數求值。數值轉符號，就是函數內插。

符號轉數值，函數必須取樣、擇鄰。稱作「離散化discretization」。

regular discretization

符號轉數值。取樣：等距方格取樣。擇鄰：左右上下前後。

數值轉符號。傾向使用monotone cubic interpolation。

const int X = 8, Y = 8, Z = 8; // 取樣數目 float f[X]; // 一元函數（等間距） float f[X][Y]; // 二元函數（正方形網格） float f[X][Y][Z]; // 三元函數（正方體網格） const int N = 8; // 取樣數目 const float Δx = 1.0; // 取樣間距 const float x0 = 2.0; // 取樣起點 void discretize(float (*function)(float), float f[N]) { // 從x=2開始往右抓8點 float x = x0; for (int i=0; i<N; ++i) { f[i] = function(x); x += Δx; } } // monotone Catmull–Rom spline interpolation float continuize(float f[N], float x) { int i = (x - x0) / Δx; // 邊界外面補充函數點。一次內插。 float fL = f[0] - (f[1] - f[0]); float fR = f[N-1] + (f[N-1] - f[N-2]); // 取四個點 float f0 = (i-1 >= 0) ? f[i-1] : fL; float f1 = (i >= 0) ? f[i] : fL; float f2 = (i+1 < N) ? f[i+1] : fR; float f3 = (i+2 < N) ? f[i+2] : fR; // 一次導數：左右兩點斜率 float d1 = (f2 - f0) / (2 * Δx); float d2 = (f3 - f1) / (2 * Δx); float d = (f2 - f1) / Δx; // Δy/Δx // 單調 if (d == 0) d1 = d2 = 0; // 建議考慮浮點數誤差 if (d * d1 < 0) d1 = 0; // 建議寫成sign函數 if (d * d2 < 0) d2 = 0; // 建議寫成sign函數 // 多項式係數 float c3 = (d * -2 + d1 + d2) / Δx / Δx; float c2 = (d * 3 - d1 * 2 - d2) / Δx; float c1 = d1; float c0 = f1; x -= x0 + Δx * i; // return c3 * x*x*x + c2 * x*x + c1 * x + c0; return ((c3 * x + c2) * x + c1) * x + c0; }

Chebyshev discretization

符號轉數值。取樣：等距圓弧取樣。擇鄰：左右上下前後。

數值轉符號。傾向使用polynomial interpolation。

多項式內插的缺點是Runge phenomenon：左右邊界震盪大，函數曲線與函數點的走向沒有貼合。

等距圓弧取樣，中央疏、左右密，為的就是改善缺點。

專著《Spectral Methods in MATLAB》。函式庫chebfun。

const int N = 8; // 取樣數目 const float r = 1.0; // 取樣範圍的一半 const float c = 6.0; // 取樣範圍的中點 const float π = acos(0) * 2; // 圓周率 void discretize(float (*function)(float), float f[N]) { // 以x=6為中心抓8點 for (int i=0; i<N; ++i) { float x = c + r * cos(π / (N-1) * i); f[i] = function(x); } } // Lagrange interpolation float continuize(float f[N], float x) { float y = 0; for (int i=0; i<N; ++i) { float a = 1, b = 1; for (int j=0; j<N; ++j) { if (j == i) continue; float xi = c + r * cos(π / (N-1) * i); float xj = c + r * cos(π / (N-1) * j); a *= x - xj; b *= xi - xj; } y += f[i] * a / b; } return y; }

meshfree discretization

符號轉數值。取樣：隨機取樣。擇鄰：自訂距離，門檻以內都是鄰居；自訂半徑，圓內都是鄰居。

數值轉符號。傾向使用radial basis function interpolation。

mesh discretization

p = {{6.1,7.5,0.84},{0.95,0.52,0.21},{0.51,7.5,0.61},{7.0,2.1,0.86},{2.3,3.8,0.64},{5.3,5.0,0.99},{3.5,0.49,0.37},{2.1,5.7,0.77},{0.43,3.1,0.34},{4.5,7.7,0.98},{2.9,1.9,0.50},{5.2,0.096,0.41},{0.24,5.6,0.46},{5.0,2.9,0.82},{7.8,5.3,0.92},{0.42,1.7,0.25},{7.5,7.7,0.55},{7.7,0.67,0.66},{3.9,6.0,0.97},{5.8,2.0,0.75},{7.4,3.1,0.97},{0.44,4.3,0.41},{3.2,4.6,0.82},{3.9,2.4,0.66},{6.3,5.4,0.99},{3.0,7.7,0.97},{6.3,3.3,0.94},{4.6,1.6,0.62},{6.1,0.63,0.56},{1.1,5.1,0.58},{2.6,0.95,0.37},{7.0,6.5,0.85},{7.9,4.5,0.98},{7.8,6.7,0.70},{4.2,4.2,0.87},{1.2,6.6,0.68},{5.3,6.7,0.98},{5.1,4.2,0.94},{2.4,6.6,0.87},{0.19,0.020,0.12},{2.0,1.9,0.43},{1.9,4.6,0.66},{1.2,3.7,0.48},{3.6,3.6,0.77},{3.2,2.9,0.65},{2.0,3.1,0.53},{5.0,5.9,1.0},{4.3,0.018,0.34},{6.2,4.7,1.0},{1.7,0.80,0.28},{2.4,0.16,0.26},{1.8,7.6,0.86},{3.2,6.0,0.91},{7.1,5.0,0.98},{3.9,6.9,0.99},{1.2,2.6,0.40},{4.3,5.2,0.96},{2.5,2.9,0.58},{6.3,1.3,0.68},{5.6,5.9,1.0},{4.9,7.1,0.99},{0.23,0.87,0.18},{8.0,2.1,0.91},{2.5,4.4,0.72}}; g1 = ListPointPlot3D[p, PlotStyle -> {PointSize[Large], Black}, Filling -> Bottom, FillingStyle -> Thick, Boxed -> False, Axes -> False]; q = Transpose[{p[[All,1]], p[[All,2]]}]; (*remove z*) m = DelaunayMesh[q]; q = p; q[[All,3]] = 0; (*z:=0*) g2 = MeshRegion[q, MeshCells[m, 1], MeshCellStyle -> {{1, All} -> { RGBColor[79/255,129/255,189/255], Thick}, {0, All} -> None}]; Show[g2, g1]

f[{x_, y_}] = Sin[x * y] * Cos[x + y]; p = Flatten[Table[{i, j}, {i, -3, 3, 0.8}, {j, -3, 3, 0.8}], 1]; q = Transpose[{p[[All,1]], p[[All,2]], N[Map[f, p]]}]; g1 = ListPointPlot3D[q, PlotStyle -> {PointSize[Large], Black}, Filling -> Bottom, FillingStyle -> Thick, Boxed -> False, Axes -> False]; m = DelaunayMesh[p]; g2 = MeshRegion[q, MeshCells[m, 1], MeshCellStyle -> {{1, All} -> {Black, Thick, Opacity[0.5]}, {0, All} -> None}]; Show[g2, g1]

符號轉數值。取樣：隨機取樣。擇鄰：三角剖分。

數值轉符號。傾向使用natural coordinates。

加減乘除

對應位置加減乘除。

void add(float f1[X], float f2[X], float g[X]) { for (int i=0; i<X; ++i) g[i] = f1[i] + f2[i]; }

外插

陣列邊界以外，數值是多少呢？數值轉符號：所有數值做內插，計算量太大。外插：邊界附近取一部分數值做內插，計算量較低。

陣列邊界附近取值，多項式內插，建立延長線，估計數值。取零個是補零，取一個是複製邊界，取兩個是直線，取三個是拋物線。

void extrapolate(float f[X], int i) { // 複製邊界 if (i < 0) return f[0]; if (i >= X) return f[X-1]; return f[i]; }

微分、積分

數列微分a[x] - a[x-1]，函數微分[f(x) - f(x-dx)] / dx，離散化函數微分[f(x+Δx) - f(x-Δx)] / 2Δx，此處談的是最後一種。

無限微小、略大於零的數值dx，改成了取樣間距Δx。

因為捨棄了極限，所以衍生了左邊、右邊、中央三種版本。左邊版本，優點是直觀簡潔，缺點是數值往左偏移。中央版本，優點是數值不會偏移，缺點是微分兩次不會得到二次微分。

一次微分
backward [f(x) - f(x-Δx)] / Δx
 forward [f(x+Δx) - f(x)] / Δx
 central [f(x+Δx) - f(x-Δx)] / 2Δx

二次微分
backward [f(x) - 2f(x-Δx) + f(x-2Δx)] / Δx²
 forward [f(x+2Δx) - 2f(x+Δx) + f(x)] / Δx²
 central [f(x+Δx) - 2f(x) + f(x-Δx)] / Δx²

一次積分
backward [... + f(x-Δx) + f(x)] ⋅ Δx
 forward [... + f(x-2Δx) + f(x-Δx)] ⋅ Δx
 central [... + f(x-3Δx) + f(x-Δx)] ⋅ 2Δx

void differentiate(float f[X], float g[X]) { // f′(x) ≈ [f(x+Δx) - f(x-Δx)] / 2Δx // f[-1] = f[0] - (f[1] - f[0]); // 外插：直線延伸 // g[0] = (f[1] - f[-1]) / 2Δx; // 微分 g[0] = (f[1] - f[0]) / Δx; // 兩式合併 g[X-1] = (f[X-1] - f[X-2]) / Δx; for (int i=1; i<X-1; ++i) g[i] = (f[i+1] - f[i-1]) / (2 * Δx); } void integrate(float f[X], float g[X]) { // ∫f(x) ≈ ∫f(x-2Δx) + f(x-Δx) ⋅ 2Δx g[0] = 0; g[1] = f[0] * Δx; for (int i=2; i<X; ++i) g[i] = g[i-2] + f[i-1] * 2 * Δx; }

散度、旋度、梯度、梯度的散度

都是微分來的。Δx = Δy，得以簡化算式。

// 添加邊框，方便微分 void divergence(float Fx[X+2][Y+2], float Fy[X+2][Y+2], float g[X+2][Y+2]) { extrapolate_boundary(Fx, Fy); // 邊框外插 for (int i=1; i<=X; ++i) for (int j=1; j<=Y; ++j) { g[i][j] = (Fx[i+1][j] - Fx[i-1][j]) + (Fy[i][j+1] - Fy[i][j-1]) g[i][j] /= (2 * Δx); } } void curl(float Fx[X+2][Y+2], float Fy[X+2][Y+2], float g[X+2][Y+2]) { extrapolate_boundary(Fx, Fy); // 邊框外插 for (int i=1; i<=X; ++i) for (int j=1; j<=Y; ++j) { g[i][j] = (Fy[i+1][j] - Fy[i-1][j]) - (Fx[i][j+1] - Fx[i][j-1]); g[i][j] /= (2 * Δx); } } void gradient(float f[X+2][Y+2], float Gx[X+2][Y+2], float Gy[X+2][Y+2]) { extrapolate_boundary(f); // 邊框外插 for (int i=1; i<=X; ++i) for (int j=1; j<=Y; ++j) { Gx[i][j] = f[i+1][j] - f[i-1][j]; Gx[i][j] /= (2 * Δx); Gy[i][j] = f[i][j+1] - f[i][j-1]; Gy[i][j] /= (2 * Δx); } } void divgrad(float f[X+2][Y+2], float g[X+2][Y+2]) { extrapolate_boundary(f); // 邊框外插 for (int i=1; i<=X; ++i) for (int j=1; j<=Y; ++j) { g[i][j] = f[i+1][j] + f[i-1][j] + f[i][j+1] + f[i][j-1] - f[i][j] * 4; g[i][j] /= (Δx * Δx); } }

反梯散（一）【目前稱作Poisson equation】

高斯消去法O(N³)，鬆弛法O(N²K)。

對稱正定矩陣。可以採用專門的演算法。

為了得到唯一解，必須設定邊界條件。邊框外插即可！

1D field. N = 5. zero boundary.
    ⎡ -2  1  0  0  0 ⎤ ⎡ f[1] ⎤   ⎡ g[1] ⎤
 1  ⎢  1 -2  1  0  0 ⎥ ⎢ f[2] ⎥   ⎢ g[2] ⎥
——— ⎢  0  1 -2  1  0 ⎥ ⎢ f[3] ⎥ = ⎢ g[3] ⎥
Δx² ⎢  0  0  1 -2  1 ⎥ ⎢ f[4] ⎥   ⎢ g[4] ⎥
    ⎣  0  0  0  1 -2 ⎦ ⎣ f[5] ⎦   ⎣ g[5] ⎦

1D field. N = 5. duplicate boundary.
    ⎡ -1  1  0  0  0 ⎤ ⎡ f[1] ⎤   ⎡ g[1] ⎤
 1  ⎢  1 -2  1  0  0 ⎥ ⎢ f[2] ⎥   ⎢ g[2] ⎥
——— ⎢  0  1 -2  1  0 ⎥ ⎢ f[3] ⎥ = ⎢ g[3] ⎥
Δx² ⎢  0  0  1 -2  1 ⎥ ⎢ f[4] ⎥   ⎢ g[4] ⎥
    ⎣  0  0  0  1 -1 ⎦ ⎣ f[5] ⎦   ⎣ g[5] ⎦

2D field. N = 4×4. zero boundary.
    ⎡-4 1 0 0|1 0 0 0|0 0 0 0|0 0 0 0 ⎤ ⎡ f[1][1] ⎤   ⎡ g[1][1] ⎤
    ⎢ 1-4 1 0|0 1 0 0|0 0 0 0|0 0 0 0 ⎥ ⎢ f[1][2] ⎥   ⎢ g[1][2] ⎥
    ⎢ 0 1-4 1|0 0 1 0|0 0 0 0|0 0 0 0 ⎥ ⎢ f[1][3] ⎥   ⎢ g[1][3] ⎥
    ⎢ 0 0 1-4|0 0 0 1|0 0 0 0|0 0 0 0 ⎥ ⎢ f[1][4] ⎥   ⎢ g[1][4] ⎥
    ⎢--------+-------+-------+--------⎥ ⎢---------⎥   ⎢---------⎥
    ⎢ 1 0 0 0-4 1 0 0|1 0 0 0|0 0 0 0 ⎥ ⎢ f[2][1] ⎥   ⎢ g[2][1] ⎥
    ⎢ 0 1 0 0|1-4 1 0|0 1 0 0|0 0 0 0 ⎥ ⎢ f[2][2] ⎥   ⎢ g[2][2] ⎥
    ⎢ 0 0 1 0|0 1-4 1|0 0 1 0|0 0 0 0 ⎥ ⎢ f[2][3] ⎥   ⎢ g[2][3] ⎥
 1  ⎢ 0 0 0 1|0 0 1-4|0 0 0 1|0 0 0 0 ⎥ ⎢ f[2][4] ⎥   ⎢ g[2][4] ⎥
——— ⎢--------+-------+-------+--------⎥ ⎢---------⎥ = ⎢---------⎥
Δx² ⎢ 0 0 0 0|1 0 0 0-4 1 0 0|1 0 0 0 ⎥ ⎢ f[3][1] ⎥   ⎢ g[3][1] ⎥
    ⎢ 0 0 0 0|0 1 0 0|1-4 1 0|0 1 0 0 ⎥ ⎢ f[3][2] ⎥   ⎢ g[3][2] ⎥
    ⎢ 0 0 0 0|0 0 1 0|0 1-4 1|0 0 1 0 ⎥ ⎢ f[3][3] ⎥   ⎢ g[3][3] ⎥
    ⎢ 0 0 0 0|0 0 0 1|0 0 1-4|0 0 0 1 ⎥ ⎢ f[3][4] ⎥   ⎢ g[3][4] ⎥
    ⎢--------+-------+-------+--------⎥ ⎢---------⎥   ⎢---------⎥
    ⎢ 0 0 0 0|0 0 0 0|1 0 0 0-4 1 0 0 ⎥ ⎢ f[4][1] ⎥   ⎢ g[4][1] ⎥
    ⎢ 0 0 0 0|0 0 0 0|0 1 0 0|1-4 1 0 ⎥ ⎢ f[4][2] ⎥   ⎢ g[4][2] ⎥
    ⎢ 0 0 0 0|0 0 0 0|0 0 1 0|0 1-4 1 ⎥ ⎢ f[4][3] ⎥   ⎢ g[4][3] ⎥
    ⎣ 0 0 0 0|0 0 0 0|0 0 0 1|0 0 1-4 ⎦ ⎣ f[4][4] ⎦   ⎣ g[4][4] ⎦

// Gauss–Seidel iteration void divgrad⁻¹(float g[X+2][Y+2], float f[X+2][Y+2]) { memset(f, 0, sizeof(f)); // 意思意思遞推5回合，省時間。 for (int k=0; k<5; k++) for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { float s = f[i+1][j] + f[i-1][j] + f[i][j+1] + f[i][j-1] + g[i][j] * Δx * Δx; f[i][j] = s / 4; } } // successive over-relaxation void divgrad⁻¹(float g[X+2][Y+2], float f[X+2][Y+2]) { // successive over-relaxation的權重 float w = 2.0 / (1.0 + 3.1416 / X); float temp[X+2][Y+2]; memset(f, 0, sizeof(f)); // 遞推有限回合，避免發散而永不結束。 for (int k=0; k<100; k++) { // 暫存上回合結果，以便計算每回合的進展。 for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) temp[i][j] = f[i][j]; // successive over-relaxation for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { float s = f[i+1][j] + f[i-1][j] + f[i][j+1] + f[i][j-1] + g[i][j] * Δx * Δx; s /= 4; f[i][j] = f[i][j] * (1-w) + s * w; } // 統計前後回合差距，差距足夠小就立即結束。 float diff = 0; for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) diff += f[i][j] - temp[i][j]; if (diff < 1e-5) break; } }

UVa 199

反梯散（二）【目前稱作Poisson equation】

兩種解法：一次方程組（時域）、傅立葉轉換（頻域）。

循環矩陣＝循環卷積＝頻域乘法。時間複雜度O(NlogN)。

邊界必須循環（邊界條件被迫是零）。邊界輪廓必須是矩形，邊長是2的次方（必須補零）。

傅立葉轉換無法處理其他邊界條件，於是大家發明了各種歪招，我也不清楚是否合乎邏輯。比方說連續傅立葉轉換：

          fourier            |            fourier
     f(x) <—————> g(θ)       |     d/dx f <—————> iθ 2π g
   f(x,y) <—————> g(θ,φ)     |   d²/dx² f <—————> -θ² 4π² g
af₁ + bf₂ <—————> ag₁ + bg₂  |  divgrad f <—————> -(θ² + φ²) 4π² g

   divgrad(f) = b
=> -(θ² + φ²) 4π² fft(f) = fft(b)         等式兩邊同做傅立葉轉換
=> fft(f) = fft(b) / [-(θ² + φ²) 4π²]     移項，注意到是數列除法
=> f = ifft( fft(b) / [-(θ² + φ²) 4π²])   等式兩邊同做逆向傅立葉轉換

註：時域常數項＝頻域首項。
　　反梯散擁有常數項，因此頻域首項可以是任意值，也就不必除以零。

https://math.stackexchange.com/questions/1809871/
https://math.stackexchange.com/questions/1657756/
http://www.damtp.cam.ac.uk/user/naweb/ii/poisson_equation/poisson_equation.php
https://people.eecs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html

反梯度（一）

簡易的做法：切成一條一條，各自積分。

反梯度必須存在，才能這樣做。

void gradient⁻¹(float Gx[X+2][Y+2], float Gy[X+2][Y+2], float f[X+2][Y+2]) { extrapolate_boundary(Gx, Gy); // 邊框外插（補零） memset(f, 0, sizeof(f)); // zero boundary for (int i=1; i<=X; ++i) { if (i == 1) f[i][1] = Gy[i][0] * 2 * Δy; else f[i][1] = (f[i-2][1] + Gy[i-1][1]) * 2 * Δx; for (int j=2; j<=Y; ++j) f[i][j] = (f[i][j-2] + Gx[i][j-1]) * 2 * Δy; } }

反梯度（二）

合理的做法：一次方程組，聯立每個地點的橫向微分與縱向微分。等式太多因而無解，求平方誤差最小的解。

1D field. N = 5. zero boundary.
    ⎡  0  1  0  0  0 ⎤ ⎡ f[1] ⎤   ⎡ g[1] ⎤
 1  ⎢ -1  0  1  0  0 ⎥ ⎢ f[2] ⎥   ⎢ g[2] ⎥
——— ⎢  0 -1  0  1  0 ⎥ ⎢ f[3] ⎥ = ⎢ g[3] ⎥
2Δx ⎢  0  0 -1  0  1 ⎥ ⎢ f[4] ⎥   ⎢ g[4] ⎥
    ⎣  0  0  0 -1  0 ⎦ ⎣ f[5] ⎦   ⎣ g[5] ⎦

2D field. N = 3×3. zero boundary.
    ⎡ 0 1 0 | 0 0 0 | 0 0 0 ⎤ ⎡ f[1][1] ⎤   ⎡ Gy[1][1] ⎤
    ⎢-1 0 1 | 0 0 0 | 0 0 0 ⎥ ⎢ f[1][2] ⎥   ⎢ Gy[1][2] ⎥
    ⎢ 0-1 0 | 0 0 0 | 0 0 0 ⎥ ⎢ f[1][3] ⎥   ⎢ Gy[1][3] ⎥
    ⎢-------+-------+-------⎥ ⎢---------⎥   ⎢----------⎥
 1  ⎢ 0 0 0 | 0 1 0 | 0 0 0 ⎥ ⎢ f[2][1] ⎥   ⎢ Gy[2][1] ⎥
——— ⎢ 0 0 0 |-1 0 1 | 0 0 0 ⎥ ⎢ f[2][2] ⎥ = ⎢ Gy[2][2] ⎥
2Δx ⎢ 0 0 0 | 0-1 0 | 0 0 0 ⎥ ⎢ f[2][3] ⎥   ⎢ Gy[2][3] ⎥
    ⎢-------+-------+-------⎥ ⎢---------⎥   ⎢----------⎥
    ⎢ 0 0 0 | 0 0 0 | 0 1 0 ⎥ ⎢ f[3][1] ⎥   ⎢ Gy[3][1] ⎥
    ⎢ 0 0 0 | 0 0 0 |-1 0 1 ⎥ ⎢ f[3][2] ⎥   ⎢ Gy[3][2] ⎥
    ⎢ 0 0 0 | 0 0 0 | 0-1 0 ⎥ ⎣ f[3][3] ⎦   ⎢ Gy[3][3] ⎥
    ⎢-----------------------⎥               ⎢----------⎥
    ⎢ 0 0 0 | 1 0 0 | 0 0 0 ⎥               ⎢ Gx[1][1] ⎥
    ⎢-1 0 0 | 0 0 0 | 1 0 0 ⎥               ⎢ Gx[1][2] ⎥
    ⎢ 0 0 0 |-1 0 0 | 0 0 0 ⎥               ⎢ Gx[1][3] ⎥
    ⎢-------+-------+-------⎥               ⎢----------⎥
    ⎢ 0 0 0 | 0 1 0 | 0 0 0 ⎥               ⎢ Gx[2][1] ⎥
    ⎢ 0-1 0 | 0 0 0 | 0 1 0 ⎥               ⎢ Gx[2][2] ⎥
    ⎢ 0 0 0 | 0-1 0 | 0 0 0 ⎥               ⎢ Gx[2][3] ⎥
    ⎢-------+-------+-------⎥               ⎢----------⎥
    ⎢ 0 0 0 | 0 0 1 | 0 0 0 ⎥               ⎢ Gx[3][1] ⎥
    ⎢ 0 0-1 | 0 0 0 | 0 0 1 ⎥               ⎢ Gx[3][2] ⎥
    ⎣ 0 0 0 | 0 0 1 | 0 0 0 ⎦               ⎣ Gx[3][3] ⎦

反梯度（三）

常見的做法：先算散度，再算反梯散。

一樣可以得到平方誤差最小的解。證明有點長。

已知G，求F。統計每個地點的平方誤差，令所有地點的誤差總和盡量小：

min ‖ F - G ‖²

min sum  ( Fx(x,y) - Gx(x,y) )² +
   (x,y) ( Fy(x,y) - Gy(x,y) )²


已知G，求f。場F改成了勢f：

min ‖ grad(f)forward - G ‖²

min sum  ( [f(x+1,y) - f(x,y)] / Δx - Gx(x,y) )² +
   (x,y) ( [f(x,y+1) - f(x,y)] / Δy - Gy(x,y) )²


誤差總和e分別對每個未知數f(x,y)偏微分，一次偏微分等於零的地方是極值：

  ∂e          總共三個地點有f(x,y)。
——————— = 0 = -2 ( [f(x+1,y) - f(x,y)] / Δx - Gx(x,y)   ) / Δx
∂f(x,y)     + -2 ( [f(x,y+1) - f(x,y)] / Δy - Gy(x,y)   ) / Δy
            +  2 ( [f(x,y) - f(x-1,y)] / Δx - Gx(x-1,y) ) / Δx
            +  2 ( [f(x,y) - f(x,y-1)] / Δy - Gy(x,y-1) ) / Δy

              重新整理。前兩行是divgrad，後兩行是div。
            = -2 [ f(x+1,y) - 2f(x,y) + f(x-1,y) ] / Δx / Δx
            + -2 [ f(x,y+1) - 2f(x,y) + f(x,y-1) ] / Δy / Δy
            +  2 [ Gx(x,y) - Gx(x-1,y) ] / Δx
            +  2 [ Gy(x,y) - Gy(x,y-1) ] / Δy
           
            = -2 divgrad(f)central (x,y)
            +  2 div(G)backward (x,y)

移項得到：
divgrad(f)central = div(G)backward

平方誤差最小，就是一次微分相等。橫向縱向加總，就是散度相等。梯度的平方誤差最小，就是梯度的散度相等。

右邊版本的梯度、中央版本的梯度的散度、左邊版本的散度，儘管亂七八糟，但是計算學家還是默許了──總之Δx足夠小就行了。

void gradient⁻¹(float Gx[X+2][Y+2], float Gy[X+2][Y+2], float f[X+2][Y+2]) { extrapolate_boundary(Gx, Gy); // 邊框外插（補零） // 散度 float div[X+2][Y+2] = {}; for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { div[i][j] = (Gx[i][j] - Gx[i-1][j]) + (Gy[i][j] - Gy[i][j-1]); div[i][j] /= Δx; } // 反梯散 memset(f, 0, sizeof(f)); // zero boundary for (int k=0; k<5; k++) for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { float s = f[i+1][j] + f[i-1][j] + f[i][j+1] + f[i][j-1] + div[i][j] * Δx * Δx; f[i][j] = s / 4; } }

散旋諧分解、反散度、反旋度

反散度，即是散場。先算反梯散，再算散度。

反旋度，即是旋場。添上負號，先算反梯散，再算旋度。

void DCH2D(float Fx[X+2][Y+2], float Fy[X+2][Y+2], float Dx[X+2][Y+2], float Dy[X+2][Y+2], float Cx[X+2][Y+2], float Cy[X+2][Y+2], float Hx[X+2][Y+2], float Hy[X+2][Y+2]) { extrapolate_boundary(Fx, Fy); // 邊框外插 // 散度 float div[X+2][Y+2] = {}; for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { div[i][j] = (Fx[i][j] - Fx[i-1][j]) + (Fy[i][j] - Fy[i][j-1]); div[i][j] /= Δx; } // 散勢 memset(f, 0, sizeof(f)); for (int k=0; k<5; k++) for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) f[i][j] = (f[i+1][j] + f[i-1][j] + f[i][j+1] + f[i][j-1] + div[i][j] * Δx * Δx ) / 4; // 散場（反散度） for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { Dx[i][j] = (f[i+1][j] - f[i][j]) / Δx; Dy[i][j] = (f[i][j+1] - f[i][j]) / Δy; } // 旋度 float curl[X+2][Y+2] = {}; for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { curl[i][j] = (Fy[i][j] - Fy[i-1][j]) - (Fx[i][j] - Fx[i][j-1]); curl[i][j] /= Δx; curl[i][j] = -curl[i][j]; // 提早變號 } // 旋勢 memset(f, 0, sizeof(f)); for (int k=0; k<5; k++) for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) f[i][j] = (f[i+1][j] + f[i-1][j] + f[i][j+1] + f[i][j-1] + curl[i][j] * Δx * Δx ) / 4; // 旋場（反旋度） for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { Cx[i][j] = (f[i][j+1] - f[i][j]) / Δx; Cy[i][j] = -(f[i+1][j] - f[i][j]) / Δy; } // 諧場 for (int i=1; i<=X; i++) for (int j=1; j<=Y; j++) { Hx[i][j] = Fx[i][j] - Dx[i][j] - Cx[i][j]; Hy[i][j] = Fy[i][j] - Dy[i][j] - Cy[i][j]; } } typedef float ScalarField[X+2][Y+2][Z+2]; typedef ScalarField VectorField[3]; void DCH3D(VectorField F, VectorField D, VectorField C, VectorField H) { ScalarField d, dp; divergence(F, d); divgrad⁻¹(d, dp); gradient(dp, D); VectorField c, cp; curl(F, c); divgrad⁻¹(c[0], cp[0]); divgrad⁻¹(c[1], cp[1]); divgrad⁻¹(c[2], cp[2]); negate(cp); curl(cp, C); VectorField T; subtract(F, D, T); subtract(T, C, H); }

shape function

Lagrange interpolation + delta function。

https://www.fidelisfea.com/post/what-are-shape-functions-in-fea-and-how-are-they-derived
https://teachbooks.tudelft.nl/computational-modelling/introduction/shape.html

basis interpolation

基底函數的規範。一、中央格點是1、其餘格點是0。二、常數場1的內插結果是常數場1。

⎰ Nᵢ(xⱼ) = 1   if i = j
⎱ Nᵢ(xⱼ) = 0ᵢ  if i ≠ j

sum Nᵢ(x) = 1   for all x
 ⁱ

linear basis function:
Nᵢ(x) = ⎰ 1 - (x - xᵢ) / h  , 0 ≤ |x - xᵢ| < h
        ⎱ 0                 , otherwise

derivative of basis function:
d
—— Nᵢ(x) = ⎰ -sign(x - xᵢ) / h  , 0 ≤ |x - xᵢ| < h
dx         ⎱ 0                  , otherwise

FEM basis function
  linear/quadratic ---> rbf
  signal interpolation (sinc) ---> rbf

level set（contour line）（isoline）

「等高線」。函數曲線一樣高的地方。習慣畫成俯瞰平面圖。

連續可微函數的等高線，可以是一條封閉曲線、一條開放曲線、多條上述曲線。

甚至等高線粗細不一。

函數等於定值，形成方程式，即是等高線。

f(x,y)是一群等高線。f(x,y) = 0是一個等高線。f(x,y) ≥ 0是一群等高線，形成區域，區分內部外部。

等高線，無限微小、無限綿密、無限多，填滿整個空間。

一、函數：等高線彼此不相交。

二、連續：等高線沒有端點。等高線升高降低，宛如膨脹消泄。

三、可微：等高線沒有尖角。等高線每一點只有唯一一條切線。

極值位於兩個等高線的交集。注意到極值可能不只一個。

一、函數：兩個等高線，若等高、又最高，則交集是極值。

二、連續：兩個等高線，只會內接或外接，不會交叉。

三、可微：兩個等高線，交會之處，每一點只有唯一一條切線。

gradient

梯度線：等高線的垂直方向。等高線（圓形）與梯度線（水系）相交於鞍點。

稜線、谷線：若內高外低，等高線們一齊凸起之處是稜線、一齊凹陷之處是谷線。若內低外高，則相反。【查無定義】

稜線、谷線：各種尺寸的同心圓（甚至其他形狀）的極值，連成一線。【查無名稱】

各式各樣的等高線

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

容易計算的函數類型

monotone function：單調函數。遞增、遞減。

convex function：凸函數。斜率是單調函數，遞增是凸、遞減是凹。

unimodal function：單峰函數。先遞增（減）再遞減（增）。

quasiconvex function：擬凸函數。等高線是凸的，遞增是擬凸、遞減是擬凹。

https://mjo.osborne.economics.utoronto.ca/index.php/tutorial/index/1/qcc/t

f(a,x) = ax  在第一三象限是擬凹函數、第二四象限是擬凸函數
https://www.desmos.com/calculator/vkru83pare

[A,X] = meshgrid(-10:0.5:10,-10:10);
Z = A .* X;
surf(A,X,Z)

f(a,x) = (ax)^2 四個象限都是擬凹函數，合起來看像個碗（假的凸函數），截面是拋物線
https://www.desmos.com/calculator/roky6zr9lf

[A,X] = meshgrid(-10:0.5:10,-10:10);
Z = A .* X;
Z = Z .* Z;
surf(A,X,Z)

冰淇淋。螺旋上升函數。

level set

多變量函數是箭矢圖。據我所知似乎沒有等高線的概念。

容易計算的多變量函數類型

遞增函數、凸函數、平緩函數都必須重新定義。

我找不到相關資料，事情相當棘手。以下的定義通通尚待確認。

increasing function ⮕ positive semidefinite matrix

多變量函數的遞增函數，定義成一次微分是半正定矩陣。

單變量函數：斜率大於零（分子分母同號）
                           f(b) - f(a)
function f is monotone iff ——————————— ≥ 0
                              b - a

多變量函數：點積大於零（分子分母同號）（轉角介於±90˚）
multivariate function f is monotone iff (b - a) ∙ (f(b) - f(a)) ≥ 0

https://encyclopediaofmath.org/wiki/Monotone_operator

遞增數列，就是數字越來越大。遞增陣列，就是數字往右往下越大、往左往上越小。

遞增函數，就是數字越來越大。遞增二元函數，就是數字往右上的各種方向都越來越大。

輸入輸出推廣成向量，事情變得複雜：

一、採用輸出向量的長度來比大小。然而向量長度沒有負數，沒辦法越來越小，不存在嚴格遞減。這個方式不行。

二、改用輸入向量與輸出向量的點積來比大小。然而大小的意義整個扭曲了，轉角介於±90˚叫做變大，否則叫做變小。

單變量函數是遞增函數：處處斜率都是正數或零。
多變量函數是遞增函數：處處梯度都是半正定矩陣（特徵值都是正數或零）。
線性函數Ax是遞增函數：處處梯度都是Aᵀ，Aᵀ是半正定矩陣。

increasing function ⮕ nonnegative divergence

多變量函數的遞增函數，定義成散度大於等於零。

單變量函數：斜率大於零
function f is monotone iff d/dx f(x) ≥ 0

多變量函數：散度大於零
multivariate function f is monotone iff [ ∂/∂x ] ∙ [ f(x) ] ≥ 0
                                        [ ∂/∂y ]   [ f(y) ]

一維箭矢圖有個特性：遞增函數（處處斜率大於零），箭矢自零散開。遞減函數（處處斜率小於零），箭矢聚集至零。

高維箭矢圖則是：處處散度大於零，箭矢自零散開。處處散度小於零，箭矢聚集至零。

散開的地方稱作源，聚集的地方稱作匯。當函數起起伏伏、而且很多地方是零，就有許多源匯。正負無限大的地方也是源匯。

儘管源匯很吸睛，然並卵。

Lipschitz function ⮕ Johnson–Lindenstrauss lemma

(1-δ) ‖x-y‖² ≤ ‖f(x)-f(y)‖² ≤ (1+δ) ‖x-y‖²

半正定矩陣：向量經過變換，角度差異在九十度以內。

JL矩陣：向量經過變換，長度平方差異在δ倍以內。

Lipschitz function ⮕ diagonally dominant matrix

多變量函數的平緩函數定義成對角優勢矩陣。

monotone: 0 ≤ (f(b) - f(a)) / (b - a)
Lipschitz: -1 ≤ (f(b) - f(a)) / (b - a) ≤ +1
positive Lipschitz => monotone

|f(a) + f(b)| ≤ |a + b|   其中一個變數變號

diagonally dominant matrix = 每個橫條都是Lipschitz function
(對於橫條的非主對角線元素來說。)
(例如Gauss–Seidel，更新x，方程式5x + y = 5是陡峭的，很近就撞到，所以收斂。
(更新式x = 1 - 0.2y，對y而言，斜率-1 ≤ -0.2 ≤ 1是平緩的。)
positive diagonally dominant => positive definite => monotone
https://books.google.com.tw/books?id=KVfSBwAAQBAJ&pg=PA105

convex function ⮕ submodular function

《Learning Submodular Functions》
Submodularity is in many ways similar to concavity of functions
defined on Rn. For example, concavity of differentiable functions
is equivalent to gradient monotonicity, and submodularity is equivalent
to monotonicity of the marginal values:

• Non-negative if f(S) ≥ 0 for all S.
• Monotone (or non-decreasing) if f(S) ≤ f(T) for all S ⊆ T.
• Submodular if f(A) + f(B) ≥ f(A ∪ B) + f(A ∩ B) ∀A,B ⊆ [n].
  or equivalently f(A ∪ {i}) − f(A) ≥ f(B ∪ {i}) − f(B)
  for all A ⊆ B ⊆ [n] and i ∉ B.
• Subadditive if f(S) + f(T) ≥ f(S ∪ T) for all S,T ⊆ [n].
• L-Lipschitz if |f(S + x) − f(S)| ≤ L for all S ⊆ [n] and x ∈ [n].

isocline

函數圖形當中，當輸入輸出是相同集合，則輸入輸出可以重疊。

不重疊繪圖、重疊繪圖，這兩種繪圖方式，目前尚無正式名稱。

兩種繪圖方式當中，函數運算擁有不同的幾何意義。比方來說，函數加法是上升和位移，函數微分是斜率和密度。教科書幾乎只談不重疊繪圖，然而現實世界幾乎都是重疊繪圖。

isocline

time series graph   函數曲線。以時間為主軸。
gradient field plot 移動速度。(dx/dt, dy/dt)向量場。
phase portrait      移動路線。初始條件(x₀,y₀)沿著向量場跑。

nodal set：特徵函數，高度為零的等高線。
nullcline：梯度場，向量為零的等高線。
streamline：梯度場的向量的切線方向連線；勢能場的等高線的垂直方向。

isocline：導數相等之處。
nullcline：導數為零之處。
steady state：導數皆為零之處。走不動而停下來了。
stable：停在穩態、不斷循環。
unstable：走向無限。

電力線（流線）是另一種製圖方式，電流大小是電力線多寡，電流方向是電力線走向。