拉普拉斯方法

来自testwiki
跳转到导航 跳转到搜索

Template:Roughtranslation数学上,以皮埃尔-西蒙·拉普拉斯命名的拉普拉斯方法是用于得出下列积分形式的近似解的方法:

abeMf(x)dx

其中的 ƒ(x) 是一個二次可微函数M 是一個很大的數,而積分邊界點 ab 則允許為無限大。此外,函數 ƒ(x) 在此積分範圍內的 全域極大值 所在處必須是唯一的並且不在邊界點上。則它的近似解可以寫為:

abeMf(x)dx2πM|f(x0)|eMf(x0) as M.

其中的 x0 為極大值所在處。這方法最早是拉普拉斯在 (1774, pp. 366–367) 所發表。(待考查)

拉普拉斯方法的想法概論

Template:Cot

首先,我們得先知道的是,這裡所謂的近似解是以 相對誤差 在講,而不是指 絕對誤差 。因此,令

s=2πM|f(x0)|

s 很顯然的當 M 很大時為一個微小的數,則上述的積分可以寫為

abeMf(x)dx=seMf(x0)1sabeM(f(x)f(x0))dx=seMf(x0)(ax0)/s(bx0)/seM(f(sy+x0)f(x0))dy

因此,我們的相對誤差很顯然的為

|(ax0)/s(bx0)/seM(f(sy+x0)f(x0))dy1|

現在,讓我們把積分分為 y[Dy,Dy] 的與餘下的部分。

Template:Cob

Template:Cot

我们现在来研究 M(f(x)f(x0))泰勒展開 。這裡展開的原點取 x0 ,並將變數由 x 換為 y 因為我們是在對 y 做積分,再加上 x0 為函數極大值所在處,因此,

M(f(x)f(x0))=Mf(x0)2s2y2+Mf(x0)6s3y3+=πy2+O(1M).

會發現高於二次微分的部分當 M 增大時,會被壓抑掉,而此 指數函數 則會越來越接近 高斯函數 eπy2eπy2dy=1 。如下圖所示:

圖為 eM[f(sy+x0)f(x0)] ,其中M為 1, 2 與 3。紅色線為 eπy2

Template:Cob

Template:Cot

y[Dy,Dy] 的部分由上面論述可知當 M 越大則越接近高斯函數;並且有趣的是,相對應的 x 區間則越來越小,因為  x[sDy,sDy]M 成反比。

Template:Cob

Template:Cot

所以,即便 Dy 取非常大, sDy 最後還是會變得非常小,那麼,我們怎麼保證餘下的部分的積分也會趨近到0?

簡單的想法是,設法找一個函數 m(x) ,使它在餘下的區域內 m(x)f(x) ,並且隨著 M 的增大,eMm(x) 的積分會趨近於0 即可。而指數函數只要指數部分為實數,它就一定大於0,並且當指數變大時,指數函數也會跟著變大。基於這兩個理由, eMf(x) 的積分也會趨近於0。 在此,我們可以選擇此 m(x) 為過 x=sDy 的切線函數,如下圖所示:

m(x) 以斜直線代表,為通過位於 x=sDy+x0 的切線,當 sDy 越小,涵蓋的範圍越大

若為有限積分區間,我們可以看出,只要 M 夠大,這切線的斜率就會趨近於0 ,因此,不管 f(x) 在其餘部分是否連續,它們總會在這切線以下,並且很容易證明(後面會證) eMm(x) 的積分當 M 很大時,就會趨近於0。

然而,若積分範圍是無限大的話,那麼,m(x) 就有可能總是與 f(x) 有相交,若有相交,則表示無法保證 eMf(x) 的積分最後會趨近於0。比方 f(x)=sin(x)x 的情形下,不管 M 為多少,它的 0eMf(x)dx 總是會發散。所以,對於積分範圍為無限大的情形,還得要求 deMf(x)dx 會收斂才行,那麼,只要 d 越來越大,該積分就會趨近於0。而這 d 可選為m(x)f(x) 相交之處即可,即保證所有餘下的部分的積分會趨近於0。

也許您會問說,為何不是要求 def(x)dx 會收斂?舉個例子,若 f(x) 在餘下區域內為 lnx ,則 ef(x)=1x ,積分會發散,然而,當 M=2eMf(x)=1x2 的積分卻是收斂的。所以,有的函數在 M 小的時候會發散,然而,當 M 夠大時卻會是收斂的。

