根据彩色树理论，构造了两种求解Stratonovich型随机微分方程的半隐式三阶随机Runge-Kutta方法，给出了这两种方法的稳定性分析，其稳定区域比现有方法的稳定区域大；数值模拟的结果表明两个方法都具有较高的精度。
朱晓临,王子洁,李井刚:两种半隐式三阶随机 Runge-Kutta方法 2012,48(15) 局部截断误差,,表示随机 runge-Kuta方法在tn的 4 666 p=(B B2 B,) 整休截断误差。对n=0,1,…,N,如果 1/2 取f=(164616),由式(6)和式(7)解得: n)=o(b)且Ec)-O0 b=-1,c=2或a=1,b 43 10 35 7 那么E--o(n)。 相应地,得到下列两个半隐式三阶随机 Runge 由定理1知,若数值方法满足均方条件: utt格式,并分别称之为D1和D2方法。 12 Elym-y D1方法: 均值条件:Eyn-(切=0(n2.则随机 Runge-Kutta 2=,+4(0)+(x+号42() 方法1阶(p=1)强收敛到随机微分方程(1)。 根据彩色树理论,由上述均方条件,树t,12t有 E(h-hae)=0,E(-Je)|=0 x=y,+“(x)+4()+( (x)+48(2)+g( 0 D2方法: 化简得: ae=l,B'e=l, p Be=2 6 Y2=y+[()+(x)+11 由上述均值条件,树41,120,t2有: Y1=y,+5(x)+(+ El(J2-h, Be a'e+la B g()-439(x) E(Jo-hBAe pe+(e ae h y=+(+4()+( JiB(Be)= (x)+4()+( 3f(be)15 4 h+ ah(Be) 4稳定性分析 EUJm-JiB B(Be) 对线性检验随机微分方程 dy= ay dt by o dw(t) 36-36(Be)l5h2+15h(/B(Be) 应用随机 Runge-Kutta方法,设得到的迭代格式为 化简上面四个方程,得 R(h, a,b,D) B Ae 其中J=△wn=w(tn+)-w(n)N0,b),h为步长 SBe) pB(Be)=I (7)b∈Cl 定义3称R(,a,b)=E(b为数值方法 在上述过程中用到° 的均方稳定函数,如果对给定步长h,R(h,a,b)<1 E(2+)=0,E(2) (2k 2k! 则称此方法是均方稳定的。 取 Butcher点阵为 若令p=mh,q=bh,则记RVh,a,b)为R(p,q, 000 称S={(D,引)|(p,9)<l为该数值方法的稳定区 44 域,简称稳定域。 0 定理2对于方程(5),D1方法的均方稳定域为 (P,9)|R(p,9)<l},即 000 (x)+(y2)+37)1672)+27y3+67274 <1(9) 则 ,B=a00 8-36D+6p 48-36p+6p 0 其中 48+12p Y2=12pq-2pq+ 2 2012,48(15) Computer Engineering and Applications计算机工程与应用 y3=24q+8pq,4=8q 对式(12)求均方,得: 证明将D1方法应用于方程(8),得: E(R(P, q) (x)2+()2+37) R(p,q)=1+kh(1+4R2+1R3)+ 680-1260+20y2) 41(+4R2+R3) (10) 16(74)+23+62<1 其中J 4+p+ 680-1260p+210p N(0,1),R 4 于是D2方法的均方稳定的条件为 R=1+5(R2+R3)+(2R2-1 16(4)+2y3+672 + <1 8+21p2+8wh+8y/+8() (1680-1260+210m2)(1680-1260+2l02) (4-p)(2-p) 下面将D1、D2方法与比较常见的两种半隐式三 令1=48+12-612-2p3,y2=12-22q+48,阶随机 Runge-Kutta方法S1M1方法和Ts1方法的稳 y3=24g2+8p2,y4=8q,对式(10)求均方,得 定区域进行比较,结果见图1。 ()+(y2)+3(y3 3.0 DI D2 48-36p+6 DI D2 SIMI l6(y)+23+672y 1.5 48-36+6p 于是D1方法的均方稳定的条件为 0.5 ()+(2)+3(72)16()+23+6y -10-9-8-7-6-5-4-3-2-10 (48-362+57)(48-367+67) 图1D1D2、SIⅥ1和TS1方法的稳定 定理3对于方程(8),D2方法的均方稳定城为 SIM1方法: S={(D,y)R(p,)<1} Y2=yn+|/( (X)-f(x)+18() + 680-1260p+210p y=,+(m)+(y)+2( 6(2)+273+672y (11) 号小(x)+g() 680-1260p+210p yx;=y,+4/(x)+3/(x+ 其中,=1680+420p-210p2-70p3,72=420pq 41:(x)+38(2) 178p2q+1680q,73=8402+172,4=2W SM1方法的均方稳定域为 证明将D2方法应用于方程(8),有: p2+4+6)894+723-640n2-624p+1152 R(P,q)=1+h(+4R2+1R3)+ 4(3-p) 643-p)(1-p) q 1(+4R2+R3 (12) 213p-264p+1008 20+5p+14q 64(3-p)(1-p 其中R2=-20-5 2(3-p)(1-p) TS1方法 R=1+2(8+)+1(982-3) 280+70+35y2+284xy150971+2807(4 =+4(0++0 (4-p)(70-35p) (4-p)(70-35p y=n+b/(x)+2(x2)-8(1 令n1=1680+420D-210p2-702,72=420g n1=y+(x)+4(x)+(x)+ 178p2q+1680g,3=8402+172py2,4=280q3 g(r)+4g(Y2)+g(rs)] 朱晓临,王子洁,李井刚:两种半隐式三阶随机 Runge-Kutta方法 2012,48(15) 仁+精俯解 1.7 D方 1.7 1.6 上n2方法 14 1 65432 l111 543 2 1.1 1.0b 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 00.10.20.30.40.50.60.70.80.91.0 00.10.20.30.40.50.60.70.80.91.0 00.10.20.30.40.50.60.70.80.91.0 图2精确解的一次路径 图3D1方法模拟的一次路径 图4D2方法数模拟的一次路径 1.7 SlM方法 765432 →TSI方法 11111110 65432 l1ll11 8 11000 098 701020304050607080910 7 00.10.20.30.40.50.60.70.80.91.0 图5SIMl方法模拟的一次路径 图6TS1方法模拟的一次路径 IS1方法的均方稳定域为: 表1D1、D2算法的最大误差和平均误差 ()+(2)+3(2)6+2m2+624∠1 DI D2 h 最大误差平均误差最大误差平均误差 (24-6p)2 24-6p) 0.1135 0.0493 0.1139 0.0495 其屮%=24+18+6p2+p2,%2=18p9+6p2+24g,3= 0.0721 0.0441 0.0723 0.0442 1292+9p2,4=2q。 表2TS1SIM1算法的最大误差和平均误差 从图1可以看出,D1、D2方法的稳定区域比SM1、 TSI SIMI h TS1方法的稳定区域都要大。通过分析上述四种方 最大误差平均误差最大误差平均误差 0.11380.04940.11300.049 法,D1、D2方法的稳定冈域之所以比SIM1、TS1方法 0.07220.04410.07200.0440 的稳定区域都大,主要是 Butcher点阵和a,B的选择 更合理,使得Y1,Y2Y3在整个D1、D2方法中能更均是精确解式(14)利四种数值方法在500次样本模拟 衡地发挥作用使得稳定区域更大。另外与TS1方法其中的一次路径图,显示四种数值方法都具自较高的 相比,D1、D2方法中的Y是隐式的,这也是Dl、D2精度。 方法的稳定区域比TS1方法的稳定区域大的另一个 原因。 6结论 根据彩色树理论,构造了两种半隐式三阶随机 5数值实验 Runge-Kutta方法:D1方法和D2方法,并分析了它们 在实验方程(8)中,取=-2,=1,方程化为 的稳定性,其稳定区域比现有的半隐式三阶随机 ( dy(=-2y(t)dt+(Oodw(), tE[O (13) Runge-Kut方法的稳定区域都大,数值模拟的结果 表明本文的两个方法都具有较高的精度。 上式的解析解为 y()=exp-2.5+W(t) (14) 参考文献: 表1、表2给出对于不同步长h,分别采用500次 [1] Kloeden P E, Platen E Numerical solution of stochastic 样木模拟的四种算法的最大误差和平均误差 differential equations[M]. Berlin: Springer-Verlag, 1995 由表1和表2可知本文构造的两个三级随机Ru 161-221 ge-Kutta方法与SM1TS1方法的精度一样。图2~6 (下转128页)