劉永財(cái)，柏明花

(曲靖師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，云南 曲靖 655011)

0 引 言

隨機(jī)微分方程中隨機(jī)因素的存在使得求解方程的過程非常復(fù)雜，很多隨機(jī)微分方程很難解出精確解，故數(shù)值求解具有重要意義，國(guó)內(nèi)外涌現(xiàn)出大批學(xué)者對(duì)隨機(jī)微分方程的數(shù)值法求解進(jìn)行研究.劉濤[1]對(duì)隨機(jī)微分方程的數(shù)值方法進(jìn)行了系統(tǒng)的研究；王彩霞[2]研究了隨機(jī)微分方程數(shù)值方法的穩(wěn)定性與收斂性；李煥榮[3]不僅對(duì)隨機(jī)微分方程數(shù)值求解方法進(jìn)行研究，還研究了其應(yīng)用；謝曉波[4]對(duì)幾類隨機(jī)微分方程與偏微分方程的數(shù)值解法進(jìn)行了研究；陳晨[5]研究了隨機(jī)微分方程解的存在唯一性與其數(shù)值方法的穩(wěn)定性；張永強(qiáng),盧俊香[6]研究了一類線性隨機(jī)微分方程的解法；周迎春[7]研究了幾種隨機(jī)微分方程數(shù)值方法以及數(shù)值模擬；葉安[8]深入研究了伊藤隨機(jī)微分方程的兩種數(shù)值方法；汪勇[9]對(duì)Euler-Maruyama方法進(jìn)行了細(xì)致的研究；Li Lei等[10]研究了基于高斯混合的隨機(jī)微分方程的數(shù)值方法；Mustafa Bayram[11]研究了隨機(jī)微分方程的數(shù)值模擬方法.本文在學(xué)者們的研究基礎(chǔ)之上，對(duì)顯式Euler方法、梯形Euler方法及顯式Milstein方法進(jìn)行研究，通過觀察不同步長(zhǎng)時(shí)的誤差和運(yùn)算時(shí)間來分析數(shù)值解的精度.

1 隨機(jī)微分方程及其數(shù)值方法

1.1 It型隨機(jī)微分方程

布朗(Brown)通過顯微鏡對(duì)液體中的花粉做出觀察，發(fā)現(xiàn)它們作不規(guī)則運(yùn)動(dòng).下面給出布朗運(yùn)動(dòng)的定義：

定義1 一個(gè)在[0，T]上的隨機(jī)過程W(t),如果它連續(xù)的依賴于t∈[0,T]并且滿足以下的三個(gè)條件：

(1)W(0)=0；

(2)對(duì)于0≤s

(3)對(duì)于0≤s

標(biāo)準(zhǔn)布朗運(yùn)動(dòng)的性質(zhì)：

(1)P(W(0)=0)=1;

(2)?t>0,h>0,△W(t)=W(t+h)-W(t)～N(0,h)

(3)?t∈T,W(t)不可微.

由上可知，Brown運(yùn)動(dòng)處處連續(xù)，處處不可微，有著曲折的運(yùn)動(dòng)軌跡.

Wi=Wi-1+dWi,i=1,2,…,n,

圖1 布朗運(yùn)動(dòng)軌跡圖

常微分方程的一般形式是:

dX(t)=X(t)dt,t∈[t0,T]

(1)

隨機(jī)微分方程與常微分方程最大的不同之處在于隨機(jī)微分方程里面存在著隨機(jī)因素，因此方程的求解就有不確定性，常微分方程并不能有效地描述這種不確定性系統(tǒng)，所以要對(duì)方程(1)進(jìn)行適當(dāng)?shù)恼{(diào)整，即加入隨機(jī)因素，使它能夠把帶隨機(jī)干擾的現(xiàn)象反映出來，It型隨機(jī)微分方程因?yàn)榫哂绪毙再|(zhì)而讓其計(jì)算更加簡(jiǎn)便.

定義2 假設(shè){X(t,w),0≤t≤T}是一個(gè)隨機(jī)過程，{W(t),t≥0}是標(biāo)準(zhǔn)的Brown運(yùn)動(dòng)，對(duì)任一個(gè)劃分

