随机微分方程数值解法.ppt

《随机微分方程数值解法.ppt》由会员分享，可在线阅读，更多相关《随机微分方程数值解法.ppt（37页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、随机微分方程数值解随机微分方程数值解法法现在学习的是第1页，共37页随机微分方程数值解法1.随机微分方程概述随机微分方程概述1.1 布朗运动介绍布朗运动介绍1.2 随机积分随机积分1.3 两种形式的随机微分方程两种形式的随机微分方程2.随机微分方程数值方法介绍随机微分方程数值方法介绍2.1 随机随机Taylor展开展开2.1 Euler方法方法2.2 Milstein方法方法3.数值试验数值试验3.1 精度数值试验精度数值试验3.2 稳定性数值试验稳定性数值试验 现在学习的是第2页，共37页1.随机微分方程概述随机微分方程概述 布布朗朗运运动动是是历历史史上上最最早早被被认认真真研研究究过过的

2、的随随机机过过程程。1827 年年，英英国国生生物物学学家家布布朗朗(Robert Brown)首首先先观观察察和和研研究究了了悬悬浮浮在在液液体体中中的的细细小小花花粉粉微微粒粒受受到到水水分分子子连连续续撞撞击击形形成成的的运运动动情情况况，布布朗朗运运动动也也因因此此而而得得名名。1905 年年爱爱因因斯斯坦坦(Einstein)对对它它做做出出了了合合理理的的物物理理解解释释并并求求出出了了微微粒粒的的转转移移密密度度，1918 年年维维纳纳(Norbert Wiener)在在数数学学上上严严格格地地定定义义了了布布朗朗运运动动(因因此此它它有有时时也也称称为为维维纳纳过过程程)。现现

3、在在布布朗朗运运动动已已经经成成为为了了描描述述随随机机现象的基石。现象的基石。1.1 布朗运动介绍布朗运动介绍现在学习的是第3页，共37页 物理上理解，布朗运动的起因是液体的所有分子都处在运动物理上理解，布朗运动的起因是液体的所有分子都处在运动中，而且相互碰撞，从而微粒周围有大量的分子以微小但起伏不中，而且相互碰撞，从而微粒周围有大量的分子以微小但起伏不定的力共同作用于它，使它被迫作不规则运动。如果用定的力共同作用于它，使它被迫作不规则运动。如果用 表示表示微粒在时刻微粒在时刻 所处位置的一个坐标，由于液体是均匀的，自然设所处位置的一个坐标，由于液体是均匀的，自然设想从时刻想从时刻 到到 的

4、位移的位移 是许多几乎完全独立的小位移是许多几乎完全独立的小位移之和，因而根据中心极限定理，可以合理的假定之和，因而根据中心极限定理，可以合理的假定 服从服从正态分布，而且对于不同时间段的位移应该是相互独立的。因此正态分布，而且对于不同时间段的位移应该是相互独立的。因此，布朗运动有如下定义：，布朗运动有如下定义：现在学习的是第4页，共37页 定义定义1.1 一个随机过程一个随机过程 ,它在一个微小时间间隔它在一个微小时间间隔 之间内的变化为之间内的变化为 。如果。如果1);2)，其中，其中 为一常数。为一常数。3)对于任何两个不同时间间隔对于任何两个不同时间间隔,的值相互独立的值相互独立,即独

5、立增量。即独立增量。称随机变量称随机变量 的运动遵循的运动遵循(标准标准)维纳过程或者布朗运维纳过程或者布朗运动。动。若若 ，则称，则称 为为标准布朗运动标准布朗运动或或标准标准Wiener过程。过程。现在学习的是第5页，共37页注：注：1）布朗运动是处处连续的，并且它是处处是不可微的。直观）布朗运动是处处连续的，并且它是处处是不可微的。直观上来看，这意味着它的运动轨迹相当曲折。上来看，这意味着它的运动轨迹相当曲折。2）对于标准布朗运动，）对于标准布朗运动，即，即若记随机变量若记随机变量 则有则有 形式上看，当形式上看，当 时，如同普通微积分中的情形，有：时，如同普通微积分中的情形，有：由于布

6、朗运动是处处不可微的，此处的由于布朗运动是处处不可微的，此处的 只能视为一种简单记只能视为一种简单记法。法。现在学习的是第6页，共37页布朗运动的模拟布朗运动的模拟 以下对一维的布朗运动进行随机模拟。一维的布朗运动可以以下对一维的布朗运动进行随机模拟。一维的布朗运动可以看做质点在直线上作简单随机游动，则看做质点在直线上作简单随机游动，则 表示质点在时刻表示质点在时刻 时时在直线上的位置。利用在直线上的位置。利用Matlab模拟布朗运动的程序代码如下：模拟布朗运动的程序代码如下：%布朗运动的模拟布朗运动的模拟 randn(state,100)%设置随机数发生器的状态设置随机数发生器的状态 T=1

