李 煥 榮

(重慶工商大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，重慶 400067)

0 引 言

近十年以來(lái)，隨機(jī)數(shù)學(xué)模型[1]的研究在很多歐美國(guó)家已經(jīng)得到了高度重視和大力發(fā)展, 吸引了大批應(yīng)用數(shù)學(xué)家和計(jì)算數(shù)學(xué)家的關(guān)注, 隨機(jī)數(shù)學(xué)模型的數(shù)值求解方法[2-4]和理論分析也得到了快速的發(fā)展。目前, 隨機(jī)數(shù)學(xué)模型在美國(guó)已經(jīng)成為最重要的應(yīng)用數(shù)學(xué)和計(jì)算數(shù)學(xué)的研究方向之一。美國(guó)很多部門(mén)比如能源部、空軍和國(guó)家實(shí)驗(yàn)室，均設(shè)立了專項(xiàng)基金來(lái)支持和發(fā)展隨機(jī)數(shù)學(xué)模型的數(shù)值計(jì)算、理論分析和相關(guān)應(yīng)用等研究工作。

隨機(jī)常微分方程(Stochastic Ordinary Differential Equation, 簡(jiǎn)寫(xiě)為SODE)是近幾十年來(lái)才興起的熱門(mén)邊緣學(xué)科，是概率論與常微分方程(Ordinary Differential Equation, 簡(jiǎn)寫(xiě)為ODE)相結(jié)合的產(chǎn)物，自日本數(shù)學(xué)家伊藤清[5](It)先生于20世紀(jì)40年代創(chuàng)立隨機(jī)積分(It積分)以來(lái)，隨機(jī)微分方程很快就被廣泛應(yīng)用到金融經(jīng)濟(jì)、自然科學(xué)和工程技術(shù)等很多領(lǐng)域，如股市預(yù)測(cè)、資料同化、油藏模擬、天氣預(yù)測(cè)、流體動(dòng)力學(xué)和生物遺傳學(xué)等。湯濤院士等[1]從不確定性量化的角度討論了隨機(jī)數(shù)學(xué)模型的發(fā)展近況； Markos[2]討論了隨機(jī)相場(chǎng)模型Allen-Cahn的噪聲正則性和相關(guān)計(jì)算；伊藤清[5]討論了隨機(jī)微分方程的推導(dǎo)以及它和確定性微分方程的關(guān)系。但求隨機(jī)微分方程的精確解卻是非常困難的，因此本文將討論求解隨機(jī)微分方程數(shù)值解(近似解)的方法，并應(yīng)用實(shí)例驗(yàn)證比較幾種數(shù)值方法的優(yōu)劣。

先從一階常微分方程解析解的求解出發(fā)，對(duì)比了常微分方程幾種經(jīng)典的數(shù)值求解方法：歐拉方法、改進(jìn)的歐拉方法、三階Runge-Kutta公式和四階Runge-Kutta公式，并舉例編程可視化對(duì)比了精確解(或解析解)和相應(yīng)的數(shù)值解；然后從It隨機(jī)微分方程[5]的產(chǎn)生出發(fā)，舉例圖示對(duì)比了在不同布朗運(yùn)動(dòng)影響下的隨機(jī)常微分方程精確解和相應(yīng)的確定性(deterministic)常微分方程精確解的差別；最后給出了隨機(jī)常微分方程的3種數(shù)值求解方法[6]：Euler-Maruyama方法、Milstein方法和 Runge一Kutta方法，并舉例圖示對(duì)比了不同布朗運(yùn)動(dòng)下的3種數(shù)值方法的求解結(jié)果。

1 常微分方程的數(shù)值方法

先介紹常微分方程的解析解求解，再給出經(jīng)典的數(shù)值求解方法。

1.1 解析解

