李鑫宇, 陳旭梅, 岳 華

(1.吉林大學(xué) 數(shù)學(xué)學(xué)院, 長春 130012; 2.江蘇大學(xué) 理學(xué)院, 江蘇 鎮(zhèn)江 212013)

在實際應(yīng)用中, 許多物理、 化學(xué)等領(lǐng)域的微分方程系統(tǒng)都會受隨機擾動的影響, 因此考慮具有隨機項的微分方程對實際應(yīng)用有重要意義.文獻[9]基于隨機向量場的分裂法, 對于一般Stratonovich意義下的具有守恒量的隨機微分方程提出了一種新的保守方法, 在求解過程中能保持系統(tǒng)中的守恒量, 且該方法在噪聲為可交換時可達到均方1階收斂; 文獻[10]對于Stratonovich意義下的具有守恒量的隨機微分方程應(yīng)用分裂方法求解, 不僅使計算過程中保持守恒量, 而且可以達到與已有的Milstein方法相同的收斂階, 同時證明了隨機微分方程隱格式的收斂性；文獻[11]給出了標(biāo)準(zhǔn)布朗運動要滿足的3個條件, 同時基于Euler-Maruyama方法和Milstein 方法等給出了數(shù)值實例.

近年來, 將分裂方法應(yīng)用于求解具有一些特殊性質(zhì)的物理系統(tǒng)受到廣泛關(guān)注, 如Hong等[12]對于一類具有能量耗散性質(zhì)的二維Maxwell方程組提出了分裂有限差分時域法, 該方法能保持系統(tǒng)能量耗散的特性, 同時是無條件穩(wěn)定的; Cai等[13]對于三維時域Maxwell方程組提出了一種指數(shù)算子分裂方法, 該方法具有保能量、 高精度、 無條件穩(wěn)定且節(jié)約計算資源等優(yōu)點; Li等[14]對于非線性Dirac方程提出了4種時間分裂方法, 這些方法都能滿足離散能級上的電荷守恒, 且節(jié)省了運算時間.本文對一類有守恒量的高維隨機Hamilton系統(tǒng)提出一種具有降維和簡化運算效果的分裂求解算法.該方法主要利用中點公式構(gòu)造一種對稱的分裂算法, 借助Milstein公式和局部收斂與整體收斂的關(guān)系定理得到該方法在均方意義下的整體收斂階.數(shù)值算例結(jié)果表明了該方法的有效性.

1 隨機Hamilton系統(tǒng)的分裂算法

考慮如下2n維Stratonovich型隨機Hamilton系統(tǒng):

Y: dy(t)=SH(y)dt+TH(y)°dw(t)，

(1)

其中:y(t)=(y(1)(t),y(2)(t))T∈n×n;S和T是兩個斜對稱矩陣;w(t)是標(biāo)準(zhǔn)Wiener過程;H(y)是能量函數(shù), 也稱為Hamilton函數(shù), 其為滿足下述方程的守恒量:

HTSH=0,HTTH=0.

為表示方便, 記

iH∶=y(i)H,iS∶=y(i)S,iT∶=y(i)T,i=1,2.

從而可將方程(1)寫成如下形式:

(2)

其中S11,S22,T11和T22均為斜對稱矩陣, 且(S12)T=-S21, (T12)T=-T21.

基于對稱分裂算法, 可將原系統(tǒng)的向量場Y分解為如下兩個子向量場:

(3)

(4)

考慮將方程(1)在時間區(qū)間[t0,T]上進行均勻剖分, 步長為h, 剖分節(jié)點為tn=t0+nh,n=0,1,…,N, 并記Δw=w(tn+1)-w(tn).在單個時間區(qū)間[tn,tn+1]上, 假設(shè)y(tn)是方程(1)在tn處的精確解, 用yn表示解y(tn)的數(shù)值近似.對式(3)用中點格式離散, 得數(shù)值近似格式為

其中：

(5)

對式(4)用中點格式離散得數(shù)值近似格式為

其中：

(6)