7、；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不不 允许允许 for j=2:N dW(j)=sqrt(dt)*randn；现在学习的是第7页，共37页 W(j)=W(j-1)+dW(j);endplot(0:dt:T,0,W,r-)%绘图绘图 xlabel(t,FontSize,16)ylabel(W(t),FontSiz

8、e,16,Rotation,0)现在学习的是第8页，共37页图图1 布朗运动布朗运动现在学习的是第9页，共37页 还可以如下进行模拟：还可以如下进行模拟：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)现在学习的是

9、第10页，共37页111.2 随机积分随机积分 随机积分分为随机积分分为It型随机积分型随机积分和和Stratonovich型随机积分型随机积分。以。以下假设下假设Wiener过程过程 定义在概率空间定义在概率空间 上，上，为为 的上升滤子（即的上升滤子（即 且对且对 )，对任意，对任意 ，关于关于 可测，且满足可测，且满足 此外，对随机过程此外，对随机过程 引入以下三个条件：引入以下三个条件：现在学习的是第11页，共37页 关于关于 可测；可测；(1)即即 为为 可测的；可测的；(2)(3)以下是以下是It型随机积分的定义：型随机积分的定义：定义定义1.2 设设 为标准布朗运动，随机过程为标

10、准布朗运动，随机过程 满足条件满足条件(1)-(3)。对。对 ，将，将 作划分，任取作划分，任取 令令 若随机变量序列若随机变量序列 现在学习的是第12页，共37页 (4)均方收敛于唯一极限，则称均方收敛于唯一极限，则称 (5)为为 关于关于 在在 上的上的It积分。上述定积分。上述定义中，作和式（义中，作和式（4）时不能像通常积分那样，）时不能像通常积分那样，在在 中任取中任取，否则可能导致均方极限不存在。（，否则可能导致均方极限不存在。（5）中取的是的）中取的是的 的的左左端点端点 ，得到，得到It型随机积分型随机积分 。现在学习的是第13页，共37页 若取区间若取区间 的中点的中点 时，

11、就得到时，就得到Stratonovich型积分型积分，记为，记为 。现在学习的是第14页，共37页1.3 两种形式的随机微分方程两种形式的随机微分方程 随机微分方程亦分为随机微分方程亦分为It型随机微分方程型随机微分方程和和Stratonovich型型随机微分方程随机微分方程。目前研究的较多的。目前研究的较多的It型随机微分方程的一般形型随机微分方程的一般形式如下：式如下：(6)其中其中 均为均为 上的上的Borel可测函数，分别被称为漂移系数和扩散可测函数，分别被称为漂移系数和扩散系数。系数。现在学习的是第15页，共37页方程方程(6)的积分形式为：的积分形式为：(7)其中的随机积分为其中的

12、随机积分为It型随机积分。型随机积分。若将若将It型随机积分替换为型随机积分替换为Stratonovich型随机积分，则型随机积分，则(7)式式变为变为 (8)对应的微分方程为对应的微分方程为 (9)现在学习的是第16页，共37页方程（方程（9）即为）即为Stratonovich型随机微分方程。型随机微分方程。注：注：1)It型随机微分方程型随机微分方程(6)和和Stratonovich型随机微分方程（型随机微分方程（9）是可以相互转换的。在标量情形下，对方程（是可以相互转换的。在标量情形下，对方程（6）令）令 在矢量情形下，令在矢量情形下，令 其中其中 则方程（则方程（6）可以转化为）可以转

13、化为Stratonovich性随机微分性随机微分方程如下：方程如下：现在学习的是第17页，共37页注：注：1)大部分随机微分方程的解析解是无法获得的，可以求得解大部分随机微分方程的解析解是无法获得的，可以求得解析解的随机微分方程多为线性随机微分方程。析解的随机微分方程多为线性随机微分方程。2)有些随机微分方程的解析解虽然可以求到，但是形式很复有些随机微分方程的解析解虽然可以求到，但是形式很复杂，处理起来很不方便。杂，处理起来很不方便。3)在实际应用中，实用的方法是在计算机上进行数值求解，在实际应用中，实用的方法是在计算机上进行数值求解，即不直接求出即不直接求出 的解析解，而是在解所存在的区间上

