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

隨機(jī)微分方程數(shù)值解法隨機(jī)微分方程數(shù)值解法1.隨機(jī)微分方程概述1.1布朗運(yùn)動介紹1.2隨機(jī)積分1.3兩種形式的隨機(jī)微分方程2.隨機(jī)微分方程數(shù)值方法介紹2.1隨機(jī)Taylor展開2.1Euler方法2.2Milstein方法3.數(shù)值試驗(yàn)3.1精度數(shù)值試驗(yàn)3.2穩(wěn)定性數(shù)值試驗(yàn)

第2頁,共37頁，2024年2月25日，星期天1.隨機(jī)微分方程概述

布朗運(yùn)動是歷史上最早被認(rèn)真研究過的隨機(jī)過程。1827年，英國生物學(xué)家布朗(RobertBrown)首先觀察和研究了懸浮在液體中的細(xì)小花粉微粒受到水分子連續(xù)撞擊形成的運(yùn)動情況，布朗運(yùn)動也因此而得名。1905年愛因斯坦(Einstein)對它做出了合理的物理解釋并求出了微粒的轉(zhuǎn)移密度，1918年維納(NorbertWiener)在數(shù)學(xué)上嚴(yán)格地定義了布朗運(yùn)動(因此它有時也稱為維納過程)?，F(xiàn)在布朗運(yùn)動已經(jīng)成為了描述隨機(jī)現(xiàn)象的基石。

1.1布朗運(yùn)動介紹第3頁,共37頁，2024年2月25日，星期天

物理上理解，布朗運(yùn)動的起因是液體的所有分子都處在運(yùn)動中，而且相互碰撞，從而微粒周圍有大量的分子以微小但起伏不定的力共同作用于它，使它被迫作不規(guī)則運(yùn)動。如果用表示微粒在時刻所處位置的一個坐標(biāo)，由于液體是均勻的，自然設(shè)想從時刻到的位移是許多幾乎完全獨(dú)立的小位移之和，因而根據(jù)中心極限定理，可以合理的假定服從正態(tài)分布，而且對于不同時間段的位移應(yīng)該是相互獨(dú)立的。因此，布朗運(yùn)動有如下定義：第4頁,共37頁，2024年2月25日，星期天

定義1.1一個隨機(jī)過程,它在一個微小時間間隔之間內(nèi)的變化為。如果1);2)，其中為一常數(shù)。3)對于任何兩個不同時間間隔,的值相互獨(dú)立,即獨(dú)立增量。稱隨機(jī)變量的運(yùn)動遵循(標(biāo)準(zhǔn))維納過程或者布朗運(yùn)動。若，則稱為標(biāo)準(zhǔn)布朗運(yùn)動或標(biāo)準(zhǔn)Wiener過程。第5頁,共37頁，2024年2月25日，星期天注：1）布朗運(yùn)動是處處連續(xù)的，并且它是處處是不可微的。直觀上來看，這意味著它的運(yùn)動軌跡相當(dāng)曲折。2）對于標(biāo)準(zhǔn)布朗運(yùn)動，，即若記隨機(jī)變量則有形式上看，當(dāng)時，如同普通微積分中的情形，有：

由于布朗運(yùn)動是處處不可微的，此處的只能視為一種簡單記法。第6頁,共37頁，2024年2月25日，星期天布朗運(yùn)動的模擬

以下對一維的布朗運(yùn)動進(jìn)行隨機(jī)模擬。一維的布朗運(yùn)動可以看做質(zhì)點(diǎn)在直線上作簡單隨機(jī)游動，則表示質(zhì)點(diǎn)在時刻時在直線上的位置。利用Matlab模擬布朗運(yùn)動的程序代碼如下：

%布朗運(yùn)動的模擬randn('state',100)%設(shè)置隨機(jī)數(shù)發(fā)生器的狀態(tài)T=1；N=500；dt=T/N；dW=zeros(1,N);%布朗增量存放位置W=zeros(1,N)；%預(yù)分配，提高效率dW(1)=sqrt(dt)*randn；%循環(huán)前的初始化W(1)=dW(1)；%Matlab中數(shù)組下標(biāo)從1開始，故W(0)=0不允許

forj=2:NdW(j)=sqrt(dt)*randn；第7頁,共37頁，2024年2月25日，星期天

W(j)=W(j-1)+dW(j);endplot([0:dt:T],[0,W],’r-’)

%繪圖xlabel(’t’,’FontSize’,16)ylabel(’W(t)’,’FontSize’,16,’Rotation’,0)第8頁,共37頁，2024年2月25日，星期天圖1布朗運(yùn)動第9頁,共37頁，2024年2月25日，星期天還可以如下進(jìn)行模擬：

randn('state',100)T=1；N=500；dt=T/N；

dW=sqrt(dt)*randn(1,N)；%向量化,提高運(yùn)算效率