17世紀(jì)末，牛頓劃時(shí)代的巨著《自然哲學(xué)的數(shù)學(xué)原理》面世了，從此微分方程便誕生了。大多數(shù)數(shù)學(xué)家們一開(kāi)始都努力去尋找微分方程的解析解(精確解)，從而試圖去解釋由微分方程所描繪的物理規(guī)律和自然現(xiàn)象。數(shù)學(xué)家們?cè)诖朔矫嬉泊_實(shí)取得了一系列較成熟的成果，比如一階常微分方程：

(3) 恰當(dāng)微分方程M(x,y)dx+N(x,y)dy=0。

還有二階常系數(shù)線性微分方程：

都得到了較成熟的求解方法。

1.2 數(shù)值解法

但是，到了18世紀(jì)60年代，數(shù)學(xué)家們漸漸意識(shí)到絕大多數(shù)的微分方程是無(wú)法求得它們的解析解的，于是人們才逐漸認(rèn)識(shí)到從另外一個(gè)角度來(lái)研究微分方程的解，即數(shù)值解，是十分有必要的。其中，瑞士數(shù)學(xué)家歐拉在數(shù)值求解方面就做出了開(kāi)創(chuàng)性的工作。1768年，歐拉有關(guān)月球運(yùn)行理論的著作出版了，創(chuàng)造了后來(lái)被廣泛應(yīng)用于求解常微分方程初值問(wèn)題數(shù)值解的方法，即被人們后來(lái)稱道的Euler(歐拉)方法。一階常微分方程的歐拉方法主要有向前的歐拉方法:

yn+1=yn+hf(xn,yn),n=0,1,2,…

向后的歐拉方法:

yn+1=yn+hf(xn+1,yn+1),n=0,1,2,…

改進(jìn)的歐拉方法:

其中，改進(jìn)的歐拉方法在解決實(shí)際問(wèn)題中最為常用。

Euler方法是求解常微分方程初值問(wèn)題的最簡(jiǎn)單數(shù)值方法，但它的收斂階偏低，向前和向后的歐拉方法僅僅一階精度，改進(jìn)的歐拉方法的精度雖然提高到了二階，但對(duì)于更復(fù)雜的實(shí)際問(wèn)題，二階精度卻還是不夠的。于是，在1895年，德國(guó)數(shù)學(xué)家Runge在改進(jìn)的歐拉方法基礎(chǔ)上發(fā)表了《常微分方程數(shù)值解法》，此文章成了高精度數(shù)值求解常微分方程Runge-Kutta方法的開(kāi)端， 在常微分方程數(shù)值方法發(fā)展史上具有里程碑意義。Runge-Kutta方法在創(chuàng)立之時(shí)并未達(dá)到完善，后來(lái)又經(jīng)過(guò)大量數(shù)學(xué)家們多年的共同努力才逐漸完善成熟，目前常用的有三階Runge-Kutta公式：

四階Runge-Kutta公式：

關(guān)于其他的常微分方程數(shù)值解法，限于篇幅有限，就不在此贅述。

1.3 應(yīng)用比較

下面分別用向前的歐拉方法、改進(jìn)的歐拉方法、三階Runge-Kutta公式和四階Runge-Kutta公式數(shù)值求解下列方程：

這是伯努力微分方程，通過(guò)變量變換可轉(zhuǎn)化為非齊次線性微分方程，從而利用常數(shù)變易法可以得到方程的解為

過(guò)點(diǎn)(0，1)的解，即滿足初始條件的精確解為

從圖1的比較可以看出，用向前的歐拉方法( 圖1(a))和改進(jìn)的歐拉方法(圖1(b))求出的數(shù)值解與精確解(The exact solution)的誤差較大，而用三階Runge-Kutta公式(圖1(c) )和四階Runge-Kutta公式(圖1(d))求出的數(shù)值解與精確解的誤差較小，尤其是四階Runge-Kutta公式求出的數(shù)值解與精確解幾乎是完全重合的(圖1(d))，這些數(shù)值結(jié)果與理論分析是完全一致的。

圖1 常微分方程的4種數(shù)值方法

2 隨機(jī)微分方程

2.1 It隨機(jī)微分方程