14、，求得一的解析解，而是在解所存在的区间上，求得一系列点系列点 上的近似值。上的近似值。现在学习的是第18页，共37页2.随机微分方程数值方法介绍随机微分方程数值方法介绍 目前随机微分方程的数值求解方法有目前随机微分方程的数值求解方法有Euler方法、方法、Milstein方方法法、Runge-Kutta方法等。方法等。Runge-Kutta方法的复杂程度比方法的复杂程度比Euler方法和方法和Milstein方法的程度要高。在实际应用中，一般情况下用方法的程度要高。在实际应用中，一般情况下用Euler方法和方法和Milstein方法来对模型进行数值模拟。由于方法来对模型进行数值模拟。由于It型

15、随机型随机微分方程与微分方程与Stratonovich型随机微分方程是可以相互相互转化的，型随机微分方程是可以相互相互转化的，以下介绍求解以下介绍求解It型随机微分方程型随机微分方程(6)的的Euler方法和方法和Milstein方法。方法。首先给出随机微分方程解的首先给出随机微分方程解的存在唯一性定理存在唯一性定理以及数值方法以及数值方法强强收敛收敛与与弱收敛弱收敛的定义如下：的定义如下：现在学习的是第19页，共37页 定理定理2.1(解的存在唯一性定理）解的存在唯一性定理）若若 满足满足(i)(线性增长条件）存在正常数线性增长条件）存在正常数 使得使得 (ii)(Lipschitz条件）条

16、件）存在正常数存在正常数 使得使得 且有且有 ，则方程则方程(6)存在唯一解且存在唯一解且 。现在学习的是第20页，共37页定义定义 2.1(强收敛性强收敛性)若存在常数若存在常数 (与与 独立独立)，使得，使得 则称该数值方法是则称该数值方法是 阶强收敛的。阶强收敛的。定义定义2.2(弱收敛性）弱收敛性）若对适当的若对适当的 次可微的多项式次可微的多项式 ，存，存在在 ，使得：，使得：则称该数值方法是则称该数值方法是 阶弱收敛的。阶弱收敛的。现在学习的是第21页，共37页 强收敛性与弱收敛性是数值方法的两种收敛性评价标准。强强收敛性与弱收敛性是数值方法的两种收敛性评价标准。强收敛性要求对随机

17、微分方程进行数值模拟时，数值近似的轨迹必收敛性要求对随机微分方程进行数值模拟时，数值近似的轨迹必须充分接近真实轨迹。弱收敛则并不关注解过程的轨迹，而仅仅须充分接近真实轨迹。弱收敛则并不关注解过程的轨迹，而仅仅是解过程的矩性质。是解过程的矩性质。现在学习的是第22页，共37页2.1 随机随机Taylor展开展开 方便起见，对如下的标量自治型随机微分方程进行讨论：方便起见，对如下的标量自治型随机微分方程进行讨论：(10)其中其中 是标准是标准Wiener过过程。程。随机随机Taylor展开式是随机微分方程数值算法的基础，展开式是随机微分方程数值算法的基础，Euler算算法和法和Milstein算法

18、都是在随机算法都是在随机Taylor展开式不同的地方截断而得到展开式不同的地方截断而得到的数值算法。的数值算法。设设 是正整数，是正整数，利用随机利用随机Taylor展开式和展开式和It公式，可以得到：公式，可以得到：现在学习的是第23页，共37页 其中其中 是余项，算子是余项，算子 和和 分别为分别为 则则(10)式可以写为：式可以写为：现在学习的是第24页，共37页 (12)求解方程求解方程(10)的的Euler方法和方法和Milstein方法均是在方法均是在(12)的基础上进行的基础上进行截断而得到的。截断而得到的。现在学习的是第25页，共37页2.2 Euler 方法方法对于方程对于方

19、程(9)，Euler方法的格式如下：方法的格式如下：(13)注：注：1）Euler 方法的强收敛阶是方法的强收敛阶是 ，弱收敛阶是，弱收敛阶是 1.2)方法方法(13)为显式的为显式的Euler方法，还有如下形式的半隐式方法，还有如下形式的半隐式Euler方法和半隐式方法和半隐式Euler方法：方法：半隐式半隐式Euler方法方法:隐式隐式Euler方法：方法：现在学习的是第26页，共37页2.3 Milstein方法方法对于方程对于方程(10)，Milstein方法的格式如下：方法的格式如下：(14)注：注：1）Milsten 方法的强收敛阶是方法的强收敛阶是 1.2)方法方法(14)为显式

20、的为显式的Milsten方法，还有如下形式的半隐式方法，还有如下形式的半隐式Milstein方法和半隐式方法和半隐式Milsten方法：方法：半隐式半隐式Milsten方法方法:现在学习的是第27页，共37页隐式隐式Milstein方法：方法：注：注：Euler方法和方法和Milstein方法方法 的形式比较简单，是求解随机微分的形式比较简单，是求解随机微分方程最常用的两种数值方法。方程最常用的两种数值方法。现在学习的是第28页，共37页3.数值试验数值试验3.1 精度数值试验精度数值试验 精度即误差，即用数值方法求出的数值解和精确解之间的差精度即误差，即用数值方法求出的数值解和精确解之间的差

21、异。对于可以求出解析解的随机微分方程，可以通过比较数值解异。对于可以求出解析解的随机微分方程，可以通过比较数值解和精确解之间轨迹的差异，也可以通过比较平均绝对误差来比较和精确解之间轨迹的差异，也可以通过比较平均绝对误差来比较。若用数值方法求解随机微分方程时，进行。若用数值方法求解随机微分方程时，进行 次样本模拟，记次样本模拟，记 和和 表示第表示第 次模拟时在点次模拟时在点 处的数值解和精确解处的数值解和精确解,则：则：即为即为平均绝对误差。平均绝对误差。现在学习的是第29页，共37页 以下对以下对Euler方法进行精度数值试验。选取线性试验方程如下方法进行精度数值试验。选取线性试验方程如下

22、(15)方程（方程（15)的解析解为：的解析解为：令令 现在学习的是第30页，共37页程序程序1(数值解与精确解轨迹比较）数值解与精确解轨迹比较）randn(state,100)lambda=2;mu=1;Xzero=1;%参数赋值参数赋值T=1;N=28;dt=1/N;dW=sqrt(dt)*randn(1,N);%布朗增量，用于模拟数值解布朗增量，用于模拟数值解W=cumsum(dW);%累加求和，用于模拟精确解累加求和，用于模拟精确解 Xtrue=Xzero*exp(lambda-0.5*mu2)*(dt:dt:T)+mu*W);%求精确解求精确解plot(0:dt:T,Xzero,Xt

23、rue,m-),hold on%绘出精确解轨迹绘出精确解轨迹现在学习的是第31页，共37页R=4;Dt=R*dt;L=N/R;%设置数值求解的步长，改变设置数值求解的步长，改变R可以改变可以改变 DtXem=zeros(1,L);%预分配，提高效率预分配，提高效率 Xtemp=Xzero;for j=1:L Winc=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-*),hold off%绘制数值解轨迹绘制数

