随机微分方程数值实验 Euler-Maruyama方法（matlab）

求解随机微分方程如下：

微分方程的真解为：

Euler-Maruyama方法的数值格式：

程序如下：

randn('state',50)
lambda=1.6;mu=2.3;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],'g-'),hold on 
 R=1;Dt=R*dt;L=N/R;
 Xem=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;
 end
 plot([0:Dt:T],[Xzero,Xem],'r--*'),hold off
 xlabel('t','FontSize',20)
 ylabel('X','FontSize',20,'Rotation',0,'HorizontalAlignment','right')
 title('R=1时数值模拟示意图')
 emerr=abs(Xem(end)-Xtrue(end))

% Mean-square and asymptotic stability test for E-M

randn('state',150)
T=15;M=40000;Xzero=1;
ltype={'g-','r--','b-.'};
subplot(211)
lambda=1.6;mu=2.3;
for k=1:3
    Dt=2^(1-k);
    N=T/Dt;
    Xms=zeros(1,N);Xtemp=Xzero*ones(M,1);
    for j=1:N
        Winc=sqrt(Dt)*randn(M,1);
        Xtemp=Xtemp+Dt*lambda*Xtemp+mu*Xtemp.*Winc;
        Xms(j)=mean(Xtemp.^2);
    end
    semilogy([0:Dt:T],[Xzero,Xms],ltype{k},'Linewidth',2),hold on
end
legend('\Delta t=1','\Delta t=1/2','\Delta t=1/4')
title('Mean-Square:\lambda=1.6,\mu=2.3','FontSize',16)
ylabel('E[x^2]','FontSize',12),axis([0,T,1e-20,1e+20]),hold off

randn('state',130)
lambda=1.6;mu=2.3;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],'g-'),hold on 
 R=4;Dt=R*dt;L=N/R;
 Xem=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;
 end
 plot([0:Dt:T],[Xzero,Xem],'r--*'),hold off
 xlabel('t','FontSize',16)
 ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
 title('R=4,randn(state,130)时数值模拟示意图')
 emerr=abs(Xem(end)-Xtrue(end))

结果：