%EM Euler-Maruyama method on linear SDE
%
% SDE is dX = lambda*X dt + mu*X dW, X(0) = Xzero,
% where lambda = 2, mu = 1 and Xzero = 1.
%
% Discretized Brownian path over [0,1] has dt = 2^(-8).
% Euler-Maruyama uses timestep R*dt.
clear all;close all;
randn('state',100)
lambda = 0.5; %2
mu = 1; Xzero = 0.0001; % problem parameters
T = 1; N = 2^8; dt = T/N;
dW = sqrt(dt)*randn(1,N); % Brownian increments
W = cumsum(dW); % discretized Brownian path
%每个元素对应的列向上(行向左)求和矩阵
Xtrue = Xzero*exp((lambda-0.5*mu^2)*([dt:dt:T])+mu*W);
% plot([0:dt:T],[Xzero,Xtrue],'m-'), hold on
R = 4; Dt = R*dt; L = N/R; % L EM steps of size Dt = R*dt
Xem = zeros(1,L); % preallocate for efficiency
Xtemp = Xzero;
for j = 1:L
Winc = sum(dW(R*(j-1)+1:R*j));
Xtemp = Xtemp + Dt*lambda*Xtemp + mu*sqrt(Xtemp)*Winc...
+ 0.25*(Winc^2 - Dt); %近似值的表达 *sqrt(Xtemp)*0.5*(Xtemp)^(-0.5)*
Xem(j) = Xtemp;
end
plot([0:Dt:T],[Xzero,Xem],'r--*'), hold off
xlabel('t','FontSize',12)
ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
grid on
set(gca,'xgrid','on','gridcolor','bl')
emerr = abs(Xem(end)-Xtrue(end));
a=sum(Xem<0);
评论2