當(dāng)忽略掉一些比較次要的影響因素，用確定性微分方程[7-8]來(lái)試圖描述自然現(xiàn)象的時(shí)候，實(shí)際上得到的解釋是不完全準(zhǔn)確的。目前，隨著高速計(jì)算機(jī)的產(chǎn)生和對(duì)科學(xué)研究的深入，以及對(duì)自然現(xiàn)象描述解釋準(zhǔn)確度要求的提高，必須重新考慮那些曾經(jīng)被忽略掉的隨機(jī)因素?？茖W(xué)家們也逐漸認(rèn)識(shí)到，隨機(jī)因素不僅僅是對(duì)確定性模型存在缺陷的一個(gè)補(bǔ)充，很多情況下更是反映了物理規(guī)律和自然現(xiàn)象的內(nèi)在本質(zhì)。于是借助于概率論這一工具, 科學(xué)家們將忽略掉的隨機(jī)因素統(tǒng)一建模成隨機(jī)變量，隨機(jī)微分方程的理論研究和數(shù)值研究就勢(shì)在必行了[1]。1951年，日本數(shù)學(xué)家It發(fā)表了影響整個(gè)數(shù)學(xué)界的關(guān)于隨機(jī)微分方程的學(xué)術(shù)論文，自此才有了對(duì)隨機(jī)因素嚴(yán)格意義上的數(shù)學(xué)描述.

引入Wiener過(guò)程的定義[5]。

定義1 如果一個(gè)實(shí)值隨機(jī)過(guò)程W(·)滿足下列3個(gè)條件:

1)W(0)=0;

2)對(duì)所有0≤s≤t，都有W(t)-W(s)服從分布N(0，t-s);

3) 對(duì)所有0

盡管Wiener過(guò)程和實(shí)驗(yàn)相符，但它只是對(duì)真實(shí)布朗運(yùn)動(dòng)的理想化表述，和牛頓力學(xué)相距甚遠(yuǎn)。1930年， Uhlenbeck和Ornsteinv 從牛頓力學(xué)的角度系統(tǒng)地解釋和發(fā)展了布朗運(yùn)動(dòng)的動(dòng)力學(xué)理論。

以y(t) 表示作布朗運(yùn)動(dòng)的粒子在時(shí)刻t的速度，m表示粒子的質(zhì)量。根據(jù)牛頓定律，可以得到如下方程:

dy(t)=-βy(t)dt+σdW(t)

(1)

這里，-β為漂移系數(shù)，σ為擴(kuò)散系數(shù)，W(t)則是Wiener過(guò)程。

若已知初始速度y0，將式(1)寫(xiě)成如下積分方程形式：

當(dāng)σ和-βy(t)均依賴于布朗運(yùn)動(dòng)的粒子速度及時(shí)間t時(shí)，就得到了如下的It隨機(jī)積分方程：

將其寫(xiě)成微分形式：

dy(t)=f(t,y(t))dt+g(t,y(t))dW(t)

(2)

2.2 舉例對(duì)比

dy(t)=-ry(t)dt+σydW(t)

(3)

y(0)=1

其中,r和σ為正常數(shù),在這里取r=σ=1。該隨機(jī)微分方程其滿足初始條件y(0)=1的解析解為

隨機(jī)微分方程式(3)對(duì)應(yīng)的確定性常微分方程為

其滿足初始條件y(0)=1的解析解為

y=e-rt

是關(guān)于時(shí)間t的一條指數(shù)型曲線。

在圖2(a)和2(b)中可以看到，在連續(xù)的布朗運(yùn)動(dòng)(Brownian motion)即隨機(jī)Wiener過(guò)程的影響下，隨機(jī)微分方程(SODE)的精確解在區(qū)間[0,3]上，偏離了相應(yīng)的確定性常微分方程(ODE)的精確解， 受布朗運(yùn)動(dòng)(虛線)影響較大，后面區(qū)間上影響較小。

