上一篇 我們了解到了物體軌跡在物理意義上符合某種三角函數關係,以及在尾篇提到的傅立葉級數形式可以表示其他很多複雜函數,傅立葉變換也可以將一個目標函數分解為一系列正弦和餘弦函數的線性組合。今天我們去探索下基於三角函數的傅立葉級數在函數表示的貢獻,以及與經典數學其他分支的聯繫,比較偏理論,可能會引起不適,之後我們用 傅立葉變換 做些有趣的嘗試。
轉動的風扇與極曲線
我想先對上一篇表示平面圓的極坐標 (polar coordinates) 做一個補充,對極坐標 (r , θ) 圓的半徑 r 是固定的,那麼如果 r = f(θ) 會發生什麼呢,如果 r 滿足以下的函數關係,
這個被稱為玫瑰線 (polar rose/rhodonea curve), 看上去像是花瓣,正好可以用來當做風扇的扇葉;這種曲線有個有趣的特性,當k是整數,且k為奇數時則曲線有k個花瓣,當k是偶數時曲線則有2k個花瓣。當然我們可以對 r 開平方根將這個特性強行改成 k 為正整數。以下是具體代碼的實現,
nPetal 是扇葉數,transformations定義了一些列曲線函數,這裡只用了玫瑰線一種,在line 16設置了每隔段時間增加一片扇葉,fanFrame函數主要是繪製風扇的扇心軸和圓形邊框,當然你可以更加調得精細美觀點,核心計算在 fan 這個函數,其實這裡用2π周期就行了,這裡計算了在每個離散角度時候 r 的值,然後通過上一篇介紹的 極坐標系 到 直角坐標系變換 繪製;這裡想簡單講下平面旋轉,
在 css 裡面,我們操作元素旋轉很方便,只要申明 transform: rotate(n deg); 就能完成,但這不足以讓我們認識旋轉的本質,這裡我們對一個點的旋轉也就是 將 (x, y) 到 (x', y') 相對 角度 Θ 進行一次線性映射,
(來自wiki)
x = r × cos( θ ), y = r × sin( θ ) , 在 極坐標的旋轉也就是 (r , a - θ), θ 就是旋轉的角度, 帶入並根據公式展開可以得到線性表達式,
如果將 x, y 記作向量表達 <x, y>. 就可以得到在 歐幾裡得空間下的矩陣表示,
R 被稱做平面旋轉矩陣 (plane rotation matrix)
在之後的 線性代數章節會介紹諸如遊戲裡常用的四元數(Quaternion) 在立體空間的物體旋轉方法
好了,運行canvas代碼,我們的風扇就做好了,
我們再來複習下三角函數
A代表振幅,函數周期是 2π / w, 頻率是周期 w / 2π, b是函數初相位,k在信號處理中稱為直流分量. 我們可以改變這些變量從而改變 風扇的外形,你也可以用 perlin noise 這一類噪音算法將扇葉做得看起來很舊,這部分留給有興趣的人去試試了。
更多類型的極曲線可以參考來自wiki 的下圖去調試 r 的值
(圖片來自wiki)
函數擬合
上一篇文章,我們了解了簡諧運動符合某種正弦/餘弦周期函數關係,也在結尾處提到了我們可以用傅立葉級數對一系列簡單周期函數進行疊加(linear superposition) 近似表示一個 複雜函數 f(x), 且對於f(x)級數表達式的收斂性,滿足閉區間上 狄利克雷條件(Dirichlet's condition) 或者 二次可積函數 (square integrable) ,以下為周期在2π 的傅立葉展開式 (Fourier expansion)[1],
(sine, cosine被稱為基函數 basis)
我們要做的是求出 an, 和 bn 的參數 (coefficients),有個簡單的定理為,如果 f(x) 為偶函數(y軸對稱 f(x) = f(-x) ), 採用傅立葉 傅立葉餘弦展開式;如果為奇函數(原點對稱 f(-x) = -f(x) ) 則採用 正弦展開式,
三角函數的正交性
我們在中學課堂中學習解析幾何中知道,兩條線相交呈直角關係,它們是正交的。在向量空間裡,若兩個n維向量的點積(Dot Product) 為零,我們則稱它們是正交的 (Orthogonal) 或線性獨立(linear independent);這種向量也叫 基向量 (basis vector),
對於函數正交意味著什麼呢,我們說 函數 f ,g正交若它們在閉區間的積分等於零,這種函數也叫 基函數,
對於線性空間,(定理)基向量可以唯一線性表示空間中任意一個向量。我們通過 gram schmidt 過程將這個子空間的這組基向量 通過投影轉換為 該空間的 正交基 (orthogonal basis).
比如三維歐幾裡得空間存在正交基<1,0,0>, <0, 1, 0>, <0, 0, 1>,那麼對於函數空間,每個連續函數可以表示為基函數的線性組合,也就是我們為什麼可以通過線性組合的正弦,餘弦函數(傅立葉基底 Fourier Basis)在周期(ex. [-π, π])內表示一個連續函數。
通過 gram schmidt 過程我們也可以通過單項函數 1, x, x^2 .... 找到內積空間的正交多項式(orthogonal polynomial)w ,
正交基中的任意兩個基的內積 (inner product) 為零,以 傅立葉 為例,
當 w1 不等於 w2, 就是一個在複平面繞原點旋轉的圓,對稱性,內積為 0, 說明 <1, f, g> 屬於正交基,<f, f> = 1, <f, g> = 0, 且此空間完備,滿足一定條件下,在此空間內的任意函數可以由 f 或 g 疊加而成。
傅立葉級數 是 一堆正交基的和 [4],看下圖
比如在物理上有名的 hermite/legendre polynomials 它們在研究. 在不同學科都對這種正交關係比較感興趣,比如統計分布中的伯努利過程(拋硬幣)。
空間與可測性 (messure)
距離怎麼定義? 在數學上 指的是一個集合(set),任意元素間的距離可定義。
定義,設X為一個集合,存在一個實數空間映射關係 d : X x X -> R, 以下 x,y, z 屬於 X, 有以下 非負性,對稱性,滿足三角不等式 的性質 [8],
滿足條件 1,3,4 的 (X, d) 為一個 度量空間 (metric space),d 為集合X的 度量 (metric)。我們中學解析幾何就是基於 2, 3維的有限維度量空間,又稱 歐幾裡得空間(Euclidean Space)。這裡的距離不僅指點與點直線距離,也包含向量(vector),函數等。
在歐幾裡得空間裡,兩個點 p, q, 在n維空間裡的距離記作,(畢達哥拉斯 或 勾股定理)
向量空間(vector space) 或線性空間 屬於 度量空間的子空間,增加了條件2 零向量 和 標量(scalar) 改變其長度不改變向量方向, 向量空間在自然科學 和工程上有廣泛的應用,
兩個向量 p, q 的n緯空間距離
上面我們熟知的距離計算空間又被稱為L2 空間,
柯西序列
如果度量空間中所有柯西序列(cauchy sequence) 都會收斂為一個點 x 且此點屬於該空間X,稱完備空間 (compact space). 柯西序列就是一個收斂的序列 {x1, x2, ....xn}; 下圖更直觀些,
(柯西數列來自wiki)
我們說任何空間X是完備的(complete), 若在X內的一點 a 滿足上式。(Z, Q, R, or C)
如果給向量空間加上 長度(norm),就被稱為 範數向量空間 (Normed Vector Space), 零向量長度為0, aV 為原向量 V 的 |a| 倍,NVS 中向量的 範數 (norm). ||v|| 表示向量v 的norm 或長度。
度量 metric d ?
在 範數向量空間 中, 有
曼哈頓度量 或 計程車距離 (manhattan metric),百度地圖用的距離, p = 1
歐幾裡得度量 (Euclidean metric), 直角距離, p = 2
切比雪夫度量 (Chebyshev metric) ,棋盤距離, p -> ∞
明可夫斯基度量 (Minkowski metric) 廣義化 generalization 上面三種。
M metric D 定義如下,
對於空間的劃分,wiki上有個圖,如下,每種空間都有不同的運算規則 (交換律,比如矩陣相乘就不具備交換 或 非阿貝爾,四則運算,研究這些運算規則,向量空間 公理系統,拓撲等的在數學上屬於 抽象代數的範疇 Abstract Algebra).
內積空間 加上柯西完備 在數學上叫 希爾伯特空間 (Hilbert space),是有限維 歐幾裡得空間的推廣,使之不局限於實數和有限維,又不失完備性(compactness) ,(完備的範數向量空間 叫 巴拿赫空間 Banach Space,任何一個 希爾伯特空間是巴拿赫空間,反之不成立)
比如在 [0, 2𝜋] f 和 f 的積分 存在此空間,在 希爾伯特空間內
這就是為什麼我們能對函數做傅立葉分析,將已知函數展開成傅立葉級數的運算屬於 調和分析 (Harmonic Analysis),其理論應用中 還有著名的小波 (wavelet ), 是一種類 局部化的傅立葉變換。
上面提到 sine,cosine 函數在 [-π, π] 平方可積 (square integrable),在
之後會有一章線性代數(linear algebra) 專門去探討空間中 向量間的關係, 我們會把 css transform 的絕大多數函數用 基於線性代數實現一遍,這樣我們會更加清楚瀏覽器在處理這些規則時候經歷了什麼)
傅立葉級數 例子,
如果我們有這樣一個函數,, f(x) 是偶函數,那麼找出他的傅立葉cosine表達式,
那麼最終我們得到的 f(x) 表達式為,
為了驗證我們寫一段代碼去將這個近似函數(approximated function) 可視化,看看隨著級數增加函數形狀有怎樣的變化,下面是簡單代碼實現,
以上代碼實現了 我們推算出來的傅立葉餘弦展開式 (fourierCosineFunctor) 與 目標函數 (func) 的比較,傅立葉級數初始項數為1,每500毫秒增加1,
通過對級數項n的逐步遞增, 我們可以上圖看到逼近函數擬合的過程,我們可以看見在 n = 10 時候函數幾乎收斂到目標函數。更加吹毛求疵,我們可以計算出級數隨n的均方差(mean squared error), 這是統計上常用的測距方法,
在控制臺打出來,
忘記關這個模擬,去上了個廁所回來發現,n 在354時候,MSE = 0.00009810487466497763,真是肉眼可見的收斂。
(補充) 對於不連續,比如我們上一篇結尾提到的傅立葉級數收斂方波(square wave), 具體看下面的可視化小節。傅
(來自wiki)
我們可以看見對於不連續點的周期函數逼近,點附近出現衰減震蕩現象,且隨著級數項增加,峰值無法下降,對這種帶稜角的信號,稱之為吉布斯現象 (gibbs phenomenon) [4]. 在jpeg等圖像壓縮算法中,由於gibbs現象的影響,圖像中的邊緣(edge) 看起來像是糊掉的, 這種也被稱為 振鈴效應 (ringing artifacts).
級數 - 函數逼近的經典方法
在大學微分課中,我們學過 泰勒展開(Taylor Series), 用無限多項式級數和表示一個函數,這些多項式由這個函數在某一點求導得到,關於自變量在零點 (a = 0) 的求導的泰勒級數又被稱為 麥克勞林級數 (Maclaurin series).
冪函數 power series
表示將函數在 c 處展開成冪級數,不是每一個函數都能被展開成冪級數。
對於冪級數的收斂性,最簡單的方法有根測試 (root test),
r < 1, 數列收斂 (converge)
r > 1, 數列發散 (diverge)
r = 1, 則無法 根測試 判斷
對於
|x| < R, R如果存在,R 叫做 收斂半徑 (radius of convergence).
泰勒級數 taylor series
是一種冪級數,且收斂。(其公式在前面已標記), 其可以由微積分基本定理推導 與 分部積分法 得到,
繼續重複分部積分 可以得到泰勒展開的更多高階項 n
* 接著也可通過 數學歸納法(mathematical induction) 再次證明 這個式子是否 在 k + 1 也成立。
常用的泰勒級數,
我們可以在傅立葉級數到變換中,使用泰勒級數推導出歐拉公式,但是歐拉大神是如何想到這個精妙的式子呢。
應該對 計算機無法模擬函數連續或無窮等理想理論,所以必須將其離散化,泰勒級數無疑提供了很好的離散運算結構。
組合數學:母函數 generating function
排列組合是算法編程經常遇到問題,比如這樣一個問題。
「有n階梯子,一個人從最下面爬,這個人可以一次爬一階梯子或者一次爬兩階梯子,問有多少種方法爬上頂」
解
如果 n = 1, f(1) = 1, 只有一種爬法
n = 2, f(2) = 2, {(1, 1) , 2}
n = 4, f(4) = 5, {(1,1,1), (1,1,2), (2,1,1), (1,2,1), (2,2)}
憑經驗猜 題的解為,
因在 最後一次走的時候,必須是只能走 1 或者 2 步,只用走一步 為 F(n - 1) 走 兩步 為 F(n - 2). 所以走到 F(n) 則為 其和。
這個數列又叫 斐波納西數列 (Fibonacci sequence), 在經典算法中 我們用動態規劃 (dynamic programming ) 將遞歸的時間複雜(time complexity) O(e^n) 優化到 O(n); 但是缺在空間複雜度 (space complexity) 增加了O(n),
function fibo(n) { if (n <= 0) return 0; var f = Array(n + 2).fill(0); f[0] = 0; f[1] = 1; for(let i = 2; i <= n; i++) { f[i] = f[i - 1] + f[i - 2]; } return f[n];}有沒有一種方法使時間與空間複雜為 O(1) 呢,通過 母函數 推導可以得到 [7] , 證明過程太長,有興趣可以看參考連結,
至此我們了解了級數在組合學的一些應用, 對母函數 感興趣的朋友可深入學習下, 必定受益匪淺。
計算 ?
阿蘭圖靈奠定了現代計算機理論基礎,馮諾依曼則創建了通用圖靈計算機物理結構 (Von Neumann architecture) 或 普林斯頓結構,參與計算的值是以二進位為單位的比特,硬體由控制器(control unit),存儲(memory),輸入輸出(input/output),運算器組成(ALU: arithmetic logic unit)。下圖來自wiki,
基本運算規則(減乘除)其實是通過不同的邏輯門基於加法和位移完成的。這個網站可以看見這些基本運算規則算法 在 寄存器級別的模擬 【10】,老實說挺複雜的,光一個加法就有 這麼多種算法。計算機不理解數學中的連續性(continuity),無窮(infinity), 無理數(irrational),計算集合必須是有限可數(countable)且離散(discrete)。計算機很擅長算整數因為能以1,0表示,對於浮點數這種逼近數 (approximation), 每一次細微的誤差會向前疊加(forward error), 很多對計算精確性要求高的一般會使用大數(big number) 的函數庫,使用模擬ALU運算,想要理解 big number 這類函數庫,請參考 【10】.
計算機如何計算 sine,cosine
泰勒級數小節提到了 正餘弦函數的泰勒級數展開 或者我們在正交性部分提到的切比雪夫多項式,都屬於使用多項式函數近似逼近三角函數的構造法。計算過程會大量使用浮點運算(flops), 有些設備可能缺乏硬體乘法器,用級數方法會很耗時。在1959年,Volder提出的CORDIC算法僅用到加法和位移就可計算三角函數,支持了在純硬體的實現,有興趣的程式設計師朋友可以在javascript實現以下此算法。
* 逼近論 (approximation theory) 研究用函數,多項式,三角多項式等代替或逼近複雜函數的數學分支,其內容證明逼近函數的存在性,構造法,定量,近似程度,與插值(interpolation)法聯繫密切。
傅立葉級數可視化 - 頻率分量
傅立葉級數 可以用不同的正弦 餘弦波組成在 周期內的不同函數 f(x), 但是我們並不能很好從視覺上理解。
傅立葉級數表示方波,上一篇和在吉布斯現象提到了方波 (square wave), 我們試著用傅立葉級數表示
實現代碼如下,
直接跑,
n = 3,我們可以看見三個圓的旋轉線性疊加 簡單接近於 方波的形狀,繼續增加 n 會有更多的小圓(高頻的震動)逐步逼近目標函數,
傅立葉 n = 30 的方波級數擬合,看到我們的30個旋轉圓組成的 Spirograph 畫出來的周期方波已經很接近目標函數了,但不管n多大,在函數不連續出的 吉布斯現象非常明顯。既然傅立葉級數能表示任意函數,那我們可以用它來表示任意圖形麼?
無線維 - 傅立葉級數到變換
傅立葉變換(Fourier Transform) 通過線性積分變換,從時域 到 頻域的變換 得到複變函數 F(w), 它包含了組成該信號的各個頻率(freq),振幅(ampl),相位(phase);通過傅立葉逆變換,我們可以還原目標信號 f(t).
通俗地講就是我們把蘋果汁,梨汁,西瓜汁混合倒進傅立葉變換這個黑盒子,分離出蘋果, 梨,西瓜汁。或者我們常用的流行音樂,合成了伴奏,歌手演唱等音軌的,我們可以通過傅立葉變換分離這些頻段進行編輯。
對於周期為 (-a/2, a/2) 函數 f(x) 的傅立葉級數形式寫作 [4],
引入 歐拉公式(Euler Formula),
得到傅立葉級數的複數表達式,我們在級數小節 知道了 sine,cosine 在真正計算機計算中是會被展開成 泰勒級數對於有限正整數n的擬合,相比 指數複數形式,複雜度增加了不少;
將 a → ∞. 及傅立葉級數不再只能用來表示周期性函數,拓展可適用於非周期函數 f,
兩邊乘,(∆k/2π)/(∆k/2π) = 1, 得到,
根據積分的定義我們得到 傅立葉級數的積分形式, (integral formation), 「正向 和 逆向 變換」
二維傅立葉變換
之後文章會用傅立葉變換處理二維的圖片信號,有必要了解其雙積分形式,它的形式寫作,
對於二維傅立葉變換,高頻部分分布在頻域圖邊界,中心和一維FT一樣為低頻主要信號頻段。
複數 - 實分析的拓展
複數(complex number) 寫做 z = a + ib, i 是 -1 的平方根,常用於工程數學和物理。[9]
(註:數學虛數記作i,物理則為 j )
Re(a) = a 為實數部分(real variable) ; img(b) = ib 虛數部分 (imaginary variable),每一個複數都對應了複平面上一個向量。z 的模長 |z| 為 a 平方 加 b 平方 的平方根。複數運算同樣滿足結合律(associative),分配律(distributive),交換律(commutative)。
通過之前的 極坐標轉換, 在複平面(complex plane) 中 橫軸 為實屬部分,縱軸為虛數部分。
前面我們提到 歐拉公式 將複平面的點 z 表示為 極坐標 (polar coord), |z| = 1, 正交基模長為1 [9]
向量的表示exp(iφ) = cos φ + i sin(φ) , (更多向量之後會在線性代數中介紹)
它的共軛複數 (complex conjugate) 為 exp(-iφ), 在下面單位圓相對其橫軸對稱,隨φ增加,其向量作逆時針旋轉。
(來自wiki)
φ ∈ [-2π, 2π], 其為單位圓上的旋轉, 就這樣簡單滴介紹了複數以及表達 :(
我們在大學學的微分積分 屬於是實分析(real analysis) ,那麼在對於複變函數 f(z) , z = x + iy, 微分(differentiation) , 積分(integral) 又有和不同呢,在對 f(t) 進行傅立葉變化得到的 F(w) 是複變函數,對 F(w) 進行逆變換,實分析的工具就不能用了,
微分
複數函數 f(z), z 為復變量,可以寫成,{z = x + iy } [9]
他在點z0 且若 f(z0) 存在 的復導數,
全純函數(holomorphic function)定義在複平面的開子集,在內每一點復可微(complex differentiable) 的函數 或 解析函數 (analytic function), 在整個複平面全純的函數為 整函數(entire funciton).
根據柯西-黎曼方程
積分
光滑曲線長度 對於任何有界參數曲線(parametric curve) y,他的長度,
複變函數積分最重要的 柯西積分定理 (Cauchy's Theorem), 關於複平面全純函數的在閉合路徑(連續,不自交能定義長度)路徑積分為零,說明了 任何閉合區域的全純函數內部值取決於區域邊界的值。
如果 單連通區域 (simply connected domain) 有洞(極點 pole),簡單來說就是被積複函數分母為零的因子,所以滿足柯西定理純函數的條件,對於這種情況,
柯西積分公式 (Cauchy's Integral Formula), 滿足 f 是全純函數,所有 ω 都被包含在一個閉合,正向,不自交的光滑曲線內,則存在以下關係,
留數定理 (residue theorem) 如果z1,.... z5 是正向,簡單,閉合光滑曲線γ 的孤立奇點(isolated singularity) {z | 0 < | z - zi | < R}, R > 0,則除了奇點的圍線積分(contour integral) 為,
解析函數的級數展開 - Laurent Series
假設 f 在 A 空間內是全純函數,及 f 可以被表示成 A 中的 勞倫特級數形式,A = { z ∈ C : R1 < | z - z0| < R2}.
γ 是任意中心在 z0 的圓
定理,如果 f(z) 是周期2π解析函數, 則複變函數 f 有傅立葉級數展開。
設 f(z) = g(exp(iz)),
|ω| = 1, ω = exp(iθ), 0 ≤ θ ≤ π, dω = iω dθ, 勞倫特級數參數 cn 可為,
傅立葉級數與勞倫特級數有一定聯繫。
複變函數 簡單介紹至此,回到 傅立葉級數。
我們可以把 cosine, 和 sine 寫成復指數函數形式,推導出以下,
最終我們可以得到他的積分形式,
自變量 x(t) 表示時間 (s),復變量 w 表示頻率 (hz)
X(w) 為傅立葉變換(FT), x(t) 為傅立葉逆變換 (Inverse FT)
傅立葉變換(Fourier Transform), 將一個函數從時域(time domain) 轉換 頻域 (frequency domain),頻域也就是我們在 傅立葉級數可視化 裡的那些一個個轉動的圓,
在頻域 F(w) 的 實數部分 包含 振幅(amplitude) = A ; 頻率(frequency) = w / 2π (上面周期為 2π / w) ; 虛數部分則包含頻率分量的 相位(phase) b 的信息,也就是該頻率波的起始位置;
我們回頭看看前面提到的 吉布斯現象,當我們對 矩形函數(方波) 進行傅立葉變換後,會發現 頻域 為辛格函數(sinc function = sin (x) / x ),不連續或突變的反映在圖像上就是 圖像的邊緣部分,所以經過DCT壓縮(compression), 重建(reconstruction) 或 傅立葉餘弦逆變換 打開圖片發現邊緣部分有 偽影存在, 特別是經過邊緣銳化後 振鈴現象更為明顯,由於計算機計算不存在無限項,計算過程中用有限項近似無限項從而丟失原始信號中的高頻成分 也就是sinc 的兩側長尾部分或者高頻波段。一般可以用一個窗口函數或地通濾波器 (low pass filter) 過濾掉高頻波段,再用IDCT重建信號,具體之後可能會有簡單在web前端實現 photoshop 一些功能 去具體理解 圖像振鈴現象,
另外 我們也可以簡單總結出,若一個函數越光滑(smooth), (n階可導函數,n 越大 越光滑;若函數為無限可導, 稱之為 C∞ 函數,比如 sine,cosine;當然正餘弦函數無限可導,無限可積,也是其相比多項式正交基表示函數的好處) , 則它在頻域的高頻modes就越少,且收斂更快 [2 頁880]。黎曼-勒貝格定理(Riemann-Lebesgue lemma).
隨機漫步與熱方程
概率論上研究了很多分布,比如前面提到伯努利分布,離散的二項式概率分布。用斯特林公式(stirling formula) 我們可以將伯努利分布轉換成著名的高斯分布(Gaussian Distribution) [5]. 我們生活中絕大多數的現象都符合高斯分布,比如白噪(white noise),考試成績分布,所以在工程,統計,社科中 高斯分布有很大的分析價值,比如卡爾曼濾波器(Kalman filter).
隨機漫步的粒子
假設有個粒子在直線上,它可以以50%概率向左或向右移動,在 τ 時間間隔,它每次移動的距離都是 h 。假設在初始時間 t = 0,這個粒子在 x = x0 這個位置, 我們想知道這個粒子在時間 t 它的位置的函數 u(x, t).
那麼在 t + τ 時候可以寫成以下等式,
用之前提到的 泰勒級數,可以得到 u(x, t + τ ) 和 u(x -+ h, t ) 展開式,
帶入後,我們可到我們 粒子運動位置模型 的導數形式,
h 和 τ 都是常量,如果假設它們等於 2D,D是有限正數。只取到泰勒二階導數,我們可以更加簡化這個式子, 𝛿 是 delta 函數,又叫脈衝函數,是個分段函數,我們可以用來表示 粒子的初始位置的函數 (初始條件 IC: initial condition),粒子活動的範圍 (-∞,∞)又叫邊界條件 (BC: boundary condition),
(dirac delta 來自wiki)
然後得到以下式子,
(*)
這個在數學中屬於偏微分方程 (partial differential equation: PDE) 中的 一維熱動力方程 (heat equation) 因其描述隨時間演化的現象,又歸類於 拋物線PDE (parabolic pde) 其餘兩大類為 雙曲形 和 橢圓PDE,(如果方程中僅含有一個自變量的微分方程,我們叫它常微分方程 (ordinary differential equation: ODE))。ODE的解析解(analytical solution or exact solution) 可以用 拉普拉斯轉換成線性代數方程解出來;但現實是,很多微分建模通常包含至少 時間 和 位置的變量,所以只能用 泛函(functional analysis) 分析工具,比如我們可以用 前面推導出的 傅立葉變換 來解我們的粒子運動方程,但是對於真實場景中複雜的邊界條件 (飛行中飛機機身氣流強度熱力圖 (navier stokes equations)),可能不存在解析解,在工程中逼近(approximatiojn) 解更加普遍及實用。
傅立葉變換求解粒子運動分布PDE
在約瑟夫・傅立葉的著作 熱的解析理論 (The Analytical Theory of Heat), 第一次使用了三角函數解熱量怎麼在立體管道中怎麼傳播,經過後人的努力將傅立葉變換理論逐漸完善化。我們這兒直接用看看怎麼使用 傅立葉變換解出這個1維熱力方程,
對PDE (*) 等式兩側分別進行傅立葉變換,
得到 得到了一個簡單的一階常微分方程 :)直接對它不定積分,
表示 U(ω,t) = F(ω)G(ω,t) , 對它進行 傅立葉逆變換找到解 u(x,t) ,
假設函數可積,可得到 F*G 的卷積 (convolution), (**)
根據卷積定理(convolution theorem), "函數卷積的傅立葉變換是函數傅立葉變換的乘積。即一個域中的卷積對應於另一個域中的乘積 ", 快速傅立葉變換利用此定理可以將計算複雜度從O(n^2) 降低到 O(n log n).
我們對 g 進行傅立葉逆變換,可以將 g(ω) 轉換成一階常微分方程(1st Order ODE), 得到 g(x) 關於 t,通過傅立葉變換將微分方程轉化到頻域空間簡化了其複雜性,
g 帶入 (**),最終得到解 u(x, t),
假設我們的初始條件IC 為 德爾塔(dirac delta)函數, 它滿足 不定積分 = 1,
那麼 s = 0, 最終,u(x, t) = h(x, t),
這個不是挺像 正態分布 的概率密度函數麼 (Gaussian PDF), 我們找到了 熱力方程 與 統計間的聯繫,
方差 , 期望值 為 0。
我們編碼看下它怎麼隨時間變化呢,為了簡便,假設 k = 1, 只觀察 2π 周期,
運行可以得到觀察到以下運動,隨著時間增減,粒子的位置分布會漸漸均勻。
至此,展示極簡情況下的熱方程以及用傅立葉變換簡化了求解。對於其它複雜的偏微分方程需要更多數學分析技巧, 在求解這些方程的解析解中,總結了一系列具有特定性質的函數,在物理和工程上被統一歸納為 特殊函數(special functions), 比如概率上或複變函數積分重要的 Γ 伽馬函數(gamma), 波動方程求解常用的 貝塞爾函數 (bessel) 或者前面提到的 辛格函數 (sinc) 等 [2].
下一篇我們玩一玩有趣的傅立葉變換。畢竟在康納爾大學評出的十大計算機算法中[3],計算離散傅立葉變換的 快速傅立葉 (FFT) 鶴立其中,可見一斑 (不然 美顏磨皮相機 和 變聲算法 都要卡成翔),據說每天FFT在全世界會運行幾十億次。 FFT 又稱 Cooley-Tukey 算法(沒錯就是算法導論的作者) 是使用了算法中分而治之(divide & conquer)的思想 將 離散傅立葉級數 轉化為子問題遞歸求解。
通過三角函數和傅立葉級數的探索我們找到了經典數學分支間的某種聯繫,以及前端工程裡它的身影:)
* 關於連續性,緊緻,封閉這一系例概念屬於點集拓撲的概念,由於太抽象,就不展開了。
在時域空間 紛繁複雜 的世界,在 頻域 卻如此單調,像上帝譜好的譜子。
非數學系碼農一名,若有紕漏,請多包含;篇幅有限,只能廣度優先;在以後的篇幅裡,會結合應用深入介紹,謝謝閱讀,希望對您有所啟發。
參考
1. CONVERGENCE OF FOURIER SERIES, http://www.math.utk.edu/~freire/m435f07/m435f07fourierconvergence.pdf
2. Paul Goldbart's, Mathematics for physics. Book.
3. Top10Algorithms, http://pi.math.cornell.edu/~ajt/presentations/TopTenAlgorithms.pdf
4. FOURIER SERIES AND INTEGRALS, https://math.mit.edu/~gs/cse/websections/cse41.pdf
5. derivation of gaussian distribution from binomial,
https://people.bath.ac.uk/pam28/Paul_Milewski,_Professor_of_Mathematics,_University_of_Bath/Past_Teaching_files/stirling.pdf
6. Albert. B, First course in wavelets with fourier analysis.
http://links.uwaterloo.ca/amath391w13docs/boggess_narcowich.pdf
7. Mihai, N, generating function, http://www.math.toronto.edu/mnica/csplash.pdf
8. Dietmar. S, Theo, B. Functioanl Analysis, book
9. B.V. Shabat. Introduction to Complex Analysis - excerpts
10. http://www.ecs.umass.edu/ece/koren/arith/simulator/