W=cumsum(dW);%累加和計算命令，W(j)=dW(1)+dW(2)+…+dW(j)；j=1,…N

plot([0:dt:T],[0,W],’r-’)

%繪圖

xlabel(’t’,’FontSize’,16)

ylabel(’W(t)’,’FontSize’,16,’Rotation’,0)

第10頁,共37頁，2024年2月25日，星期天111.2隨機(jī)積分

隨機(jī)積分分為Itó型隨機(jī)積分和Stratonovich型隨機(jī)積分。以下假設(shè)Wiener過程定義在概率空間上，

為的上升濾子（即且對)，對任意，關(guān)于可測，且滿足

此外，對隨機(jī)過程引入以下三個條件：第11頁,共37頁，2024年2月25日，星期天

關(guān)于可測；(1)即為可測的；(2)(3)

以下是Itó型隨機(jī)積分的定義：

定義1.2設(shè)為標(biāo)準(zhǔn)布朗運(yùn)動，隨機(jī)過程滿足條件(1)-(3)。對，將作劃分，任取

令若隨機(jī)變量序列

第12頁,共37頁，2024年2月25日，星期天

(4)均方收斂于唯一極限，則稱

(5)為關(guān)于在上的Itó積分。上述定義中，作和式（4）時不能像通常積分那樣，在中任取，否則可能導(dǎo)致均方極限不存在。（5）中取的是的的左端點(diǎn)，得到Itó型隨機(jī)積分。

第13頁,共37頁，2024年2月25日，星期天

若取區(qū)間的中點(diǎn)時，就得到Stratonovich型積分，記為。

第14頁,共37頁，2024年2月25日，星期天1.3兩種形式的隨機(jī)微分方程

隨機(jī)微分方程亦分為Itó型隨機(jī)微分方程和Stratonovich型隨機(jī)微分方程。目前研究的較多的Itó型隨機(jī)微分方程的一般形式如下：(6)其中

均為上的Borel可測函數(shù)，分別被稱為漂移系數(shù)和擴(kuò)散系數(shù)。第15頁,共37頁，2024年2月25日，星期天方程(6)的積分形式為：(7)

其中的隨機(jī)積分為Itó型隨機(jī)積分。若將Itó型隨機(jī)積分替換為Stratonovich型隨機(jī)積分，則(7)式變?yōu)?8)對應(yīng)的微分方程為(9)

第16頁,共37頁，2024年2月25日，星期天方程（9）即為Stratonovich型隨機(jī)微分方程。注：1)Itó型隨機(jī)微分方程(6)和Stratonovich型隨機(jī)微分方程（9）是可以相互轉(zhuǎn)換的。在標(biāo)量情形下，對方程（6）令

在矢量情形下，令

其中則方程（6）可以轉(zhuǎn)化為Stratonovich性隨機(jī)微分方程如下：

第17頁,共37頁，2024年2月25日，星期天注：1)

大部分隨機(jī)微分方程的解析解是無法獲得的，可以求得解析解的隨機(jī)微分方程多為線性隨機(jī)微分方程。2)有些隨機(jī)微分方程的解析解雖然可以求到，但是形式很復(fù)雜，處理起來很不方便。3)在實(shí)際應(yīng)用中，實(shí)用的方法是在計算機(jī)上進(jìn)行數(shù)值求解，即不直接求出的解析解，而是在解所存在的區(qū)間上，求得一系列點(diǎn)上的近似值。第18頁,共37頁，2024年2月25日，星期天2.隨機(jī)微分方程數(shù)值方法介紹

目前隨機(jī)微分方程的數(shù)值求解方法有Euler方法、Milstein方法、Runge-Kutta方法等。Runge-Kutta方法的復(fù)雜程度比Euler方法和Milstein方法的程度要高。在實(shí)際應(yīng)用中，一般情況下用Euler方法和Milstein方法來對模型進(jìn)行數(shù)值模擬。由于Itó型隨機(jī)微分方程與Stratonovich型隨機(jī)微分方程是可以相互相互轉(zhuǎn)化的，以下介紹求解Itó型隨機(jī)微分方程(6)的Euler方法和Milstein方法。首先給出隨機(jī)微分方程解的存在唯一性定理以及數(shù)值方法強(qiáng)收斂與弱收斂的定義如下：

第19頁,共37頁，2024年2月25日，星期天

