随机微分方程数值实验 平衡法 Euler方法(matlab)

题目:

计算带有乘性噪声的一维Ito型积分,即

$\begin{aligned} & d X(t)=\sigma X(t) d W(t) \\ & X(0)=x .\end{aligned}$

其真解为

$X(t)=x e^{-\frac{1}{2} \sigma^2 t+\sigma W(t)}$

一类平衡法数值格式为

$X_{k+1}=X_k+\sigma X_k \Delta_k W(h)+\sigma\left(X_k-X_{k+1}\right)\left|\Delta_k W(h)\right|$

Euler方法格式为

$\begin{aligned} & X_{K+1}=X_k+\sigma X_k \Delta_k \omega(h) \\ & X_0=x .\end{aligned}$

我们计算当$\sigma=4, x=1$时,比较两种方法的数值解逼近真解的情况。

程序:

randn('state',100)
sigma = 4;T = 1;Xzeros = 1;
N = 1.5e+5;
dt = 1/N;
dW = sqrt(dt)*randn(1,N);%Brownian increments
W = cumsum(dW);
%----------------true solution
Xtrue = Xzeros*exp(sigma*W-sigma^2/2*[dt:dt:T]);
R = 5;
Dt = R*dt;
L = N/R;
XE = zeros(1,L);
Xtemp = Xzeros;
%----------------Euler explicit method
for j = 1:L
    Winc = sum(dW(R*(j-1)+1:R*j));
    Xtemp = Xtemp+sigma*Xtemp*Winc;
    XE(j) = Xtemp;
end
%---------------The balanced method
XB = zeros(1,L);
Xtemp = Xzeros;
for j = 1:L
    Winc = sum(dW(R*(j-1)+1:R*j));
    Xtemp = ((1+sigma*abs(Winc))*Xtemp+sigma*Xtemp*Winc)/(1+sigma*abs(Winc));
    XB(j) = Xtemp;
end
figure
plot([(0:Dt:T)],[Xzeros,XE],'r-*');
hold on;
plot([0:dt:T],[Xzeros Xtrue],'k--');
legend('显式Euler数值解','真解');
figure
plot([(0:Dt:T)],[Xzeros,XB],'y+'),
hold on
plot([0:dt:T],[Xzeros Xtrue],'k--');
legend('平衡法数值解','真解');
erro=max(abs([Xzeros,XE]-[Xzeros,XB]))

结果:

平衡方法和Euler方法的数值解图像和真解图像基本吻合,通过这两种方法得到数值解相差0.0830。

  • 13
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

菜鸟粥粥

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值