Function (ℝⁿ ⇨ ℝ)

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)   輸入變數剛好兩個

加減乘除

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

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

f + g

微分

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

請讀者參考本站文件「Sequence」提及的離散版本。此處介紹的是連續版本,只多了個dx:一個無限微小、略大於零的數值。

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

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

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

積分

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

∫ f dx

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

函數積分最簡單的演算法是Rectangle Rule:按照定義來,將面積切割成數條寬度為dx的矩形。

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

對計算學家來說,定積分就是區間和啦。前綴和改成區間和,就這樣而已。

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

當函數是曲線,仍有縫隙。只好改用Parabolic Rule又稱Simpson's Rule:梯形的斜邊改成拋物線。

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

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

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

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

最後介紹MCMC Integration:使用「Rejection Sampling」,以函數曲線下方的取樣次數比例,推定面積。取樣時,自訂x的分布,可達到自訂dx的效果。

UVa 1356 1280 12528 ICPC 3001

Function (ℝⁿ ⇨ ℝᵐ)

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⃗

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

多個函數顯然可以融合成一個多變量函數,然而一個多變量函數一定可以拆解成多個函數嗎?我也不知道。我從未見過有人討論。

multivariable:多變數。輸入變數有許多個。
multivariate:多變量。輸入輸出推廣成向量。簡單想成是輸出變數有許多個。

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/

各種運算如圖所示:

散度

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

                ∂Fx
div(F) = ∇·F = ―――               1D
                ∂x 

                ∂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,加總兩個差異。以朝外為正值。

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

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

旋度

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

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三個平面上的旋轉。

實際應用,例如渦度。

散度定理、旋度定理

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

散度定理:邊界內的向量的散度的總和、邊界上的向量的法線分量的總和(通量),兩者相等。

旋度定理:邊界內的向量的旋度的總和、邊界上的向量的切線分量的總和(環量),兩者相等。

原理是(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」。

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

梯度的散度【目前稱作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。

梯度的旋度

一定是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   |     
             場  ―――> 散度   

Exterior Algebra與Differential Form

https://www.cs.cmu.edu/~kmcrane/Projects/DGPDEC/

數學家建立了「外代數」這一套世界觀,闡釋點積、叉積、梯度、散度、旋度之間的關聯。讀者若有興趣再來學吧。

Function資料結構

數值、符號

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

                  |            1    
 1 ÷ 3 = 0.333333 |   1 ÷ 3 = ———   
                  |            3    
      numeral     |     symbol      
    無法精準記錄   |   橫槓代表分數 
  f(0.01) = 0.01  | f(x) = { x   if x > 0 
  f(0.02) = 0.02  |        { 0   otherwise
      :       :   |                       
      numeral     |     symbol            
    無法盡數記錄   |   以文字描述           

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

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

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

符號轉數值,函數必須取樣、擇鄰。

符號轉數值,計算學家稱作「離散化Discretization」,數學家稱作「有限元素法Finite Element Method」。

Regular Discretization

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

數值轉符號。傾向使用Monotone Cubic Interpolation。

Chebyshev Discretization

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

數值轉符號。傾向使用Polynomial Interpolation。

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

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

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

Random Discretization(Meshfree Discretization)

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

數值轉符號。傾向使用Radial Basis Function Interpolation。

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

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

數值轉符號。傾向使用Weighted Average Coordinates。

Function資料結構: Regular Discretization

加減乘除

對應位置加減乘除。

外插

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

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

微分、積分

數列微分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

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

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

反梯散(一)【目前稱作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] ]

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

反梯度(一)

簡易的做法:切成一條一條,各自積分。

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

反梯度(二)

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

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足夠小就行了。

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

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

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

Function資料結構: Triangular Discretization(Under Construction!)

梯度的散度

事情變得很複雜。一種運算有多種定義方式。

Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
各種定義

Discrete Laplace operators: No free lunch
如果在二維空間 三個願望一次滿足  四個願望就不行了

曲面之梯度的散度【目前稱作Laplace-Beltrami Operator】

多變量函數用在曲面:函數輸入是曲面的參數座標。

https://graphics.stanford.edu/courses/cs468-13-spring/schedule.html
http://webcourse.cs.technion.ac.il/236629/Winter2015-2016/ho_Lectures.html
http://www.cs.cmu.edu/~kmcrane/
http://fernandodegoes.org/
https://graphics.tudelft.nl/Publications-new/2016/VCDPBHB16/
http://www.geometrysummit.org/summerschool/presentations.html
http://www.geometrie.tugraz.at/sgp2015/gradschool.php

Mathematical Morphology

Mathematical Morphology

Histogram3D[{{0,0},{0,0},{0,0}},{1},Boxed->False,Axes->None]

【註:我製圖技術差,只好介紹陣列版本,而非函數版本。】

調整函數形狀的學問。因為不是顯學,所以名稱矯揉造作,不像一般的數學名詞那樣地簡潔有力。關鍵字grayscale morphology。

基本操作是dilation和erosion,進階操作由基本操作組成。

            dilation:鄰近格子取最大值。功效是補厚。
             erosion:鄰近格子取最小值。功效是削薄。
             closing:先 dilation 再 erosion。
             opening:先 erosion 再 dilation。
   top-hat transform:原函數減掉 opening。
bottom-hat transform:closing 減掉原函數。

鄰近格子可以自訂樣式。另外,大樣式多半可以改為小樣式,效果相同而且時間複雜度更低。

例如5×5,改為兩波3×3,亦得取得5×5範圍內的最小值(最大值)。本來讀取X×Y×5×5格,變成只讀取X×Y×3×3×2格。

運算不可逆、不可抵銷,沒有反運算。個人推測這些運算可以減少亂度。

額外使用過濾函數

進階版本。額外設計一個過濾函數。設計不同的過濾函數,得到不同的效果。

窮舉原函數的每個位置;針對一個位置,令過濾函數的中心對準該位置,各個位置點對點相加(相減)後,才取最大值(最小值)。

邏輯運算的版本

不知道為什麼,影像處理教科書非常喜歡討論邏輯運算的版本。數值改成false和true,最大值(最小值)改成OR(AND),功效是增厚(削薄)圖形的外緣。

UVa 12702