Method of lines is a numerical technique used to solve a large class of partial differential equations (PDEs). Method to lines uses finite differences methods by discretising all spatial variables of a PDE while time remains continuous. This results in a system of ordinary differential equations (ODEs) which can be solve numerically using solvers such as ode15s() in Matlab.

One prime example of a PDE encountered by many undergraduate students is the diffusion equation. One way to express the diffusion equation is

\[ \frac{\partial \rho(x,t)}{\partial t} = D\frac{\partial^2 \rho(x,t)}{\partial x^2}, \]

with boundary and initial condition as

\[ \rho(0,t) = 0, \quad \rho(100,t) = 0, \]

\[ \rho(x,0) = 0.5[\mathcal{H}(x-40)-\mathcal{H}(x-60)] \]

Here \(\rho\) can be the density, \(t\) is time, \(x\) is spatial domain and is limited between \(0\) and \(100\), \(D\) is the rate of diffusion and \(\mathcal{H}\) is the Heaviside function. Note that we choose \(\rho\) in this context for density of cells as seen in a previous blog.

This well known equation can model many natural phenomenons such as the temperature changing through a medium or tracking cells in biological systems. For example, in the later case researchers have shown the averaged density of cells moving randomly over a lattice can be approximated by the diffusion equation (Binder, Ross, and Simpson 2012).

In general, PDE models have countless of application to the model and better understand everything around us. Hence, it is important to be able to solve such equations and drive meaning from them. In a majority of cases, analytical solutions cannot be found for these complex PDEs models and we have rely on computers to solve it for us. In this case, we will solve it numerically using Matlab as seen below.

Discretisation

Discretisation is the process of a transforming continuous variable(s) into discrete quantities. This technique allows us to solve a substantial class of PDEs numerically. As we want to discretise all variables except time in the method of lines, we can rewrite the diffusion equation as follows

\[ \frac{d\rho_j}{dt} = D\frac{\rho_{j+1} - 2\rho_j + \rho_{j-1}}{\Delta x^2}, j = 1, \dots, M, \]

where \(\Delta x\) is the grid point spacing.

We can use this same approach to solve other more complex PDEs. Wikipedia has a nice page on the finite difference coefficients.

Note the above discretisation can be written in matrix form as

\[ \underbrace{\begin{bmatrix} 1 \quad 0 \quad \cdots \quad 0 \quad 0 \\ 0 \quad 1 \quad \cdots \quad 0 \quad 0 \\ \vdots \quad \vdots \quad \ddots \quad \vdots \quad \vdots \\ 0 \quad 0 \quad \cdots \quad 1 \quad 0 \\ 0 \quad 0 \quad \cdots \quad 0 \quad 1 \end{bmatrix}}_{\text{Left Hand Side (LHS)}} \begin{bmatrix} \frac{d\rho_1}{dt} \\ \frac{d\rho_2}{dt} \\ \vdots \\ \frac{d\rho_{M-1}}{dt} \\ \frac{d\rho_{M}}{dt} \end{bmatrix} = \underbrace{\begin{bmatrix} 0 \\ D\frac{\rho_3-2\rho_2+\rho_1}{\Delta x^2} \\ \vdots \\ D\frac{\rho_{M}-2\rho_{M-1}+\rho_{M-2}}{\Delta x^2} \\ 0 \end{bmatrix}}_{\text{Right Hand Side (RHS)}} \]

Now we have a system of ODEs, we can solve them numerically using MATLAB.

Implementation in MATLAB

We will use MATLAB’s built in ode15s() to solve the system of ODEs. But first we need to set up the discretisation. See code below.

For the first line of code, I general put clear for clearing the workspace to avoid unnecessary bugs in my code and close all figures to minimise confusion when viewing various plots.

clear
close all

Initialize time step to 500. This will be how long in time MATLAB will solve the diffusion equation.

T = 500;

Set up the space variable \(x\) from 1 to 100 with grid spacing of 0.1. Here \(\Delta x = 0.1\).

x = 1:0.1:100;

Assuming all the grid spacing are even, we can evaluate \(\Delta x\) as dx. \(\Delta x\) is simply spacing of one grid point from the neighboring one.

dx = x(2)-x(1);

Also, we will need the number of grid spaces later on, so we assign the number of the grid points as M.

M = length(x);

Furthermore, we will assign j to be the grid points as the boundary conditions implemented differently.

j = 2:M-1;

We now need to code the initial condition. The initial condition is set as the Heaviside to match the simulation from the cellular automata model.

p0 = 0.5*(heaviside(x - 40) - heaviside(x - 60));

Now we set up the discretisation from Discretisation section. Note that we have set \(D = 0.25\) as the diffusion coefficient.

D = 0.25;
dpdt = @(t,p) D*(p(j-1)-2*p(j)+p(j+1))/dx^2;

Now the set up a vector in MATLAB called RHS. This also incorporates the boundary conditions p(1) for \(\rho(0) = 0\) and u(M) for \(\rho(100,t) = 0\).

RHS = @(t,u) [0;...
              dpdt(t,p);...
              0];

As we want to solve the whole system of ODEs we setup as the identity matrix,

LHS = eye(M);

We now need to tell MATLAB to set LHS as the Mass matrix.

opts = odeset('Mass',LHS);

It’s time to solve it! Using ode15s() we input RHS as the discretisation. Choose the time to solve it over, here T = 500. p0 is the initial condition. Finally, we input opts to tell MATALB to we are using the Mass matrix.

[t,p] = ode15s(RHS,linspace(0,T,T),p0,opts);

Plotting it.

plot(x,p(round(t)==10,:),'LineWidth',2)
hold on
plot(x,p(round(t)==96,:),'LineWidth',2)
plot(x,p(round(t)==500,:),'LineWidth',2)
xlabel('$x$', 'interpreter', 'latex', 'FontSize', 20)
ylabel('$\rho(x,t)$', 'interpreter', 'latex', 'FontSize', 20)
legend("$t = 10$", "$t = 100$", "$t = 500$", 'FontSize',15,...
    'interpreter', 'latex')

The solution for the diffusion equation for different time steps is presented above from our simulation in Matlab. At early times we see the shape close following the Heaviside function. However, at later time we see a curved spreading out smoothly. Note the zero boundary conditions are satisfied which is a good sanity check.

References

Binder, Benjamin J, Joshua V Ross, and Matthew J Simpson. 2012. “A Hybrid Model for Studying Spatial Aspects of Infectious Diseases.” The ANZIAM Journal 54 (1-2): 37–49.