Template:Cob 基於上述四點,就有辦法證明拉普拉斯方法的可靠性。而 Fog(2008) 又將此方法推廣到任意精確。 ***待考查***

此方法的正式表述與證明:

假設 f(x) 是一個在 x0(a,b) 這點滿足 (1) f(x0)=max[a,b]f(x) ,(2)唯一全域最大,(3) x0 附近為二階可微且 f(x0)<0 (4) 當拉普拉斯方法的積分範圍為無限大時,此積分會收斂,

則,

limn+(abenf(x)dx(enf(x0)2πn(f(x0))))=1

Template:Cot 下限

f 要求連續的條件底下,給定任意大於0的ε 就能找到一個 δ>0 使得在 |x0x|<δ 這區間內所有的 f(x)f(x0)ε. 。 由 泰勒展開 我們可以得知,在 x(x0δ,x0+δ) 區間內, f(x)f(x0)+12(f(x0)ε)(xx0)2

這樣,我們就得到本方法的積分下限了:

abenf(x)dxx0δx0+δenf(x)dxenf(x0)x0δx0+δen2(f(x0)ε)(xx0)2dx=enf(x0)1n(f(x0)+ε)δn(f(x0)+ε)δn(f(x0)+ε)e12y2dy

其中最後一個等號來自於變數變換: y=n(f(x0)+ε)(xx0) 。請記得 f(x0)<0 ,所以我們才會對它的負號取開根號。

接著讓我們對上面的不等式兩邊同除以 enf(x0)2πn(f(x0)) 並且對 n取極限,則

limn+(abenf(x)dx(enf(x0)2πn(f(x0))))limn+12πδn(f(x0)+ε)δn(f(x0)+ε)e12y2dyf(x0)f(x0)+ε=f(x0)f(x0)+ε

既然上式是對任意大於0的 ε 都對,所以我們得到:

limn+(abenf(x)dx(enf(x0)2πn(f(x0))))1

請注意,上面的證明也適用於 a=b= 或雙雙跑到正負無限大的情況

上限

證明上限的部分其實和證明下限的部分很像,但是會較麻煩。再一次,我們取ε>0 ,不過,不再是任意而是多了一個要求 ε 得夠小以致於 f(x0)+ε<0 。接著,就如同之前的證明,因著 f 被要求連續,並且根據 泰勒定理 我們會發現總存在一個 δ>0 以至於在 |xx0|<δ 區間裡, f(x)f(x0)+12(f(x0)+ε)(xx0)2 總是可以成立。 最後,在我們的假設裡 (假設 a,b 都是有限值) ,由於 x0 是全域最大所在處,我們總可以找到一個 η>0 使得所有在 |xx0|δ 這區間裡, f(x)f(x0)η 總是成立。

現在,萬事俱備,東風就在下面啦:

abenf(x)dxax0δenf(x)dx+x0δx0+δenf(x)dx+x0+δbenf(x)dx(ba)en(f(x0)η)+x0δx0+δenf(x)dx
(ba)en(f(x0)η)+enf(x0)x0δx0+δen2(f(x0)+ε)(xx0)2dx(ba)en(f(x0)η)+enf(x0)+en2(f(x0)+ε)(xx0)2dx
(ba)en(f(x0)η)+enf(x0)2πn(f(x0)ε)

如果我們對上面的不等式兩邊皆除以 enf(x0)2πn(f(x0)) 並且順便取極限的話,會得到:

limn+(abenf(x)dx(enf(x0)2πn(f(x0))))limn+((ba)eηnn(f(x0))2π+f(x0)f(x0)ε)=f(x0)f(x0)ε

再一次,因為 ε 可以取任意大於0的值,所以我們得到了上限了:

limn+(abenf(x)dx(enf(x0)2πn(f(x0))))1

把上限與下限兩個證明同時考慮,整個證明就完成了。

注意,關於上限的證明,很明顯的當我們把它應用在 ab 為正負無限大時,該上限證明會失敗。那怎麼辦呢?我們需要再多假設一些東西。一個充分但非必要的假設是:當 n=1 時,此積分 abenf(x)dx 為有限值,並且上面所說的 η 是存在的 (注意,當 [a,b] 區間是無限的時候,這假設是必要的) 。整個證明過程就如同先前所顯示的那樣,只不過下列的積分部分要做點改變:

ax0δenf(x)dx+x0+δbenf(x)dx

必須利用上述的假設,而改為:

ax0δenf(x)dx+x0+δbenf(x)dxabef(x)e(n1)(f(x0)η)dx=e(n1)(f(x0)η)abef(x)dx

以取代先前會得到的 (ba)en(f(x0)η) ,這樣的話,當我們除以 enf(x0)2πn(f(x0)) ,就會改得到如下的結果:

e(n1)(f(x0)η)abef(x)dxenf(x0)2πn(f(x0))=e(n1)ηnef(x0)abef(x)dxf(x0)2π

這樣的話,當我們取 n 時,上式的值就會趨近於 0 。而剩下的部分的證明就還是如同原先的證明,不做改變。

再強調一次,這裡我們多加給無限大積分範圍的情形的條件,是充分,但非必要。不過,這樣的條件已經可以適用在許多情形了(但非全部)。 這考慮條件簡單來講就是積分區間得是被良好定義的(即不能是無限大的),並且被積函數在 x0 必須是真的極大 (意即 η>0 必須真的存在) ;如果這積分區間是無限大的話,要求 n=1 時的此拉普拉斯方法所用的積分值要為有限並非必要的,其實只要當 n 大於某數時,此積分值會是有限的即可。

Template:Cob

其他形式

有時拉普拉斯方法也會被寫成其他形式,如:

abh(x)eMg(x)dx2πM|g(x0)|h(x0)eMg(x0) as M

其中 h為正 (好像不必要)。

重要的是,這方法精確度與函數 g(x)h(x) 有關。 [1] ***待考查***

Template:Cot 首先,由於極大值所在設為 x0=0 並不影響計算結果,所以以下的推導皆如此假設。此外,我們想要的是相對誤差 |R|,所以

abh(x)eMg(x)dx=h(0)eMg(0)sa/sb/sh(x)h(0)eM[g(sy)g(0)]dy1+R

其中 s2πM|g(0)| ,所以,若我們令 Ah(sy)h(0)eM[g(sy)g(0)]A0eπy2 ,則

|R|=|a/sb/sAdyA0dy|

所以,只要我們求得上式的上限為何,即可得相對誤差。注意,這裡的推導還可以再優化,不過,由於此處我只想顯示它會收斂到0,因此,不再細推。

由於 |A+B||A|+|B| ,因此

|R|<|DyA0dy|(a1)+|Dyb/sAdy|(b1)+|DyDy(AA0)dy|(c)+|a/sDyAdy|(b2)+|DyA0dy|(a2)

其中 (a1)(a2) 雷同,所以只算 (a1) ,經過 zπy2 的變換後,得到

(a1)=|12ππDy2ezz1/2dz|<eπDy22πDy

也就是說,只要 Dy 取得夠大,它就會趨近於0。

(b1)(b2) 也雷同,所以也只算一個:

(b1)|Dyb/s[h(sy)h(0)]maxeMm(sy)dy|

其中

m(x)1Mlnh(x)h(0)+g(x)g(0)asx[sDy,b]

並且 h(x) 在此範圍內要與 h(0) 同號。 這裡讓我們選過 x=sDy 的切線為 m(x) 吧!即 m(sy)=g(sDy)g(0)+g(sDy)(sysDy) ,如圖所示:

m(x) 以斜直線代表,為通過位於 x=sDy 的切線

由圖可以看出,當 s 或者 Dy 變小時,滿足上列不等式的區間會變大,因此,為了涵蓋整個區間, Dy 有了上限。此外,eαx 的積分也相對容易,因此,很適合用來預估誤差。

由泰勒展開我們得知,

M[g(sDy)g(0)]=M[g(0)2s2Dy2+g(ξ)6s3Dy3]asξ[0,sDy]=πDy2+(2π)3/2g(ξ)Dy36M|g(0)|3/2,

Msg(sDy)=Ms(g(0)sDy+g(ζ)2s2Dy2),asζ[0,sDy]=2πDy+2M(π|g(0)|)3/2g(ζ)Dy2