P={0=t0

對(duì)于如下形式的隨機(jī)微分方程：

t∈[t0,T]

(2)

其中f:R+×Rn→Rn稱為漂移系數(shù)，矩陣函數(shù)g:R+×Rn→Rn×m稱為擴(kuò)散系數(shù)，W(t)是m維Brown運(yùn)動(dòng)，X0是隨機(jī)變量.

根據(jù)布朗運(yùn)動(dòng)的連續(xù)性可以把(2)寫為下面的形式：

0≤t0≤t≤T.

1.2 Euler方法和Milstein方法

大部分隨機(jī)微分方程的解析解是無法獲得的，在實(shí)際應(yīng)用中，實(shí)用的方法是在計(jì)算機(jī)上進(jìn)行數(shù)值求解，即不直接求出y(t)的解析解，而是在隨機(jī)微分方程的解存在的區(qū)間上，求得一系列點(diǎn)xn(n=1,2,…)上的近似值[5].到目前為止，有多種數(shù)值方法來求解隨機(jī)微分方程，Euler方法以及Milstein方法是對(duì)隨機(jī)微分方程數(shù)學(xué)模型進(jìn)行數(shù)值模擬的兩種便捷、高效的方法.

(1)顯式Euler方法：

Xn+1=Xn+f(Xn)h+g(Xn)△Wn

(3)

其中，Xn是在tN=nh處的近似解，△Wn=Wtn+1-Wtn表示定義在區(qū)間[tn,tn+1]上的Wiener過程的增量，是服從N(0,△t)分布的相互獨(dú)立的隨機(jī)變量.

(2)梯形Euler-Maruyama方法：

(4)

而Milstein方法是求解隨機(jī)常微分方程的具有一階強(qiáng)收斂的方法[3]，顯式Milstein方法如下：

(5)

2 數(shù)值算例

以如下齊次線性隨機(jī)微分方程為例：

(6)

其中，X(t)是隨機(jī)過程，λ和μ都是常數(shù)，W(t)是標(biāo)準(zhǔn)布朗運(yùn)動(dòng).f(X)=λX(t)為漂移系數(shù)；g(X)=μX(t)為擴(kuò)散系數(shù).隨機(jī)微分方程可能有無窮多個(gè)解，為了限定解，需要加入一個(gè)初值條件：X(0)=X0.

令λ=1，μ=1，X0=1，來求此方程的解析解，求解步驟如下：

(2)使用伊藤引理，得到d lnX(t)的隨機(jī)微分方程(SDE)：

綜上可以得出，類似方程(6)的隨機(jī)微分方程的解析解通式為：

對(duì)于方程(6)，令λ=1，μ=1，X0=1，取N=2-11來估計(jì)方程的解析解，取時(shí)間步長(zhǎng)h分別為2-2，2-5，2-8，2-11，利用MATLAB編碼可以得到方程的解析解以及顯式Euler方法求得的數(shù)值解，如圖2所示：

圖2 顯式Euler法數(shù)值運(yùn)算結(jié)果 (步長(zhǎng)為2-2，2-5，2-8，2-11)

由圖2可以得到:時(shí)間步長(zhǎng)h依次縮小8倍，顯式Euler方法得到的數(shù)值解越來越接近解析解.還可以看到t的范圍在0到0.75之間時(shí)，模擬出來的數(shù)值解與解析解誤差是很大的，但是t的范圍在0.75到1時(shí)，解析解與數(shù)值解相比于之前較為接近，誤差越小.

通過計(jì)算可以得到：步長(zhǎng)為2-2，2-5，2-8，2-11時(shí)，平均絕對(duì)誤差分別為0.0181，0.0144，0.0128，0.0103，這說明隨著步長(zhǎng)逐漸減小，得出的數(shù)值解與解析解的誤差越來越小，逼近效果也漸漸變好.

對(duì)于方程(6)，取λ=0.5，μ=1，X0=1，取T=1,h=2-8模擬數(shù)值解，分別對(duì)顯式Euler方法和梯形Euler方法得到的數(shù)值解與方程的解析解進(jìn)行對(duì)比，結(jié)果如圖3和圖4所示.

由圖3和圖4可以發(fā)現(xiàn)，顯式Euler方法與梯形Euler方法相比，用梯形Euler方法得到的數(shù)值解更加接近解析解，且用梯形Euler方法模擬數(shù)值解時(shí)，t的范圍在0.7至1之間模擬效果最好，數(shù)值解與解析解差距最小.另外顯式Euler方法與梯形Euler方法的誤差分別是0.4715和0.0406，表明梯形Euler方法比顯式Euler方法得到的數(shù)值解更加精確.

圖3 顯式Euler數(shù)值解與解析解 圖4 梯形Euler數(shù)值解與解析解

令λ=1，μ=0.3，X0=2，N=2-11，取時(shí)間步長(zhǎng)h分別為2-2，2-5，2-8，2-11，模擬方程(6)的數(shù)值解，結(jié)果如圖5所示.

可以看出，隨著步長(zhǎng)從2-2逐漸減小到2-11，數(shù)值解與解析解的誤差也越來越小.另外，t在0.6至1之間時(shí)，數(shù)值解與解析解的擬合效果最好.

下面將對(duì)Milstein方法與Euler方法進(jìn)行比較.對(duì)于隨機(jī)微分方程(6)，令λ=2，μ=1，X0=1，分別計(jì)算步長(zhǎng)為2-3，2-6，2-9，2-12時(shí)，隨機(jī)微分方程解析解、Euler方法模擬數(shù)值解以及Milstein方法模擬數(shù)值解的運(yùn)算時(shí)間，統(tǒng)計(jì)結(jié)果如表1所示；而Euler方法、Milstein方法的平均絕對(duì)誤差，計(jì)算結(jié)果如表2所示.

由表1可以看出，時(shí)間步長(zhǎng)越大，運(yùn)算時(shí)間越短，但是模擬的結(jié)果越不精確.同時(shí)由表2可以看出：隨著步長(zhǎng)逐漸縮小，Euler方法和Milstein方法的平均絕對(duì)誤差也慢慢減小，且Milstein方法的平均絕對(duì)誤差比Euler方法的平均誤差要小，說明了Milstein方法的數(shù)值解更加精確.

3 總 結(jié)

本文對(duì)齊次線性隨機(jī)微分方程數(shù)值求解的Euler方法與Milstein方法進(jìn)行研究分析，我們發(fā)現(xiàn)，顯式Euler方法與Milstein方法都能夠準(zhǔn)確模擬出隨機(jī)微分方程的解，而 Milstein方法得出的數(shù)值解的精度最高，梯形Euler方法精度次之，顯式Euler方法精度最低；通過計(jì)算能夠得出顯式Euler方法與Milstein方法的運(yùn)算時(shí)間，且運(yùn)算時(shí)間越長(zhǎng)，得到的數(shù)值解越精確.

圖5 顯式Milstein法數(shù)值運(yùn)算結(jié)果 (步長(zhǎng)為2-2，2-5，2-8，2-11)

表1 運(yùn)算時(shí)間計(jì)算結(jié)果

表2 Euler方法與Milstein方法的平均絕對(duì)誤差