UA MATH575B 數值分析下VI 統計物理的隨機模擬方法1

• 模擬SDE的解
• 用蒙特卡羅方法估計SDE解的自相關函數

模擬SDE的解

以最簡單的隨機微分方程為例，考慮隨機游走
$$d{x}_t={\xi }_{t}dt$$
其中${\xi }_t$是白噪聲，假設$?{\xi }_{t_1},{\xi }_{t_2}?=2D\delta (t_1?t_2)$，其中$D>0$，$\delta$是Dirac函數。假設隨機過程$x_t$的概率密度函數為$\rho (t,x)$，它滿足Fokker-Planck方程：
$$\frac{?\rho }{?t}=D\frac{{?}^2\rho }{?x^2}$$
這個是一維的擴散方程，$D$叫做擴散系數。這個方程的一個解是
$$\rho (t,x)=\frac{1}{\sqrt{4\pi Dt}}\exp \left(?\frac{x^2}{4Dt}\right)$$
即$x_t$是高斯過程，分布可以記為$N(0,2Dt)$。

現在對時間做離散化，記time step為$\tau$，簡單回顧一下解數值ODE的Runge-Kutta方法 (RK)：對于動力系統$\dot{x}=f(t,x)\in {\mathbb{R}}^m$，定義RK的系數為$A\in {\mathbb{R}}^{s\times s},b\in {\mathbb{R}}^s,c\in {\mathbb{R}}^s$
$$\begin{gathered} k_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}_j \end{gathered}$$
其中$s$是RK方法的階數，$h$是步長，根據上述遞推關系，可以計算出$x(t)$的序列。

對于隨機微分方程$d{x}_t={\xi }_{t}dt$，考慮RK2：
$$A=\left[\begin{array}{cc}0& 0\\ \frac{1}{2}& 0\end{array}\right],b=[0,1{]}^T,c=[0,\frac{1}{2}{]}^T,h=\tau$$
則
$$\begin{gathered} k_1=\tau {\xi }_i,k_2=\tau {\xi }_{i+\frac{1}{2}} \\ {x}_{i+1}=x_i+\tau {\xi }_{i+\frac{1}{2}} \end{gathered}$$
$$x_{i+1}=x_i+\tau {\xi }_{i+\frac{1}{2}}$$
其中$?{\xi }_i,{\xi }_j?=2D{\delta }_{ij}/\tau$，其中${\delta }_{ij}$為Kronecker符號，所有${\xi }_i$都是獨立同分布的。因此可以將上式寫成：
$$\begin{gathered} x_{i+1}=x_i+\sqrt{2D\tau }\zeta \\ \zeta ～N(0,1) \end{gathered}$$

要解這個隨機過程可以用隨機模擬的方法：

第一步，從標準正態分布中抽一個樣本$\zeta$；

第二步，根據$x_{i+1}=x_i+\sqrt{2D\tau }\zeta$更新隨機過程的值；

重復這兩步，直到達到設置的時間上限。

由此我們可以總結出模擬一個隨機微分方程的一般性做法：

用數值ODE的方法導出SDE解的遞推關系

生成正態分布的隨機數，根據遞推關系模擬這個SDE

例1下面的SDE表示幾何隨機游走
$$d{x}_t=x_t{\xi }_{t}dt$$
其中${\xi }_t$是白噪聲，假設$?{\xi }_{t_1},{\xi }_{t_2}?=2D\delta (t_1?t_2)$，$D>0$。用這個例子介紹兩種離散化的規則：

Ito規則：$x_{i+1}=x_i+\sqrt{2D\tau }{x}_i\zeta$

Stratonovich規則：$x_{i+1}=x_i+\sqrt{2D\tau }\frac{x_i+x_{i+1}}{2}\zeta$

在這個例子中，Ito規則導出的結果是$E{x}_{i+1}=E{x}_i=x_0$，這個比較明顯。但Stratonovich規則導出的是
$$\begin{gathered} 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 ) \end{gathered}$$
因此
$$E{x}_{i+1}=E{x}_i(1+D\tau )$$
因此這兩個規則模擬出來的結果會有明顯差別。

用蒙特卡羅方法估計SDE解的自相關函數

隨著時間流逝，$x_t$信息減少的規律可以用自相關函數來表示，
$$K_{XX}(\tau )=?(x_{t+\tau }?x_t{)}^2?={\int }_t^{t+\tau }d{t}_1{\int }_t^{t+\tau }d{t}_2?{\xi }_{t_1}{\xi }_{t_2}?=2D\tau$$
隨機游走的自相關函數顯然只與兩個時點的差有關，這種隨機過程叫平穩隨機過程。如果SDE的解比較復雜，甚至根本寫不出解析解，很難直接計算自相關函數，我們可以用Monte Carlo方法估計自相關函數：
$$K_{XX}(\tau )=E{x}_t{x}_{t+\tau }\approx \frac{1}{T}\sum _{i=1}^T{x}_i{x}_{i+\tau }$$
其中$x_i,x_{i+\tau }$可以用隨機模擬得到。

例2 一個簡單的單位根過程，假設擴散系數$D=\frac{1}{2}$
$$d{x}_t=(?x_t+{\xi }_t)dt$$
參考SDE的第一講，這個方程定義了一個Ornstein-Uhlenbeck過程，它的解是
$$x_t={\int }_{?\infty }^t{e}^{?(t?s)}{\xi }_{s}ds$$
根據這個解計算自相關函數，第一行到第二行用到的性質是Ito Isometry：
$$\begin{gathered} K_{XX}(\tau )=E{x}_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 \ge 0 \end{gathered}$$

如果我們不知道Ito Isometry，要計算自相關函數的解析解就會比較麻煩，這時就可以用Monte Carlo方法來估計自相關函數。先用RK2得到SDE的數值算法：
$$x_{i+1}=e^{?\tau }{x}_i+\sqrt{\tau }\zeta \approx (1?\tau )x_i+\sqrt{\tau }\zeta$$
根據這個遞推關系模擬SDE的解，然后根據
$$K_{XX}(\tau )=E{x}_t{x}_{t+\tau }\approx \frac{1}{T}\sum _{i=1}^T{x}_i{x}_{i+\tau }$$
估計自相關函數。

總結

以上是生活随笔為你收集整理的UA MATH575B 数值分析下VI 统计物理的随机模拟方法1的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。