UA MATH575B 数值分析下VI 统计物理的随机模拟方法1
UA MATH575B 數值分析下VI 統計物理的隨機模擬方法1
- 模擬SDE的解
- 用蒙特卡羅方法估計SDE解的自相關函數
模擬SDE的解
以最簡單的隨機微分方程為例,考慮隨機游走
dxt=ξtdtdx_t = \xi_t dtdxt?=ξt?dt
其中ξt\xi_tξt?是白噪聲,假設?ξt1,ξt2?=2Dδ(t1?t2)\langle \xi_{t_1},\xi_{t_2} \rangle = 2D\delta(t_1-t_2)?ξt1??,ξt2???=2Dδ(t1??t2?),其中D>0D>0D>0,δ\deltaδ是Dirac函數。假設隨機過程xtx_txt?的概率密度函數為ρ(t,x)\rho(t,x)ρ(t,x),它滿足Fokker-Planck方程:
?ρ?t=D?2ρ?x2\frac{\partial \rho}{\partial t} = D \frac{\partial^2 \rho}{\partial x^2}?t?ρ?=D?x2?2ρ?
這個是一維的擴散方程,DDD叫做擴散系數。這個方程的一個解是
ρ(t,x)=14πDtexp?(?x24Dt)\rho(t,x) = \frac{1}{\sqrt{4\pi Dt}}\exp \left(- \frac{x^2}{4Dt} \right)ρ(t,x)=4πDt?1?exp(?4Dtx2?)
即xtx_txt?是高斯過程,分布可以記為N(0,2Dt)N(0,2Dt)N(0,2Dt)。
現在對時間做離散化,記time step為τ\tauτ,簡單回顧一下解數值ODE的Runge-Kutta方法 (RK):對于動力系統x˙=f(t,x)∈Rm\dot{x}=f(t,x)∈\mathbb{R}^mx˙=f(t,x)∈Rm,定義RK的系數為A∈Rs×s,b∈Rs,c∈RsA∈\mathbb{R}^{s×s}, b∈\mathbb{R}^s,c∈\mathbb{R}^sA∈Rs×s,b∈Rs,c∈Rs
ki=hf(t+cih,x(t)+∑j=1sAijkj)x(t+h)=x(t)+∑j=1sbjkjk_i=hf(t+c_i h,x(t)+\sum_{j=1}^s A_{ij} k_j) \\ x(t+h)=x(t)+\sum_{j=1}^s b_j k_jki?=hf(t+ci?h,x(t)+j=1∑s?Aij?kj?)x(t+h)=x(t)+j=1∑s?bj?kj?
其中sss是RK方法的階數,hhh是步長,根據上述遞推關系,可以計算出x(t)x(t)x(t)的序列。
對于隨機微分方程dxt=ξtdtdx_t = \xi_t dtdxt?=ξt?dt,考慮RK2:
A=[00120],b=[0,1]T,c=[0,12]T,h=τA = \left[ \begin{matrix} 0 & 0 \\ \frac{1}{2} & 0 \end{matrix} \right],b=[0,1]^T,c=[0,\frac{1}{2}]^T,h=\tauA=[021??00?],b=[0,1]T,c=[0,21?]T,h=τ
則
k1=τξi,k2=τξi+12xi+1=xi+τξi+12k_1 = \tau \xi_i,k_2 = \tau \xi_{i+\frac{1}{2}} \\ x_{i+1} = x_i + \tau \xi_{i+\frac{1}{2}}k1?=τξi?,k2?=τξi+21??xi+1?=xi?+τξi+21??
xi+1=xi+τξi+12x_{i+1} = x_{i} + \tau \xi_{i+\frac{1}{2}}xi+1?=xi?+τξi+21??
其中?ξi,ξj?=2Dδij/τ\langle \xi_{i},\xi_{j} \rangle = 2D\delta_{ij}/\tau?ξi?,ξj??=2Dδij?/τ,其中δij\delta_{ij}δij?為Kronecker符號,所有ξi\xi_iξi?都是獨立同分布的。因此可以將上式寫成:
xi+1=xi+2Dτζζ~N(0,1)x_{i+1} = x_{i} + \sqrt{2D\tau} \zeta\\\zeta \sim N(0,1)xi+1?=xi?+2Dτ?ζζ~N(0,1)
要解這個隨機過程可以用隨機模擬的方法:
重復這兩步,直到達到設置的時間上限。
由此我們可以總結出模擬一個隨機微分方程的一般性做法:
例1下面的SDE表示幾何隨機游走
dxt=xtξtdtdx_t =x_t \xi_t dtdxt?=xt?ξt?dt
其中ξt\xi_tξt?是白噪聲,假設?ξt1,ξt2?=2Dδ(t1?t2)\langle \xi_{t_1},\xi_{t_2} \rangle = 2D\delta(t_1-t_2)?ξt1??,ξt2???=2Dδ(t1??t2?),D>0D>0D>0。用這個例子介紹兩種離散化的規則:
在這個例子中,Ito規則導出的結果是Exi+1=Exi=x0Ex_{i+1}=Ex_i = x_0Exi+1?=Exi?=x0?,這個比較明顯。但Stratonovich規則導出的是
xi+1=xi1+2Dτζ21?2Dτζ2=xi(1+2Dτζ2)(1+2Dτζ2+2Dτζ24+o(ζ2))=xi(1+2Dτζ+Dτζ2+o(ζ2))→xi(1+Dτ+2Dτζ)x_{i+1} = x_i \frac{1+ \sqrt{2D\tau} \frac{\zeta}{2}}{1- \sqrt{2D\tau} \frac{\zeta}{2}} = x_i(1+ \sqrt{2D\tau} \frac{\zeta}{2})(1+\sqrt{2D\tau} \frac{\zeta}{2}+2D\tau \frac{\zeta^2}{4}+o(\zeta^2)) \\ = x_i(1+\sqrt{2D\tau}\zeta + D\tau \zeta^2 + o(\zeta^2)) \to x_i(1+D\tau+\sqrt{2D\tau}\zeta )xi+1?=xi?1?2Dτ?2ζ?1+2Dτ?2ζ??=xi?(1+2Dτ?2ζ?)(1+2Dτ?2ζ?+2Dτ4ζ2?+o(ζ2))=xi?(1+2Dτ?ζ+Dτζ2+o(ζ2))→xi?(1+Dτ+2Dτ?ζ)
因此
Exi+1=Exi(1+Dτ)Ex_{i+1} = Ex_i (1+D \tau)Exi+1?=Exi?(1+Dτ)
因此這兩個規則模擬出來的結果會有明顯差別。
用蒙特卡羅方法估計SDE解的自相關函數
隨著時間流逝,xtx_txt?信息減少的規律可以用自相關函數來表示,
KXX(τ)=?(xt+τ?xt)2?=∫tt+τdt1∫tt+τdt2?ξt1ξt2?=2DτK_{XX}(\tau) = \langle (x_{t+\tau}-x_t)^2 \rangle = \int_t^{t+\tau} dt_1 \int_t^{t+\tau} dt_2 \langle \xi_{t_1} \xi_{t_2}\rangle =2D \tauKXX?(τ)=?(xt+τ??xt?)2?=∫tt+τ?dt1?∫tt+τ?dt2??ξt1??ξt2???=2Dτ
隨機游走的自相關函數顯然只與兩個時點的差有關,這種隨機過程叫平穩隨機過程。如果SDE的解比較復雜,甚至根本寫不出解析解,很難直接計算自相關函數,我們可以用Monte Carlo方法估計自相關函數:
KXX(τ)=Extxt+τ≈1T∑i=1Txixi+τK_{XX}(\tau) = Ex_tx_{t+\tau} \approx \frac{1}{T} \sum_{i=1}^T x_i x_{i+\tau}KXX?(τ)=Ext?xt+τ?≈T1?i=1∑T?xi?xi+τ?
其中xi,xi+τx_i,x_{i+\tau}xi?,xi+τ?可以用隨機模擬得到。
例2 一個簡單的單位根過程,假設擴散系數D=12D=\frac{1}{2}D=21?
dxt=(?xt+ξt)dtdx_t = (-x_t + \xi_t )dtdxt?=(?xt?+ξt?)dt
參考SDE的第一講,這個方程定義了一個Ornstein-Uhlenbeck過程,它的解是
xt=∫?∞te?(t?s)ξsdsx_t = \int_{-\infty}^t e^{-(t-s)}\xi_s dsxt?=∫?∞t?e?(t?s)ξs?ds
根據這個解計算自相關函數,第一行到第二行用到的性質是Ito Isometry:
KXX(τ)=Extxt+τ=E∫?∞te?(t?s)ξsds∫?∞t+τe?(t+τ?r)ξrdr=∫?∞t∫?∞t+τe?(t?s)e?(t+τ?r)Eξsξrdsdr=∫?∞t∫?∞t+τe?(t?s)e?(t+τ?r)δ(r?s)dsdr=12exp?(?τ),τ≥0K_{XX}(\tau) = Ex_{t}x_{t+\tau} = E \int_{-\infty}^t e^{-(t-s)}\xi_s ds \int_{-\infty}^{t+\tau} e^{-(t+\tau-r)}\xi_r dr \\ =\int_{-\infty}^t \int_{-\infty}^{t+\tau} e^{-(t-s)}e^{-(t+\tau-r)} E\xi_s\xi_r dsdr \\ =\int_{-\infty}^t \int_{-\infty}^{t+\tau} e^{-(t-s)}e^{-(t+\tau-r)} \delta(r-s) dsdr = \frac{1}{2} \exp(-\tau),\tau\ge0 KXX?(τ)=Ext?xt+τ?=E∫?∞t?e?(t?s)ξs?ds∫?∞t+τ?e?(t+τ?r)ξr?dr=∫?∞t?∫?∞t+τ?e?(t?s)e?(t+τ?r)Eξs?ξr?dsdr=∫?∞t?∫?∞t+τ?e?(t?s)e?(t+τ?r)δ(r?s)dsdr=21?exp(?τ),τ≥0
如果我們不知道Ito Isometry,要計算自相關函數的解析解就會比較麻煩,這時就可以用Monte Carlo方法來估計自相關函數。先用RK2得到SDE的數值算法:
xi+1=e?τxi+τζ≈(1?τ)xi+τζx_{i+1}=e^{-\tau}x_i + \sqrt{\tau} \zeta \approx (1-\tau)x_i + \sqrt{\tau} \zeta xi+1?=e?τxi?+τ?ζ≈(1?τ)xi?+τ?ζ
根據這個遞推關系模擬SDE的解,然后根據
KXX(τ)=Extxt+τ≈1T∑i=1Txixi+τK_{XX}(\tau) = Ex_tx_{t+\tau} \approx \frac{1}{T} \sum_{i=1}^T x_i x_{i+\tau}KXX?(τ)=Ext?xt+τ?≈T1?i=1∑T?xi?xi+τ?
估計自相關函數。
總結
以上是生活随笔為你收集整理的UA MATH575B 数值分析下VI 统计物理的随机模拟方法1的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: UA MATH566 统计理论7 一个例
- 下一篇: UA MATH564 概率论II 连续型