24、值解轨迹xlabel(t,FontSize,12)ylabel(X,FontSize,16,Rotation,0,HorizontalAlignment,right)现在学习的是第32页，共37页图图 2 Euler方法数值解与精确解的轨迹比较方法数值解与精确解的轨迹比较现在学习的是第33页，共37页平均绝对误差的比较：平均绝对误差的比较：randn(state,100)lambda=-2;mu=1;Xzero=1;T=1;N=211;dt=T/N;M=1000;Xerr=zeros(M,1);for s=1:M%进行进行M次循环，即进行次循环，即进行M次样本模拟次样本模拟 dW=sqrt(d

25、t)*randn(1,N);W=cumsum(dW);Xtrue=Xzero*exp(lambda*T+mu*W(end);现在学习的是第34页，共37页 Xtemp=Xzero;for j=1:N Xtemp=Xtemp*(1+lambda*dt+mu*dW(j);end Xerr(s,1)=abs(Xtemp-Xtrue);%计算每一次模拟在端点处的绝对误差计算每一次模拟在端点处的绝对误差endmean(Xerr)%计算平均绝对误差计算平均绝对误差 现在学习的是第35页，共37页不同步长求解的平均绝对误差比较：不同步长求解的平均绝对误差比较：0.16300.11850.10470.0960现在学习的是第36页，共37页3.2 稳定性数值试验稳定性数值试验 通常用均方稳定性来衡量数值算法的稳定性。假设运用数值通常用均方稳定性来衡量数值算法的稳定性。假设运用数值算法求解随机微分方程算法求解随机微分方程(15)后得到的迭代格式为：后得到的迭代格式为：其中其中 ，为步长为步长。称。称 为均方稳定函数；若为均方稳定函数；若 ，则称数值方法关于，则称数值方法关于 是均方稳定的。若令是均方稳定的。若令 记记 为为 ，记，记为为 ，则该数值方法的均方稳定区域为，则该数值方法的均方稳定区域为 现在学习的是第37页，共37页