然後將它們代回 (b1) 的計算裡,不過,您可以看到,上述的兩個餘項皆與 M 的開根號成反比,為了式子的漂亮,讓我省略它們吧!不省略只是看起來較醜陋而已,不過,那樣子較嚴謹便是。

(b1)|[h(sy)h(0)]maxeπDy20b/sDye2πDyydy||[h(sy)h(0)]maxeπDy212πDy|.

所以,原則上也是 Dy 越大則越趨近於0。只不過,在決定 m(x) 的過程,設下了 Dy 的上限。

至於靠近極大值的點的積分,我們一樣可以利用泰勒展開來求,當 h(x) 在此點的一階微分不為0時,

(c)DyDyeπy2|sh(ξ)h(0)y|dy<2πM|g(0)||h(ξ)h(0)y|max(1eπDy2)

會看到它原則上與 M 的開根號成反比,其實,當 h(x) 為常數時,積分出來的也會有如此的行為。

所以,概括的講,靠近極大點的附近的積分原則上會隨著 M 的增大而變小,而其餘的部分,只要 Dy 夠大,也會趨近於0,只是這Dy 是有上限的,取決於我們所找的上限函數 m(x) 是否在這區間內總是大於 g(x)g(0) ;不過,一旦有一個滿足條件的 m(x) 被找到,由於切線的起點為 x=sDy ,因此,Dy 可以取正比於 M 即可,所以, M 越大,Dy 也可以設越大。 Template:Cob

至於多維的情形,讓我們令 𝐱 是一個 n 維向量,而 f(𝐱) 則是一個純量函數,則此拉普拉斯方法可以寫成

eMf(𝐱)dnx(2πM)n/2|H(f)(𝐱0)|1/2eMf(𝐱0) as M

其中的 H(f)(𝐱0) 是一個在𝐱0 取值的 海森矩陣,而 || 則是指矩陣的 行列式 ;此外,與單變數的拉普拉斯方法類似,這裡的 海森矩陣 必須為 負定矩陣,即該矩陣的所有本徵值皆小於0,這樣才會是極大值所在。.[2]

拉普拉斯方法的推廣:最速下降法

此拉普拉斯方法可以被推廣到 複分析 裡頭使用,搭配 留數定理 ,以找出一個過最速下降點的 contour (翻路徑的話,會與path integral 相衝,所以,這裡還是以英文原字稱呼) 的 曲線積分 ,用來取代原有的複數空間的 contour積分。因為有時 f(z) 的一階微分為0的點未必就在實數軸上,而就算真在實數軸上,也未必二階微分在 x 方向上為小於0 的實數;此種情況下,就得使用最速下降法了。由於最速下降法中,已經利用另一條通過最速下降的鞍點來取代原有的 contour 積分,經過變數變換後就會變得有如拉普拉斯方法,因此,我們可以透過這新的 contour ,找到原本的積分的漸進近似解,而這將大大的簡化整個計算。就好像原本的路徑像是在蜿蜒的山路開車,而新的路徑就像乾脆繞過這座山開,反正目的只是到達目的地而已,留數定理已經幫我們把中間的差都算好了。請讀 Erdelyi (2012)與 Arfken & Weber (2005) 的書裡有關 steepest descents 的章節。

以下就是該方法在z 平面下的形式:

abeMf(z)dz2πMf(z0)eMf(z0) as M.

其中 z0 就是新的路徑通過的鞍點。 注意,開根號裡的負號是用來指定最速下降的方向,千萬別認為取 f(z0) 的絕對值來取代這個負號,若然,那就大錯特錯了。 另外要注意的是,如果該被積函數是 -{zh-tw:亞純函數;}- ,就有必要將被新舊 contour 包到的極點所貢獻的留數給加入,範例請參考 Okounkov 的文章 arXiv:math/0309074 的第三章。

更進一步一般化

最速下降法還可以更進一步的推廣到所謂的 非線性穩定相位/最速下降法 (nonlinear stationary phase/steepest descent method)。 這方法主要用在解非線性偏微分方程,透過將非線性偏微分方程轉換為求解柯西變換(Cauchy transform)的積分形式,就可以藉由最速下降法的想法來得到非線性解的漸進解。

艾里方程(線性) y=xy 為例,它可以寫成積分形式如下:

yj(x)=12πiγje13t3xtdt

由這條積分式,我們就可以藉由最速下降法(若 γj 指的是負實數軸,那麼就回到此拉普拉斯方法了)來得到它的漸進解了。