定理2.1(解的存在唯一性定理）若滿足(i)(線性增長條件）存在正常數(shù)使得

(ii)(Lipschitz條件）存在正常數(shù)使得

且有，則方程(6)存在唯一解且。第20頁,共37頁，2024年2月25日，星期天定義2.1(強(qiáng)收斂性)若存在常數(shù)(與獨(dú)立)，，使得

則稱該數(shù)值方法是階強(qiáng)收斂的。定義2.2(弱收斂性）若對適當(dāng)?shù)拇慰晌⒌亩囗検?，存在，使得?/p>

則稱該數(shù)值方法是階弱收斂的。第21頁,共37頁，2024年2月25日，星期天

強(qiáng)收斂性與弱收斂性是數(shù)值方法的兩種收斂性評價標(biāo)準(zhǔn)。強(qiáng)收斂性要求對隨機(jī)微分方程進(jìn)行數(shù)值模擬時，數(shù)值近似的軌跡必須充分接近真實(shí)軌跡。弱收斂則并不關(guān)注解過程的軌跡，而僅僅是解過程的矩性質(zhì)。

第22頁,共37頁，2024年2月25日，星期天2.1隨機(jī)Taylor展開

方便起見，對如下的標(biāo)量自治型隨機(jī)微分方程進(jìn)行討論：(10)

其中是標(biāo)準(zhǔn)Wiener過程。

隨機(jī)Taylor展開式是隨機(jī)微分方程數(shù)值算法的基礎(chǔ)，Euler算法和Milstein算法都是在隨機(jī)Taylor展開式不同的地方截斷而得到的數(shù)值算法。

設(shè)是正整數(shù)，利用隨機(jī)Taylor展開式和Itó公式，可以得到：第23頁,共37頁，2024年2月25日，星期天

其中是余項，算子和分別為

則(10)式可以寫為：

第24頁,共37頁，2024年2月25日，星期天

(12)求解方程(10)的Euler方法和Milstein方法均是在(12)的基礎(chǔ)上進(jìn)行截斷而得到的。第25頁,共37頁，2024年2月25日，星期天2.2Euler方法對于方程(9)，Euler方法的格式如下：

(13)注：1）

Euler方法的強(qiáng)收斂階是，弱收斂階是1.2)方法(13)為顯式的Euler方法，還有如下形式的半隱式Euler方法和半隱式Euler方法：半隱式Euler方法:

隱式Euler方法：

第26頁,共37頁，2024年2月25日，星期天2.3Milstein方法對于方程(10)，Milstein方法的格式如下：

(14)注：1）Milsten方法的強(qiáng)收斂階是1.2)方法(14)為顯式的Milsten方法，還有如下形式的半隱式Milstein方法和半隱式Milsten方法：半隱式Milsten方法:

第27頁,共37頁，2024年2月25日，星期天隱式Milstein方法：

注：Euler方法和Milstein方法的形式比較簡單，是求解隨機(jī)微分方程最常用的兩種數(shù)值方法。第28頁,共37頁，2024年2月25日，星期天3.數(shù)值試驗(yàn)3.1精度數(shù)值試驗(yàn)

精度即誤差，即用數(shù)值方法求出的數(shù)值解和精確解之間的差異。對于可以求出解析解的隨機(jī)微分方程，可以通過比較數(shù)值解和精確解之間軌跡的差異，也可以通過比較平均絕對誤差來比較。若用數(shù)值方法求解隨機(jī)微分方程時，進(jìn)行次樣本模擬，記和表示第次模擬時在點(diǎn)處的數(shù)值解和精確解,則：即為平均絕對誤差。第29頁,共37頁，2024年2月25日，星期天

以下對Euler方法進(jìn)行精度數(shù)值試驗(yàn)。選取線性試驗(yàn)方程如下

(15)方程（15)的解析解為：

令第30頁,共37頁，2024年2月25日，星期天程序1(數(shù)值解與精確解軌跡比較）randn('state',100)lambda=2;mu=1;Xzero=1;%參數(shù)賦值T=1;N=2^8;dt=1/N;dW=sqrt(dt)*randn(1,N);%布朗增量，用于模擬數(shù)值解W=cumsum(dW);%累加求和，用于模擬精確解Xtrue=Xzero*exp((lambda-0.5*mu^2)*([dt:dt:T])+mu*W);

%求精確解plot([0:dt:T],[Xzero,Xtrue],'m-'),holdon%繪出精確解軌跡第31頁,共37頁，2024年2月25日，星期天R=4;Dt=R*dt;L=N/R;

%設(shè)置數(shù)值求解的步長，改變R可以改變DtXem=zeros(1,L);%預(yù)分配，提高效率

Xtemp=Xzero;forj=1:LWinc=sum(dW(R*(j-1)+1:R*j));

%計算布朗增量

Xtemp=Xtemp+Dt*lambda*Xtemp+mu*Xtemp*Winc;Xem(j)=Xtemp;endplot([0:Dt:T],[Xzero,Xem],'r--*'),holdoff

%繪制數(shù)值解軌跡xlabel('t','FontSize',12)ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right')第32頁,共37頁，2024年2月25日，星期天圖2Euler方法數(shù)值解與精確解的軌跡比較第33頁,共37頁，2024年2月25日，星期天平均絕對誤差的比較：