圖2 It隨機(jī)常微分方程(SODE)和相應(yīng)的確定性常微分方程(ODE)的精確解

3 隨機(jī)微分方程的數(shù)值方法

在這一節(jié)重點(diǎn)介紹幾種經(jīng)典的求解隨機(jī)微分方程的數(shù)值方法。

3.1 幾種數(shù)值方法

yn+1=yn+f(yn)Δt+g(yn)ΔWn

(4)

其中，ΔWn=Wtn+1-Wtn表示定義在區(qū)間[tn,tn+1]上的Wiener過(guò)程的增量，是服從N(0，Δt)分布的相互獨(dú)立的隨機(jī)變量。該歐拉格式的強(qiáng)收斂階只有0.5。

Milstein方法：1974年， Milstein給出了求解隨機(jī)常微分方程式(3)的具有一階強(qiáng)收斂的Milstein方法 ，即

yn+1=yn+f(yn)Δt+g(yn)ΔWn+

Runge-Kutta方法:1982年, Rümelin將Runge-Kutta方法發(fā)展到隨機(jī)微分方程，構(gòu)造了求解式(3)的一階強(qiáng)收斂的隨機(jī)Runge-Kutta方法，其一般格式為

3.2 應(yīng)用比較

下面，舉例給出上述3種方法[6]數(shù)值求解的結(jié)果比較。考慮如下隨機(jī)微分方程：

(5)

該隨機(jī)微分方程式(5)滿足初始條件的精確解為

y=sinWt

在這里，只給出用Euler-Maruyama方法數(shù)值求解隨機(jī)微分方程式(5)的算法，其他兩種可類似得到。

Euler-Maruyama算法：

給定y0,n=0,1,2,…，有

圖3 隨機(jī)微分方程的數(shù)值方法與精確解的比較

在圖3的4個(gè)圖中，虛線均表示連續(xù)的布朗運(yùn)動(dòng)軌跡(Brownian motion)，即隨機(jī)Wiener過(guò)程；實(shí)線均表示隨機(jī)微分方程式(4)的精確解(The exact solution)，余下3條線分別表示用Euler-Maruyama方法、Milstein方法和Runge-Kutta方法得到的數(shù)值解。從圖3可以看到，Milstein方法和Runge-Kutta方法得到的數(shù)值解幾乎完全重合，與精確解的偏差較小，這是因?yàn)镸ilstein方法和Runge-Kutta方法都是一階強(qiáng)收斂的。而從4個(gè)圖整體來(lái)看，Euler-Maruyama方法得到的數(shù)值解偏離精確解較多，這是因?yàn)镋uler-Maruyama方法的強(qiáng)收斂階只有0.5，這些都與理論分析完全吻合。

4 小 結(jié)

先從一階常微分方程解析解的求解出發(fā)，對(duì)比了常微分方程的幾種經(jīng)典數(shù)值求解方法：歐拉方法、改進(jìn)的歐拉方法、三階Runge-Kutta公式和四階Runge-Kutta公式，并舉例編程可視化對(duì)比了精確解和相應(yīng)數(shù)值解的偏差，以此說(shuō)明他們的不同和優(yōu)勢(shì)；其次從It隨機(jī)微分方程的產(chǎn)生出發(fā)，舉例對(duì)比了隨機(jī)常微分方程精確解和相應(yīng)的確定性(deterministic)常微分方程精確解的差別，以此看到隨機(jī)布朗運(yùn)動(dòng)對(duì)微分方程解的影響；最后給出了隨機(jī)常微分方程的3種數(shù)值求解方法：Euler-Maruyama方法、Milstein方法和 Runge-Kutta方法，并舉例圖示對(duì)比了不同布朗運(yùn)動(dòng)下的3種數(shù)值方法的求解結(jié)果以及與精確解的偏差，實(shí)驗(yàn)結(jié)果與理論分析也完全吻合。本文對(duì)隨機(jī)微分方程的學(xué)習(xí)和數(shù)值解的求解及應(yīng)用都有一定的指導(dǎo)意義。