然而,若方程式如 KdV方程 yt+6yyx+yxxx=0 是個非線性偏微分方程,想要找到它的解相對應的一個複數 contour 積分的話,就沒那麼簡單,在非線性穩定相位/最速下降法裡所用到的概念主要是基於散射逆轉換 (inverse scattering transform) 的處理方式,即先將原本的非線性偏微分方程變成 Lax 對 ,其中一個像是線性的 薛丁格方程式 ,其位能障為我們要找的 y ,本徵值為 λ ,波函數為 ϕλ (不過,它並非我們所要的 y );因此可以解它的散射矩陣,若利用解析延拓將原本的波函數由實數 λ 延拓到複數空間時,就可以得到黎曼希爾伯特問題(RHP)的形式。利用這個黎曼希爾伯特問題(RHP) ,我們可以解得 ϕ 的柯西變換的積分形式,再利用此線性薛丁格方程的特性,就可以反推出 y 的複數 contour 積分 形式了。

Lax 對 的另一個偏微分方程則是決定每個 ϕλ 隨時間變化的行為,由於 yx± 時被要求為0 ,會發現整個偏微分方程會變得十分簡單,並且只決定 c(λ,t)ϕλ 裡的 c(λ,t) 的值,不過,條件是 ϕλ 必須是指由正負無限遠入射或反射波的解。這樣,我們所得到的這個只與時間與本徵值有關的係數 c(λ,t) 就可以直接被應用在上述的黎曼希爾伯特問題(RHP)裡的躍變矩陣裡了。

接著就是非線性穩定相位/最速下降法所要做的工作,即找出 鞍點 來,在該點附近基於最速下降法的精神做近似。不過,這近似因著考慮到收斂性,需要將原本的 contour 變形,與將原本的黎曼希爾伯特問題(RHP)作轉換,所以有再比原本的最速下降法多出一些步驟來。

這整個方法最早由 Deift 與 Zhou 在 1993 基於 Its 之前的工作所提出的,後續又有許多人加以改進,主要的應用則有 孤波 理論,可積模型等。

複數積分

以下的積分常用在 拉普拉斯變換#拉普拉斯逆變換 裡,

12πicic+ig(s)estds

假定我們想要得到該積分在 t>>1 時的結果(若 t 為時間,通常就是在找經過夠久時間後達穩定的結果),我們可以透過 解析延拓 的概念,先將這時間換成虛數,如 t = iu 並且一併做 s=c+ix 的變換,則我們可以將上式轉換為如下的 拉普拉斯變換#雙邊拉普拉斯變換

12πg(c+ix)euxeicudx.

這裡就可以套用此拉普拉斯方法求漸進解,最後,再利用 u = t / it 換回來,就可以得到該拉普拉斯逆變換的漸進解了。

例子1:斯特靈公式

拉普拉斯方法可以用在推導 斯特靈公式 上;當 N 很大時,

N!2πNNNeN

證明:

Γ函数 的積分定義,我們可以得到

N!=Γ(N+1)=0exxNdx.

接著讓我們做變數變換,

x=Nz

因此

dx=Ndz.

將這些代回 Γ函数 的積分定義裡,我們可以得到

N!=0eNz(Nz)NNdz=NN+10eNzzNdz=NN+10eNzeNlnzdz=NN+10eN(lnzz)dz.

經由此變數變換後,我們有了拉普拉斯方法所需要的f(x)

f(z)=lnzz

而它乃為二次可微函數,且

f(z)=1z1,
f(z)=1z2.

因此, ƒ(z) 的極大值出現在 z0 = 1 而且在該點的二次微分為 1 。因此,我們得到

N!NN+12πNeN=2πNNNeN.

例子2:貝氏網路 ,參數估計與概率推理

關於概率推理的簡介,請參考 -{R|http://doc.baidu.com/view/e9c1c086b9d528ea81c77989.html}- 。 而在 Template:Harvnb 的文章裡,則回顧了如何應用此拉普拉斯方法 (無論是單變數或者多變數) 如何應用在 概率推理 上,以加速得到系統的 後驗矩 (數學) (posterior moment) , Template:Tsl 等,並舉醫學診斷上的應用為例。

相關維基百科文章

參考文獻

Template:Reflist