隨機微分方程數值解法第一頁，共三十七頁，編輯于2023年，星期二隨機微分方程數值解法1.隨機微分方程概述1.1布朗運動介紹1.2隨機積分1.3兩種形式的隨機微分方程2.隨機微分方程數值方法介紹2.1隨機Taylor展開2.1Euler方法2.2Milstein方法3.數值試驗3.1精度數值試驗3.2穩定性數值試驗

第二頁，共三十七頁，編輯于2023年，星期二1.隨機微分方程概述

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

1.1布朗運動介紹第三頁，共三十七頁，編輯于2023年，星期二

物理上理解，布朗運動的起因是液體的所有分子都處在運動中，而且相互碰撞，從而微粒周圍有大量的分子以微小但起伏不定的力共同作用于它，使它被迫作不規則運動。如果用表示微粒在時刻所處位置的一個坐標，由于液體是均勻的，自然設想從時刻到的位移是許多幾乎完全獨立的小位移之和，因而根據中心極限定理，可以合理的假定服從正態分布，而且對于不同時間段的位移應該是相互獨立的。因此，布朗運動有如下定義：第四頁，共三十七頁，編輯于2023年，星期二

定義1.1一個隨機過程,它在一個微小時間間隔之間內的變化為。如果1);2)，其中為一常數。3)對于任何兩個不同時間間隔,的值相互獨立,即獨立增量。稱隨機變量的運動遵循(標準)維納過程或者布朗運動。若，則稱為標準布朗運動或標準Wiener過程。第五頁，共三十七頁，編輯于2023年，星期二注：1）布朗運動是處處連續的，并且它是處處是不可微的。直觀上來看，這意味著它的運動軌跡相當曲折。2）對于標準布朗運動，，即若記隨機變量則有形式上看，當時，如同普通微積分中的情形，有：

由于布朗運動是處處不可微的，此處的只能視為一種簡單記法。第六頁，共三十七頁，編輯于2023年，星期二布朗運動的模擬

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

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

forj=2:NdW(j)=sqrt(dt)*randn；第七頁，共三十七頁，編輯于2023年，星期二

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)第八頁，共三十七頁，編輯于2023年，星期二圖1布朗運動第九頁，共三十七頁，編輯于2023年，星期二還可以如下進行模擬：

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

dW=sqrt(dt)*randn(1,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)

第十頁，共三十七頁，編輯于2023年，星期二111.2隨機積分

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

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

此外，對隨機過程引入以下三個條件：第十一頁，共三十七頁，編輯于2023年，星期二

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

以下是Itó型隨機積分的定義：

定義1.2設為標準布朗運動，隨機過程滿足條件(1)-(3)。對，將作劃分，任取

令若隨機變量序列

第十二頁，共三十七頁，編輯于2023年，星期二

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

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

第十三頁，共三十七頁，編輯于2023年，星期二

若取區間的中點時，就得到Stratonovich型積分，記為。

第十四頁，共三十七頁，編輯于2023年，星期二1.3兩種形式的隨機微分方程

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

均為上的Borel可測函數，分別被稱為漂移系數和擴散系數。第十五頁，共三十七頁，編輯于2023年，星期二方程(6)的積分形式為：(7)

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

第十六頁，共三十七頁，編輯于2023年，星期二方程（9）即為Stratonovich型隨機微分方程。注：1)Itó型隨機微分方程(6)和Stratonovich型隨機微分方程（9）是可以相互轉換的。在標量情形下，對方程（6）令

在矢量情形下，令

其中則方程（6）可以轉化為Stratonovich性隨機微分方程如下：

第十七頁，共三十七頁，編輯于2023年，星期二注：1)

大部分隨機微分方程的解析解是無法獲得的，可以求得解析解的隨機微分方程多為線性隨機微分方程。2)有些隨機微分方程的解析解雖然可以求到，但是形式很復雜，處理起來很不方便。3)在實際應用中，實用的方法是在計算機上進行數值求解，即不直接求出的解析解，而是在解所存在的區間上，求得一系列點上的近似值。第十八頁，共三十七頁，編輯于2023年，星期二2.隨機微分方程數值方法介紹

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

第十九頁，共三十七頁，編輯于2023年，星期二

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

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

且有，則方程(6)存在唯一解且。第二十頁，共三十七頁，編輯于2023年，星期二定義2.1(強收斂性)若存在常數(與獨立)，，使得

則稱該數值方法是階強收斂的。定義2.2(弱收斂性）若對適當的次可微的多項式，存在，使得：

則稱該數值方法是階弱收斂的。第二十一頁，共三十七頁，編輯于2023年，星期二

強收斂性與弱收斂性是數值方法的兩種收斂性評價標準。強收斂性要求對隨機微分方程進行數值模擬時，數值近似的軌跡必須充分接近真實軌跡。弱收斂則并不關注解過程的軌跡，而僅僅是解過程的矩性質。

第二十二頁，共三十七頁，編輯于2023年，星期二2.1隨機Taylor展開

方便起見，對如下的標量自治型隨機微分方程進行討論：(10)

其中是標準Wiener過程。

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

設是正整數，利用隨機Taylor展開式和Itó公式，可以得到：第二十三頁，共三十七頁，編輯于2023年，星期二

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

則(10)式可以寫為：

第二十四頁，共三十七頁，編輯于2023年，星期二

(12)求解方程(10)的Euler方法和Milstein方法均是在(12)的基礎上進行截斷而得到的。第二十五頁，共三十七頁，編輯于2023年，星期二2.2Euler方法對于方程(9)，Euler方法的格式如下：

(13)注：1）

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

隱式Euler方法：

第二十六頁，共三十七頁，編輯于2023年，星期二2.3Milstein方法對于方程(10)，Milstein方法的格式如下：

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

第二十七頁，共三十七頁，編輯于2023年，星期二隱式Milstein方法：

注：Euler方法和Milstein方法的形式比較簡單，是求解隨機微分方程最常用的兩種數值方法。第二十八頁，共三十七頁，編輯于2023年，星期二3.數值試驗3.1精度數值試驗

精度即誤差，即用數值方法求出的數值解和精確解之間的差異。對于可以求出解析解的隨機微分方程，可以通過比較數值解和精確解之間軌跡的差異，也可以通過比較平均絕對誤差來比較。若用數值方法求解隨機微分方程時，進行次樣本模擬，記和表示第次模擬時在點處的數值解和精確解,則：即為平均絕對誤差。第二十九頁，共三十七頁，編輯于2023年，星期二

以下對Euler方法進行精度數值試驗。選取線性試驗方程如下

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

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

%求精確解plot([0:dt:T],[Xzero,Xtrue],'m-'),holdon%繪出精確解軌跡第三十一頁，共三十七頁，編輯于2023年，星期二R=4;Dt=R*dt;L=N/R;

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

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

%繪制數值解軌跡xlabel('t','FontSize',12)ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right')第三十二頁，共三十七頁，編輯于2023年，星期二圖2Euler方法數值解與精確解的軌跡比較第三十三頁，共三十七頁，編輯于2023年，星期二平均絕對誤差的比較：