The Time Evolution of a Particle’s Position

To determine the position of a particle can be a very challenging task. The obvious way is to look at it and see where it is at that time. However, as one can expect, measuring the position of a particle at every instance of time is very difficult (impossible). (A particle could have a deterministic trajectory, for example in space without any collisions, but imagine dust particles in the frequently colliding with each other.)

To overcome this obstacle, we can turn to probability theory and try to model the particle’s position in a probabilistic manner. Let’s say that we measure the position of the particle at time 0, and we want to guess where it will be at time 1, and then at time 2. For simplicity, assume that the particle can either go up, or down with equal chances. With this assumption, it is clear that at time 0, it will either be at +1 or -1. At time 2, things get a bit more complicated, because the particle can be at 3 different locations and at time 3 there are 4 different possibilities (see graph below).

As we can see, looking at it from a probabilistic point of view, we can make very significant inferences. For example, at time 2 we can deduce that the probability of the particle being at position 0 is 50%, while on the edges is 25%-25%.

Brownian motion

Now, let’s assume that the particle follows a more realistic random walk, the arithmetic Brownian motion.

Here, X denotes the position, mu denotes the drift, sigma the variance and W is the Weiner Process. The key element of this is the Weiner process, which is a standard normally distributed random variable N(0,1). In layman’s terms, at every timestep, a number is drawn from the Weiner process, which is then scaled accordingly by sigma , which will be the change in the particle’s position over that timestep.

So, if the particle follows the above Brownian motion, we can not draw out the possible positions because there are infinitely many. We can, however, look at it from the perspective of the distributions. Distributions are the generalisation of the draw out method above, so we can draw the distribution of the particle’s position at every timestep.

Fokker-Planck equation

As a result, we have transformed our problem (Where is the particle in a 1000 timesteps?) into evaluating its distribution at every timestep. But finding the distribution is analytically is also challenging, so we can use our computers to numerically simulate the particle’s position accordingly. We turn to the Fokker-Planck equation (FPE)to simulate the probability density function of this particle. The FPE describes the time evolution of a pdf, which is exactly what we need (The FPE was developed precisely to solve this, not a coincidence).

The FPE is a partial differential equation, and the Brownian motion satisfies it. Solving the FPE is technically an approximation, but for our purpose it works really well.

Numerically solving the FPE

To solve this partial differential equation numerically, we need to set up a grid. To do this, we discretize the timesteps and the possible positions by small steps. On the x-axis we have time, on the y-axis, we have the position, that is x and on the z-axis, we will have p(t,x) , which will be the probability density function(scroll down to see the pictures).

This grid allows us to simulate the PDE after specifying the boundary conditions. The top and bottom boundary conditions are that p(t,x_max) = p(t,x_min) = 0 , which means that the probability of the particle is at x_max, x_min is 0. This means that the particle will never reach that point, so we need to choose x_max, x_min accordingly (usually 8*sigma away is enough). We assume that we know the initial position of the particle, which means that at the point (0,0), the probability will be one, and everywhere else its 0:p(0,0)=1 . Note that at each time, the probabilities sum to 1.

Furthermore, we need a way to approximate the derivatives. In our model, muand sigma are constants, so we can take them out as multipliers. We approximate the derivative using the central difference method because of its low computational cost and high reliability.

On the bottom line, we now have a direct expression for the change in the probability density function with respect to time. Rearranging and introducing a,b,c:

We have reduced our PDE to a system of ODEs, which we can solve. Here, a,band c are constants, so we can rewrite the equation in matrix form.

Lambda denotes a tridiagonal matrix andp_t, q_t are vectors. Applying the Euler-Maruyama scheme with a backwards difference approximation, we can get a direct expression for p_t :

Using this expression we can finally simulate the time evolution of the particle following a Brownian motion.

Code and results

import numpy as np
import scipy.sparse
import matplotlib.pyplot as plt
import fmlib
from scipy.stats import norm
from scipy.sparse.linalg import inv
from mpl_toolkits.mplot3d import axes3d# This function calculates a,b,c.
def compute_abc(sigma, mu, dx):
a = ((sigma**2)/(2*(dx**2)) + (mu)/(2*dx))
b = -((sigma**2)/dx**2)
c = ((sigma**2)/(2*(dx**2)) - (mu)/(2*dx))
return a, b, c# This function calculates the tridiagonal matrix Lambda.
def compute_lambda( a,b,c,l ):
a_col = np.ones(l)*a
b_col = np.ones(l)*b
c_col = np.ones(l)*c
return scipy.sparse.diags( [a_col[1:],b,c_col[:-1]],offsets= [-1,0,1])

With the above helper functions, we can simulate p_t

# simulates the FPE
def simulate_fokker_planc(X0, T, mu, sigma, N, M):
dt = T/N
X_max = X0 + 8*(sigma)
X_min = X0 - 8*(sigma)dx = (X_max - X_min)/M
X = np.linspace(X_min, X_max, M+1)
t = np.linspace(0,T,N+1)P = np.zeros((N+1, M+1))
# Initial and boundary conditions
P[0,:] = np.zeros(M+1) #first column
P[:,0] = np.zeros(N+1) #bottom row
P[:,-1] = np.zeros(N+1) #top rowP[0,int(M/2+X0)] = 1a,b,c = compute_abc(sigma, mu, dx)
lmbda = compute_lambda(a,b,c,N-1)
identity = scipy.sparse.identity(N-1)for i in range(0,N):P[i+1,1:M] = (inv(identity - dt*lmbda)).dot( P[i,1:M] )return P.T, t, XP,t,X = simulate_fokker_planc(0,1,1,2,80,80)# this is to plot a 3D graph.a
t_mesh, X_mesh = np.meshgrid(t,X)
ax = plt.gca(projection='3d')
ax.plot_surface(X_mesh, t_mesh, P alpha=0.7,edgecolor='k',cmap='viridis')
ax.set_title("FPE for mu = 1, sigma = 2")

The above code generates a 3D plot of the position of the particle moving forward in time.

As we can see, the particle starts off from a fixed point and then evolves in time. Taking cross-sections allows us to see how the distribution evolves from timestep to timestep.

The distribution flattens out as time goes on, which indicates that we know less about where the particle is. It actually approaches a normal distribution, because the driving process of the random walk is a Wiener process. Another interpretation could be to think about the Fokker-Planck equation as ink diffusing in water.

Finishing comments

The Fokker-Planck equation is a very powerful tool in mathematical modelling. It is also known as the Kolmogorov forward equation and it can come up in many areas of science, such as statistical mechanics, quantum mechanics, chemistry, geography or even in mathematical finance.

I hope my article was useful. Please let me know if you have any questions, I’m happy to explain /discuss.