基于分裂法[4]的基本思想, 構(gòu)造方程(1)在時間區(qū)間[tn,tn+1]上的一步近似格式:

(7)

式(7)是一個對稱格式.

2 整體收斂性分析

定理1對于2n維Sratonovich型隨機Hamilton系統(tǒng)(1), 分裂算法格式(7)在區(qū)間[t0,T]上是整體均方一階收斂的.

證明：在區(qū)間[tn,tn+1]上, 假設(shè)yn=y(tn).設(shè)yM是以yn為初始狀態(tài)方程(1)的Milstein格式的一步迭代數(shù)值近似, 對于k=1,2, 寫成如下的分量形式:

(8)

文獻[12]給出了Milstein格式的誤差階是均方一階收斂的.為了證明數(shù)值格式(7)的誤差階也是均方一階的收斂性, 只需證明

(9)

(10)

(11)

進一步做類似地Taylor展開, 可得分裂算法格式(7)的一步近似分量方程為

(12)

EΔw=0,E(Δw)2=h,E(Δw)3=0,E(Δw)4=3h2,

可得

‖EΔi‖=O(h2), (E‖Δi‖2)1/2=O(h1/2),i=1,2.

(13)

因此, 有如下估計：

(14)

從而由式(8),(12)可以推出

(15)

進而, 式(9)可由式(15)和Minkowski不等式得到.

進一步, 由文獻[15]中定理1給出的局部收斂與整體收斂的關(guān)系, 可知該分裂算法在區(qū)間[t0,T]上的精度p=3/2-1/2=1, 即該方法在均方意義下的整體收斂階是1階的.證畢.

3 數(shù)值算例

在下面數(shù)值算例的計算中, 由于沒有精確解, 因此將所取最小步長下求得的數(shù)值解作為參考精確解y(tn).定義誤差函數(shù)

Hamilton函數(shù)的誤差為

在計算過程中, 固定hn=1.為了得到平均意義下的收斂速率, 對于每個算例, 隨機選取 1 000條樣本軌道, 然后計算平均值, 從而得到數(shù)值解相對于精確解在均方意義下的截斷誤差.

例1(隨機簡諧振蕩系統(tǒng)) 考慮一個具有一維Wiener過程的Stratonovich型隨機簡諧振蕩系統(tǒng):

(16)

在常系統(tǒng)中它是一個守恒量.易知該系統(tǒng)中S和T是兩個斜對稱矩陣.應(yīng)用上述分裂算法數(shù)值求解該系統(tǒng)的近似解流.

圖1 隨機簡諧振蕩系統(tǒng)的截斷誤差e(tn)和eH(tn)與步長h的關(guān)系

例2(隨機Kepler系統(tǒng)) 具有一維Wiener過程的隨機Kepler系統(tǒng)在Stratonovich意義下的方程為

(17)

該系統(tǒng)是一個四維的隨機Hamilton系統(tǒng), 其中Hamilton函數(shù)為

易知該系統(tǒng)中的S和T均為斜對稱矩陣.基于分裂算法的主要思想, 可以將該隨機Hamilton系統(tǒng)分解為兩個二維子系統(tǒng).

為了驗證應(yīng)用分裂算法數(shù)值求解隨機Kepler系統(tǒng)的精確性, 觀察如圖2所示的隨機Kepler系統(tǒng)的截斷誤差e(tn)和eH(tn)與步長h的關(guān)系, 其中隨機Kepler系統(tǒng)的初值為

由圖2(A)可見, 該系統(tǒng)解的誤差是一階收斂的；由圖2(B)可見, Hamilton函數(shù)的誤差也是一階收斂的.

圖2 隨機Kepler系統(tǒng)的截斷誤差e(tn)和eH(tn)與步長h的關(guān)系

綜上, 通過將本文構(gòu)造的分裂算法應(yīng)用于求解隨機簡諧振蕩系統(tǒng)和隨機Kepler系統(tǒng), 驗證了理論分析的結(jié)果, 表明該分裂算法對于數(shù)值求解一類隨機微分方程十分有效.