Function (ℝⁿ ⇨ ℝᵐ)(Under Construction!)

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

(處處有數據無法作圖,因此圖例只畫其中幾處。)

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

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

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

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

(圖例為二維純量場,省略座標軸。)

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 Operation

除了加減乘除模微積,又多了幾種運算。數學家並未命名運算名稱,只命名運算結果。

    	      result (noun)
----	----  --------------------------------------------
∇·  	div   divergence  散度
∇×  	curl  curl        旋度
∇   	grad  gradient    梯度
∇·∇ 	---   Laplacian   梯度的散度(運算符號從右往左疊)
∇×∇ 	---   ---         梯度的旋度(總是0,缺乏討論意義)
註:∇·∇經常簡寫成∇²或∆。
註: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,加總兩個差異。以朝外為正值。

物理意義是每一處的膨脹收縮多寡,strain。

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

旋度

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

用抽象的、簡潔的數學符號表達函數旋度:

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,加總兩個差異。以逆時針為正值。

物理意義是每一處的自旋多寡,vorticity。

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

實際應用,例如渦旋。

散度定理、旋度定理

任意向量場、任意封閉範圍,觀察邊界內、邊界上的向量:

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

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

原理是(b-a)+(c-b)+(d-c)+...+(z-y) = z-a,微分之後的區間和,等於頭尾兩項相減。

flux circulation

通量:邊界上的向量的法線分量的總和。朝外為正值。

環量:邊界上的向量的切線分量的總和。逆時針為正值。

梯度

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

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

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

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

X方向的左右差異,以朝右為正值。Y方向的上下差異,以朝上為正值。

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

梯度的散度

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

     d²Fx
∆F = ――――                 1D
     dx² 

     ∂²Fx   ∂²Fy
∆F = ―――― + ――――          2D
     ∂x²    ∂y² 

     ∂²Fx   ∂²Fy   ∂²Fz
∆F = ―――― + ―――― + ――――   3D
     ∂x²    ∂y²    ∂z²

物理意義有兩種解釋方式:

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

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

梯度的旋度

一定是0。不存在樓梯幻覺(Penrose stairs)。

靜態平衡、動態平衡

靜態平衡:相鄰高低差距,處處等於零,不生不滅。高者抑之,下者舉之。數學式子是∇f = 0。梯度為零,是「常數函數」。

動態平衡:相鄰出入差距,處處等於零,因果輪迴。萬物並作,吾以觀復。數學式子是∆f = 0。梯度的散度為零,梯度的旋度總是零,是「和諧函數」。和諧函數的梯度就是諧場。

慣性運動:增加了時間維度之後的靜態平衡。

簡諧運動:增加了時間維度之後的動態平衡。

源、匯

一維箭矢圖有個特性:遞增函數(斜率處處大於零),箭矢自零散開。遞減函數(斜率處處小於零),箭矢聚集至零。

多維箭矢圖則是:散度處處大於零,箭矢自零散開。散部處處小於零,箭矢聚集至零。

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

儘管源匯很吸睛,不過沒什麼用途就是了。

Helmholtz Decomposition(Under Construction!)

公式

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

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,  div(C) = 0,  div(H) = 0
          curl(D) = 0, curl(C) ≠ 0, curl(H) = 0

散旋諧分解至今沒有其他演算法,只有一道數學公式:

{ D = grad(laplace⁻¹(div(F)))
{ C = curl(laplace⁻¹(-curl(F))) = -curl(laplace⁻¹(curl(F)))
{ H = F - D - C

多變數函數之散場、旋場、諧場,線性變換之縮放、旋轉、位移,兩邊有點像。不過諧場不一定是常數。

Helmholtz Decomposition

F = grad(d) + grad(c) = D + C

將諧場隨意分攤至散場、旋場,稱之為無旋場irrotational field、無散場solenoidal field。分解方式無限多種。

Clebsch Decomposition

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

向量場拆成兩個梯度場的線性組合F=aX+bY,令其中一個係數是常數純量場a=1。分解方式無限多種。

整個式子求旋度,可以發現原向量場的旋度,垂直於B與Y。

curl(F) = curl(B) × curl(Y)

找到適當的B與Y,旋度可以畫成圖形。

https://www.youtube.com/watch?v=HUFh8cEDeII
http://www.its.caltech.edu/~achern/projects/Clebsch/
https://en.wikipedia.org/wiki/Clebsch_representation
https://www.physicsforums.com/threads/helicity-and-clebsch-decomposition.29828/

延伸閱讀:Exterior Algebra與Differential Form

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

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

Function資料結構: 等距取樣(Under Construction!)

數值、符號

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

例如實數:

                  |            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           

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

採用數值,函數必須取樣、擇鄰,有兩種方式:

一、等距取樣,形成正方形網格、正方體網格。

二、隨機取樣,連成三角形網格、四面體網格。

本篇文章只介紹等距取樣。

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]

符號轉數值

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

數值轉符號

多項式內插。

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

改進方式是Chebyshev Discretization:在半圓上面等距取樣,再投影至直線,使得中央疏、左右密。詳細介紹請參考《Spectral Methods in MATLAB》。函式庫請參考chebfun

加減乘除模

對應位置加減乘除模。

複合

錯綜複雜,查無資料和應用,姑且略過。

微分、積分

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

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

1st-order derivative
backward [f(x) - f(x-h)] / h
 forward [f(x+h) - f(x)] / h
 central [f(x+h) - f(x-h)] / 2h

2nd-order derivative
backward [f(x) - 2f(x-h) + f(x-2h)] / h²
 forward [f(x+2h) - 2f(x+h) + f(x)] / h²
 central [f(x+h) - 2f(x) + f(x-h)] / h²

1st-order integral
backward [... + f(x-h) + f(x)] ⋅ h
 forward [... + f(x-2h) + f(x-h)] ⋅ h
 central [... + f(x-3h) + f(x-h)] ⋅ 2h

公式源自泰勒展開式,捨去後續項次。當取樣間距小於1、接近0,後續項次迅速縮小,微不足道,無關緊要。

Taylor expansion
f(x+h) = f(x) / 0! + f'(x) h / 1! + f"(x) h² / 2! + ...
f(x-h) = f(x) / 0! - f'(x) h / 1! + f"(x) h² / 2! - ...

1st-order central derivative
f(x+h) - f(x-h) = 2 f'(x) h + ... ≈ 2 f'(x) h
f'(x) ≈ [f(x+h) - f(x-h)] / 2h

2nd-order central derivative
f(x+h) + f(x-h) = 2 f(x) + f"(x) h² + ... ≈ 2 f(x) + f"(x) h²
f"(x) ≈ [f(x+h) + f(x-h) - 2 f(x)] / h²

1st-order central integral 
f'(x) ≈ [f(x+h) - f(x-h)] / 2h
f(x) ≈ [∫f(x+h) - ∫f(x-h)] / 2h
f(x-h) ≈ [∫f(x) - ∫f(x-2h)] / 2h
∫f(x) ≈ ∫f(x-2h) + f(x-h) ⋅ 2h
∫f(x) ≈ [... + f(x-3h) + f(x-h)] ⋅ 2h

當超出邊界,則在邊界附近取點做內插,建立延長線,以估計數值。取一點是複製邊界,取兩點是直線延伸,取三點則是拋物線延伸。計算學家稱作「外插Extrapolation」。

最後提醒一下,陣列微分a[x] - a[x-1],函數微分[f(x) - f(x-dx)] / dx,函數離散化微分[f(x+h) - f(x-h)] / 2h。它們不一樣。

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

都是微分來的。

梯度反元素【尚無正式名稱】

梯度反元素不一定存在。沒有旋場,才有梯度反元素。

梯度的散度反元素【目前稱作Poisson Equation】

兩種解法:線性方程組、傅立葉轉換。

梯度的散度寫成線性方程組。

[ 1 -4  1  0 0  0 0 0 ... ] [ f[0][0] ]   [ g[0][0] ]
[ 0  1 -4  1 0  0 0 0 ... ] [ f[0][1] ]   [ g[0][1] ]
[ 0  0  1 -4 1  0 0 0 ... ] [ f[0][2] ] = [ g[0][2] ]
[ 0  0  0  0 1 -4 1 0 ... ] [    :    ]   [    :    ]
[ :  :                    ] [    :    ]   [    :    ]

解線性方程組,中規中矩的演算法是高斯消去法,時間複雜度O(N³)。不過大家習慣使用速度更快的「Preconditioner」,時間複雜度O(N²K)。

頻域解法。視作摺積,甚至可以推導公式解。時間複雜度O(N²logN)。

梯度廣義反元素【尚無正式名稱】

梯度反元素可能不存在。改為找平方誤差最小的梯度場,有兩種想法:

一、梯度,求最小平方誤差解。三種公式解、共軛梯度法。

二、梯度的散度,求精確解。高斯消去法、Preconditioner。

given F, find G = grad(g)

{ Gx(x,y) = g(x,y) - g(x-1,y)
{ Gy(x,y) = g(x,y) - g(x,y-1)

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

  ∂g
―――――― = g(x+1,y) - 2g(x,y) - g(x+1,y) - fx(x) = 0
∂(x,y)

「一次微分等於零」的地方是極值、鞍點。對每一處的變數偏微分,使之為零。最後得到梯度的散度。

左右兩式相減平方最小,就是左右兩式一次微分相等。縱向與橫向加總,就是左右兩式散度相等。梯度相減平方最小,就是左右兩式梯度的散度相等。

依據是左邊版本的梯度,結論卻是中間版本的梯度的散度。儘管不合邏輯,但是計算學家都默許了──總之h足夠小就行了。

梯度廣義反元素,轉化成梯度的散度反元素。問題解決了。

散場、旋場、諧場

散場:向量場求散度。解Poisson,得散勢。求梯度,得散場。

旋場:向量場求旋度。三個值變號,分別解Poisson,得三個旋勢。求旋度,得旋場。

二維旋場較好算:二維向量場求旋度,只有z維有值。一個值變號,解Poisson,得旋勢。求梯度、轉90度,得旋場。

諧場:向量場減散場、減旋場,得諧場。

想要節省儲存空間,可以改為記錄散勢、旋勢(二維)、諧勢,要用時才算梯度。補充一下散勢、旋勢、諧勢均不唯一。

散場加諧場,正是平方誤差最小的梯度場。

縱場、橫場

散場、旋場的傅立葉轉換。

http://physics.stackexchange.com/questions/1115/

散度反元素、旋度反元素【尚無正式名稱】

散度反元素、旋度反元素不唯一。

(∇×)⁻¹B = A + ∇φ
(∇·)⁻¹φ = A + ∇×B
Biot-Savart Operator
curl(A) = B   --->  A = curl⁻¹(B) = BS(B)
given B, solve A, such that curl(A) = B (and div(A) = 0)

div(curl(F)) = 0  --->  div(curl(A)) = div(B) = 0
                  --->  B is curl field + harmonic field

use DCH decomposition to get A and B ---> div(A) = 0

Function資料結構: 隨機取樣(Under Construction!)

隨機取樣

我尚未通盤理解。提供一些資料供各位參考:

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

這些資料將多變量函數用在流形表面:函數輸入是二維流形表面的參數座標。

符號轉數值

取樣:請參考本站文件「Sampling」。

擇鄰:圓半徑之內,都是鄰居。

數值轉符號

多項式內插。

多項式內插有個嚴重問題:點數過多過密,導致臺階頻繁、震盪頻繁。Runge Phenomenon。

改進方式是Radial Basis Function Interpolation:一、設計一個RBF。二、每個位置套用RBF的加權平均值,做為內插函數(想法類似Mixture Model)。三、多項式內插,求得權重。

Functional Equation(Under Construction!)

Functional Equation

普通的方程式,已知數、未知數全是實數。實數運算有加減乘除模……。

x² + 2xy + 2y = 1

「函數方程式」,已知數、未知數全是函數。函數運算又多了代入複合微積散旋,變化更多了。

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

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

Functional

普通的函數,輸入、輸出全是實數。「泛函數」則全是函數。

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

簡單來說,「泛函數」就是函數的函數。

此處的functional是一個特別的名詞。我不清楚數學家為何故意讓「泛函數(名詞)」跟「函數的(形容詞)」撞名。

範例:古典力學

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

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

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

方才的位置,是一維數線上面的位置,是一個數值。位置可以在二維平面、三維空間,而函數輸出就是二維向量、三維向量了。

物理課教過直線運動。位置是一維數線上的位置。位置的變化快慢,稱作速度,符合微分的概念。x′就是速度。

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

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

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

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

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

範例:流體力學

比方說,記錄物體速度。根據人類目前所知,每個維度互不干涉,符合分量的概念。因此物體速度可以表示成一個向量u。

數學家創造了分量,主要是為了符合人類認知。如果維度互相干涉,那麼數學家勢必要創造其他數學元件,以符合新認知。

水中的每個位置,都有一個水分子。每個水分子都有速度,符合場的概念。水分子的速度可以表示成一個三維向量場u。

水分子往周圍對流,那就是u = Δu。

水分子受重力影響,那就是u″ = g。

公式解、數值解

絕大多數的函數,難以採用符號記載。數學家的說法是:沒有封閉式(僅使用事先指定的符號)、沒有解析式(僅使用加減乘除微積乘方開方對數等常見運算,以及πe等常見符號)、沒有代數式(僅使用加減乘除乘方開方)、……。

難以採用符號記載,只好採用數值記載。函數擁有綿密的、無限多的函數點,難以計算與儲存。解法是取樣、擇鄰。選取某些x當作樣本,同時建立樣本之間的鄰居關係。

大家習慣先試符號解(公式解),再試數值解。詳細流程:查閱工程數學教科書,手工推導符號解,寫成程式碼。然而,大多數時候,函數很複雜,甚至函數不是多項式函數,難以手工推導符號解。何況目前也沒有特別好的演算法,能讓電腦自動推導符號解。最後只好運用下面段落的演算法,求得數值解。

Differential Equation(Under Construction!)

Differential Equation

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

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

f + 2 = f' + 2g   省略x的部分

註:函數微分一次兩次三次,經常簡化成f′ f″ f‴或ḟ f̈ f⃛。

函數運算多且雜。數學家先從微分運算下手,因而有了微分方程式。數學家又細分為ODE:微分變數只有一種、PDE:微分變數超過一種。或者細分為線性、非線性。

經典的微分方程式

基本款式。

Laplace Equation    Δf = 0
Poisson Equation    Δf = ∇g
Helmholtz Equation  Δf = λf

引入了時間變數。請觀賞圖片

Advection Equation  ∇f = k df/dt
Heat Equation       Δf = k df/dt
Wave Equation       Δf = k d²f/dt²

進階款式。請參考大全

Burgers' Equation                衝擊波
Korteweg-de Vries Equation       淺水波
Navier-Stokes Equation           流體力學
Hamilton-Jacobi-Bellman Equation 古典力學

UVa 199

時域解法、數值解(Runge-Kutta Method)

http://heath.cs.illinois.edu/scicomp/notes/index.html

非線性方程式難以求得公式解,目前的解法是模擬。每次只改動一個變數,輪流改動。生動的譬喻是攀岩原則:三點不動一點動。

一次走一步,每一步從最高次導數遞推到最低次導數。

Euler / Verlet
finite difference method
finite volume method

頻域解法、公式解(Galerkin Method)

線性微分方程,特徵函數通常是複數螺旋線eit。因此任意函數的微分,可以寫成特徵函數的線性組合。

http://www.thefouriertransform.com/applications/differentialequations.php
http://www.thefouriertransform.com/applications/pde.php

http://www.sd.rub.de/downloads/Galerkin_method

假設正確答案是某一套函數基底的線性組合。把問題簡化成解線性方程組。

頻域解法當中,因為離散傅立葉轉換只能用在循環數列,所以必須追加假設:在邊界外面函數不斷循環。左右邊界必須一致。

Boundary Condition

微分方程的解,通常有很多種。想要得到唯一解,常見的方式是設定「邊界條件」:妥善設定數值範圍邊界的函數值。

真實世界的應用當中,我們通常只關心一小段數值範圍,而不是從負無限大到正無限大。此時更需要邊界條件。

最知名的兩種邊界條件:

1. Dirichlet condition  已知邊界的函數值(純量場)
2. Neumann condition    已知邊界的函數梯度值(梯度向量場)

追加邊界條件之後,更難求解。幸好,當邊界條件是線性函數,一些簡單的微分方程式如Poisson Equation,還是可以找到公式解。不過此處省略這些數學知識。

範例:Vibration of Membrane

Bessel Function。Zernike Polynomial。

Vibration Modes、Dirichlet Eigenvalue。

Legendre polynomial:  laplace solution of sphere coordinate
Bessel function    :                      cylinder

Chebyshev polynomial: sine wave wrapped around a cylinder
The Sounds of Physical Shapes
integrate { (eigenfunction)^2 * dV }
sum the amplitude for all position ---> sound intensity

Differential Equation: Navier-Stokes Equation (Under Construction!)

Navier-Stokes Equation

描述流體的數學式子,諸如氣體、液體都是。

總共兩種力:

1. 擠壓(法線)
2. 黏滯(切線)

diffusion、viscosity、gravity、pressure。

兩種模擬方式:一、模擬每個分子的速度。二、模擬每個位置的分子總數量(總質量)、總速度。

計算學家稱作粒particle、格grid。數學家稱作Lagrange Description、Euler Description。

粒(Particle)(Lagrange Description)

首先建立一百個粒子,隨機散佈。

一、擠壓。粒子太近產生斥力。設定方式:彈簧力、L、Wpoly。

壓力總和具有擴散效果。

不可壓縮不好處理,當粒子開始重疊,斥力瘋狂上升。

兩層迴圈窮舉所有點對,時間複雜度O(N^2),太慢了。大家習慣採用「Uniform Grid」資料結構,檢查周圍九宮格,時間複雜度O(N)。

二、黏滯。概念類似摩擦力。方才的力,切線分量,乘以黏滯係數。

粒子移動,依照力的大小,F=ma。

需要重力的話,就加上一常數9.8。就這樣。

撞牆有兩種解法:一、牆是不動的粒子:牆給予粒子擴散力,但是牆不受力。二、彈性碰撞:位置反轉回來、速度顛倒方向(變號)。

格(Particle)(Euler Description)

http://www.intpowertechcorp.com/GDC03.pdf

一、擴散。梯度的散度。

二、重力。加上一常數9.8。就這樣。

三、黏滯。這真是個好問題。【待補文字】

四、擠壓。減去散場,保留旋場。

最後平流的時候,採用倒回線性內插。

粒與格衍生許多演算法

SPH 粒連續 http://www.cs.ubc.ca/~rbridson/fluidsimulation/GameFluids2007.pdf
SEC 粒離散 http://graphics.pixar.com/library/SEC/
PBD 粒方程 http://blog.mmacklin.com/
PIC 粒化格 http://www.cs.ubc.ca/~rbridson/fluidsimulation/2006/ParticleInCell.ppt
RBM 格八方 http://matthias-mueller-fischer.ch/realtimephysics/slidesNils.pdf
MG  格中格 https://en.wikipedia.org/wiki/Multigrid_method

這些又衍生許多演算法,比方說Affine PIC(源自PIC)。

未涉及的主題

粒子角動量。
層流、亂流。laminar、turbulence(= mean + fluctuation)
多種物質。instablity
表面張力。adhesion cohesion
溫度。temperature
液壓。
煙、火。
七大數學難題。Navier-Stokes existence and smoothness

Navier-Stokes Equation數學符號

數學式子有點長,這裡不剖析。若有興趣,請自己去看:

http://15462.courses.cs.cmu.edu/fall2015/lecture/fluids

Navier-Stokes Equation: Stream Function

streamfunction-vorticity formulation   curl curl F = - laplace F   (2D)
pressure-velocity formulation          div grad f =  laplace F

Navier-Stokes Equation等號兩邊同時求旋度。以旋度為基礎。

優點是渦流明顯。leap-frogging。微下擊暴流Microburst。

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

Schrödinger Equation

http://www.its.caltech.edu/~achern/projects/SchrodingersSmoke/

計算幾何有個實數變複數的手法:二維座標(x,y)改成複數x+yi,座標平移、縮放、旋轉改成複數加法、倍率、乘法,簡化流程。

最近有人將Navier-Stokes Equation的實數運算方式,改成了Schrödinger Equation的複數運算方式。三維實數改成二維複數,質量(密度)改成複數長度平方,速度改成複數角度,簡化流程。

詳細內容我沒弄懂。請讀者自行學習。

Integral Equation(Under Construction!)

Integral Equation

Gauss quadrature
Ewald summation
http://homerreid.dyndns.org/teaching/18.330/Notes/EwaldSummation.pdf
MCMC integration
1/e^(x^2) sqrt(pi)
1/x       ln(x) 發散
1/(x^2)   pi^2 / 6  (離散版本)
1/x!      e         (離散版本)

Curvature Equation(Under Construction!)

curvature = laplacian = stress

https://en.wikipedia.org/wiki/Brachistochrone_curve
https://en.wikipedia.org/wiki/Tautochrone_curve
https://en.wikipedia.org/wiki/Cycloid
https://en.wikipedia.org/wiki/Whewell_equation
https://en.wikipedia.org/wiki/Cesàro_equation
https://en.wikipedia.org/wiki/Catenary
https://en.wikipedia.org/wiki/Tractrix
https://en.wikipedia.org/wiki/Evolute