If you were to look at a window with sunlight shining in, you might be able to notice little particles of dust flying around. No matter how hard you try to swat the particles of dust they will continue to float around aimlessly. This phenomenon is described by Brownian motion. Trying to predict where these particles go can be extremely difficult due to the randomness of their motion. In this file exchange entry, Emma shows us how we can use the Euler?Maruyama (EM) method to simulate how the particles move in a 1 dimensional (1D) plane.
When asked about what drove Emma to write this live script on such a topic, she said that the “project was created from [her] research on Brownian particle motion while assisting Prof. Atzberger’s (UCSB) work on differentiable manifolds.” But then she adapted the work into this live script to compete in the MATLAB Online Live Editor Challenge hosted in 2018.
Main code applying EM method
Now let’s dive into some of the core elements of this live script. In this live script we can see how to apply the EM method for 1D Brownian motion but also the approximate error with this analysis. We can see this plotted below with the blue line representing the actual path and the green dotted line the predicted motion.
lambda = 2;
mu = 1;
Xnot = 1;
T = 1;
N = 2^8;
dt = 1/N;
dW = sqrt(dt)*randn(1,N);
W = cumsum(dW);
X_exact = Xnot*exp((lambda-.5*mu^2)*([dt:dt:T])+mu*W);
R = 4;
Dt = R*dt;
L = N/R;
X_EM = zeros(1,L);
X_temp = Xnot;
for j = 1:L
Winc = sum(dW(R*(j-1)+1:R*j));
X_temp = X_temp + Dt*lambda*X_temp + mu*X_temp*Winc;
X_EM(j) = X_temp;
title('Approximation error with 10^4 samples')
Solving Brownian simulation using other methods
This code also has a bonus as it shows how we would also solve this problem numerically. This directly shows the difference in the approach of tackling the problem but also shows how different it is coding these two methods.
numx = 101;
numt = 2000;
dx = 1/(numx - 1);
dt = 0.00005;
x = 0:dx:1;
C = zeros(numx,numt);
t(1) = 0;
C(1,1) = 0;
C(1,numx) = 0;
mu = 0.5;
sigma = 0.05;
C(i,1) = exp(-(x(i)-mu)^2/(2*sigma^2)) / sqrt(2*pi*sigma^2);
t(j+1) = t(j) + dt;
C(i,j+1) = C(i,j) + (dt/dx^2)*(C(i+1,j) - 2*C(i,j) + C(i-1,j));
Benefits of Live scripts
In this live script not only do we see the code in order to do this analysis, but Emma also takes us through the theory as well (pictured below). I believe this to be one of the most powerful elements in this live script. In the live script Emma shows how the EM method is derived to apply to Brownian motion. Which allows for better understanding of the code. When asked about using a live script for this project she said that,
“using live scripts is a great way to compile functions and to share them.
It makes it very easy for others to learn from since you can see the code and
run it in real time. This was actually my first time experimenting with the
live script function in MATLAB for the competition and I think it is definitely
something that I would use in the future.”
Let us know what you think here leave a comment for Emma.