%EM Euler-Maruyama method on linear SDE % % This can be modified to simulate the equation by the Milstein method by % simply changing one line of the code where the update rule over a time % step is defined % SDE is dX = lambda*X dt + mu*X dW, X(0) = Xzero, % where lambda = 2, mu = 1 and Xzero = 1. % clear * clf randn('state',100); % seed the random number generator so results are % reproducible from one run to the next. % Change this if one wants different % results (choices of random variables) % with each run. lambda = 2; mu = 1; Xzero = 1; % problem parameters T = 1; N = 2^9; % time interval, number of time steps M=100; % number of realizations(experiments) to execute dt = T/N; % time step XTsum=0; % This is where we will store the running sum of the solution X(T) % at the ending time, which will be used % to estimate the mean value of X(T) % from the simulations for s=1:M % loop over realizations (experiments) 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); % exact solution for comparison % Xem keeps track of the numerically simulated solution % Xtemp is the numerically simulated value for the solution at the current % time step Xem = zeros(1,N); % preallocate for efficiency Xtemp = Xzero; % start with the initial data for j = 1:N Winc = dW(j); % the Brownian increment over a time step Xtemp = Xtemp + dt*lambda*Xtemp + mu*Xtemp*Winc; % Numerical update: Euler-Marayama % To do Milstein simulation, simply replace this equation with the % formula for the Milstein update rule Xem(j) = Xtemp; end XTsum=XTsum+Xem(N); % Add the value X(T) from the current realization (experiment) % to the running sum over all realizations end % These plot the results of the last realization (experiment) plot(0:dt:T,[Xzero,Xtrue],'m-'), hold on plot(0:dt:T,[Xzero,Xem],'r--*'), hold off xlabel('t','FontSize',12) ylabel('X','FontSize',16,'Rotation',0,'HorizontalAlignment','right') % This is to calcluate the error between the numerical and exact solution % in the last realization. emerr = abs(Xem(end)-Xtrue(end)); disp(sprintf('The numerical error is %f at the end time T=%f.',emerr,T)) figure(1) XTmean=XTsum/M; % mean value of X(T) as estimated by numerical simulations XTmeantrue=Xzero*exp(lambda*T); % the theoretically correct average value for the mean value of X(T). disp(sprintf('The numerical estimate for the mean of X(T) is %f, while the exact answer is %f.',... XTmean,XTmeantrue))