随机系统基础

第九讲 随机系统与随机微分方程 Desmond J. Higham†, An Algorithmic Introduction to Numerical Simulation of StochasticDesmond J. Higham†, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM REVIEW 2001, Vol. 43,No . 3,pp . 525Differential Equations, SIAM REVIEW 2001, Vol. 43,No . 3,pp . 525 – –546546 C. W. Gardiner, Handbook of Stochastic s, Springer, 1996, (C. W. Gardiner, Handbook of Stochastic s, Springer, 1996, ( 中文版，中文版， 上海科技出版社上海科技出版社) ) 9.19.1 布朗运动布朗运动 十九世纪末之前，人们都认为微分方程求解解释自然现象，只要收集数据， 则我们可以预测未来。但是，量子力学基于纯粹的统计原理，二是混沌的出现发 现微小的误差可以被迅速放大，这使得人们对于预测未来的信心开始怀疑。 1827 年 Rober Brown 研究了悬浮在液体表面的颗粒无规则运动以证实分子的 存在，这种运动的数学解释是 1905 年 Einstein 给出的。 假设n个粒子悬浮在液体中，时间间隔远小于观测时间内，单个粒子x坐标增 加一个量，位移在和  d之间的粒子数目为dn，满足 dn  n()d 这里()就是概率密度，满足    ()d 1，()() 设f (x,t)为单位体积内粒子数目，那么经过时间间隔，时间t 通过点x和 x  dx之间的粒子数目为  f (x,t )dx f (x  ,t)()ddx  按照小参数展开得到 ff2f f  f()d ()d  2txx 12 令()d  D，上式化为 2 2  2 ()d  f (x,t)2f (x,t)  D tx2 这就是著名的扩散方程，偏微分方程的一种。其解为 n x2/(4Dt)f (x,t) e 4Dt 可以看出，沿 x 轴的方差为 2 x2 2Dt 是随着时间增加而增加的。 法国物理学家Langevin 提出了认为比爱因斯坦更加简单的办法描述布朗运动  x   f (t)m x 受力分成二部分，一為粘滯阻力-αv，一為隨机作用力f (t)（实际就是噪声） ，粘 滯阻力仍來自介質分子對顆粒的碰撞，將顆粒看作半徑為 a 的小球，在粘滯系 數為 η 的流體中運動，則有  6a 上式稱為斯托克斯公式(stokes’s law)。將上式對大量顆粒子求平均，即把大群顆 的運動方程相加然後用顆粒數，Langevin 推导出  x2 2kT t  这 里k是 波 尔 兹 曼 常 数 ， 而 由 扩 散 方 程 (diffusion equation) 已 经 得 到 2 x2 2Dt，那么比较可得 D  kT  这是和实际相符合的。正是由于对于布朗运动的研究，诞生了第一个随机微分方 程  x   f (t)m x 经过 Ornstein, Langevin, 和 Winer 1918 年严格的数学描述， 此后 Brown 运动 专指 Winer 过程。 标准的 Winer 过程是在[0,T]时间间隔内的连续随机变量W(t)，满足下面三 个条件： A：W(0)  0，以概率 1 B：对于时刻0  s  t  T，变量的增量W(t)-W(s)为正态分布，均值为 0， 方差为t s，即W(t)-W(s)~t sN(0,1)， C 对于时刻0  s  t  u  v  T，增量W(t)-W(s)和增量W(v)-W(u)是相互独 立的，又称为独立增量过程。 从计算的角度， 我们把 Winer 过程离散化， 采样时间t  T / N， 这里整数 N， 我们记t j  jt，那么W(t j ) W j ，按照上面条件有： A：W0 0 B：W j W j1 dW j ，j 1,2,,N CdW j ~tN(0,1)，且不同j下标下相互独立 例子：Matlab 仿真布朗运动 %BPATH1 Brownian path simulation randn( state ,100) % set the state of randn T = 1; N = 50; dt = T/N; dW = zeros(1,N); % preallocate arrays 预先分配内存 W = zeros(1,N); % for efficiency dW(1) = sqrt(dt)*randn; % first approximation outside the loop . W(1) = dW(1); % since W(0) = 0 is not allowed for j = 2:N dW(j) = sqrt(dt)*randn; % general increment W(j) = W(j-1) + dW(j); end plot([0:dt:T],[0,W], r- ) % plot W against t xlabel( t , FontSize ,16) ylabel( W(t) , FontSize ,16, Rotation ,0) ------------------ Listing 1 M-file bpath1.m. 0.2 0 -0.2 -0.4 W(t) -0.6 -0.8 -1 -1.2 -1.40 0.10.20.30.40.50.60.70.80.91 t Fig.1 布朗运动的模拟 解释：不同种子能够得到不同的运动轨迹，由于 Matlab 向量起始下标为 1，所以 W0 0应该表示为W(1)，但是画的过程中将其标记为时间 0。bpath1.m.这个程序 效率不高，因为没有利用 Matlab 向量化的特点，将其优化，我们得到程序 %BPATH2 Brownian path simulation: vectorized randn( state ,100) % set the state of randn T = 1; N = 50; dt = T/N; dW = sqrt(dt)*randn(1,N); % increments W = cumsum(dW); % cumulative sum plot([0:dt:T],[0,W], r- ) % plot W against t xlabel( t , FontSize ,16) ylabel(’W(t)’,’Fon