%EMWEAK Test weak convergence of Euler-Maruyama % % Solves dX = lambda*X dt + mu*X dW, X(0) = Xzero, % where lambda = 2, mu = 1 and Xzer0 = 1. % % E-M uses different timesteps: 2^(p-10), p = 1,2,3,4,Nres % Examine weak convergence at T=1: | E (X_L) - E (X(T)) |. % % Different paths are used for each E-M timestep. % Code is vectorized over paths. % % Uncommenting the line indicated below gives the weak E-M method. clear * randn('state',100); lambda = 2; mu = 0.1; % problem parameters Xzero = 1; % initial condition T = 1; % simulation time N= 2^10; % number of time steps at finest resolution M = 50000; % number of paths sampled Nres=5; % number of resolution levels dt=T/N; % time step at finest resolution Xem = zeros(Nres,1); % preallocate arrays for p = 1:Nres % take various Euler timesteps Dt = 2^p*dt; L = T/Dt; % L Euler steps of size Dt Xtemp = Xzero*ones(M,1); for j = 1:L Winc = sqrt(Dt)*randn(M,1); % Brownian increment. No need to use same realization % of randomness in different % simulations because only weak % error is examined %Winc = sqrt(Dt)*sign(rand(M,1)-0.5); %% use for weak E-M %% Xtemp = Xtemp + Dt*lambda*Xtemp + mu*Xtemp.*Winc; % Euler-Marayama update end Xem(p) = mean(Xtemp); % mean of X(T) sampled from the simulations at resolution level p end Xerr = abs(Xem - exp(lambda*T)); % error from exact solution: X(T) should have mean exp(lambda*T) Dtvals = 2.^(1:Nres)*dt; % Time step resolution levels subplot(222) % top RH picture loglog(Dtvals,Xerr,'b*-'), hold on loglog(Dtvals,Dtvals,'r--'), hold off % reference slope of 1 axis([1e-3 1e-1 1e-4 1]) xlabel('\Delta t'), ylabel('| E(X(T)) - Sample average of X_L |') title('Euler-Marayama Weak Error','FontSize',10) %%%% Least squares fit of error = C * dt^q %%%% A = [ones(Nres,1), log(Dtvals)']; rhs = log(Xerr); sol = A\rhs; q = sol(2); resid = norm(A*sol - rhs); figure(1) disp(sprintf('The least squares fit for the slope of the weak E-M error in the log-log plot is %f, with residual %f.',q,resid))