隨機常微分方程的幾種數(shù)值求解方法及其應用*

李 煥 榮

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

0 引 言

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

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

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

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

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

1.1 解析解

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

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

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

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

1.2 數(shù)值解法

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

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

向后的歐拉方法:

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

改進的歐拉方法:

其中，改進的歐拉方法在解決實際問題中最為常用。

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

四階Runge-Kutta公式：

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

1.3 應用比較

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

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

過點(0，1)的解，即滿足初始條件的精確解為

從圖1的比較可以看出，用向前的歐拉方法( 圖1(a))和改進的歐拉方法(圖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 隨機微分方程

2.1 It隨機微分方程

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

引入Wiener過程的定義[5]。

定義1 如果一個實值隨機過程W(·)滿足下列3個條件:

1)W(0)=0;

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

3) 對所有0

盡管Wiener過程和實驗相符，但它只是對真實布朗運動的理想化表述，和牛頓力學相距甚遠。1930年， Uhlenbeck和Ornsteinv 從牛頓力學的角度系統(tǒng)地解釋和發(fā)展了布朗運動的動力學理論。

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

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

(1)

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

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

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

將其寫成微分形式：

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

(2)

2.2 舉例對比

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

(3)

y(0)=1

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

隨機微分方程式(3)對應的確定性常微分方程為

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

y=e-rt

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

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

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

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

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

3.1 幾種數(shù)值方法

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

(4)

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

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

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

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

3.2 應用比較

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

(5)

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

y=sinWt

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

Euler-Maruyama算法：

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

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

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

4 小 結(jié)

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