Diffusion equation

The diffusion equation is widely discussed in the literature, especially in physics and mathematics. However, the diffusion equation also has relevance in other fields, for instance, biology, where it can be useful to understand some biological phenomena at different scales, including population biology. Therefore, it is crucial to know how to solve it.

The diffusion equation is a partial differential equation which involves both a time and space. In 1D it would be:

$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$

It can be discretised using numerical methods. In this case, we are going to use the finite difference method, so that the above equation can be discretised as follows:

$\frac{u^{t+1}_{i} - u^{t+1}_{i}}{\Delta t} = D \frac{u^t_{i+1} - 2u^t_{i} + u^t_{i-1}}{\Delta^2x}$

where, $t$ indicate the time set and $i$ the x position.

To obtain the steady-state, we are going to programming this in python. The below code solve the diffusion equation numerically

import numpy as np
import matplotlib.pyplot as plt

#Setting parameters
nx = 101
dx = 0.01
nt = 20000
D = .1
dt = 0.0001 #be aware setting this value

#Initial conditions
u = np.zeros(nx)
u[25:75] = 2
un = u.copy()

#Numerical scheme: Finite differences
for n in range(nt):
    for i in range(1, nx-1):
        un = u.copy()
        u[i] = un[i] + D * dt /dx**2* (un[i+1]-2*un[i]+un[i-1])

#Plot the results
plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

The solution is

Where it is shown how the distribution diffuses across the domain; however, as time passes, the distribution density tends to zero, so that the mass conservative property is not satisfied. If we wish to maintain a constant mass, we have to apply boundary conditions, in this case, we need to set-up a zero flux condition in the extremes of the domain. Then:

$\frac{\partial u}{ \partial t} = 0 $

import numpy as np
import matplotlib.pyplot as plt

#Setting parameters
nx = 101
dx = 0.01
nt = 20000
D = .1
dt = 0.0001 #be aware setting this value

#Initial conditions
u = np.zeros(nx)
u[25:75] = 2
un = u.copy()

#Numerical scheme: Finite differences
for n in range(nt):
    for i in range(1, nx-1):
        un = u.copy()
        u[i] = un[i] + D * dt /dx**2* (un[i+1]-2*un[i]+un[i-1])                                                                                                       # Zero flux condition
        u[0] = u[1]
        u[-1]= u[-2]                                                            
#Plot the results
plt.plot(np.linspace(0, 1, 101), u)
plt.ylim(0, 4)
plt.xlim(0, 1.)
plt.show()

Here, zero flux condition is included straightforward by adding:  u[0] = u[1]   (left)  and
u[-1]= u[-2] (right). Then, the mass can not go outside of the domain.

But!, be careful, because this zero flux condition only works for this simple equation. If you are working with different models, for instance, the advection-diffusion equation, zero flux condition has to be implemented differently.