Solving the diffusion equation numerically

How are computers used to solve partial differential equations? On this page the diffusion equation is solved numerically using two approaches. First the equation is solved analytically, since the known solution can be derived in this way.. This analytical solution is used to compare several numerical schemes.

Diffusion equation

The diffusion equation is given by, $$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} $$ where $u$ is a variable depending on both time ($t$) and space ($x$), and $D$ is a diffusion coefficient. To solve this equation, initial and boundary conditions are required. The solution of this equation is a function of space and time, and it depends on the diffusion coefficient. An application of this equation could be the heat distribution of a metal bar. In this example, it is assumed that the heat at both ends is zero and that the bar has a length of L ($u(x=0)=0$ and $u(x=L)=0$). One method to solve the equation is the separation of variables approach. The idea is to assume that the solution is given by, $$ u(x,t) = X(x)T(t) $$ with a contribution only dependent on space and a contribution only dependent on time. This results in the following equation, $$ X(x) T(t)^{'} = D X(x)^{''}T(t) $$ where $^{'}$ represents the derivative. Now rearranging all the terms related to space on one side and all terms related to time on the other side, gives, $$ \frac{T(t)^{'} }{DT(t)} = \frac{X(x)^{''}}{X(x)} = -\lambda $$ This can be written as two equations, $$ X(x)^{''} + \lambda X(x) = 0 $$ $$ T(t)^{'} + \lambda T(t) = 0 $$ Both equations are ordinary differential equations for which the solution is known. First, the solution for $X(x)$ is derived. It can be shown that no solution is found for $$\lambda < 0$$, except $u(x,t)=0$. Thus, the solution for $X(x)$ is given by, $$ X(x) = A \sin{\sqrt{\lambda x}} + B\cos{\sqrt{\lambda x}} $$ Using the boundary conditions $X(0)=X(L)=0$, the coefficients can be found ($B=0$ and $\sqrt{\lambda}=n \pi/L$). By substituting $X(x)=e^{rx}$ into the differential equation for $t$ we can find the solution, $$ T(t) = e^{-\lambda D t} $$ Combining both solutions, we find, $$ u(x,t) = \sum_{n=1} C_n \sin{(\frac{n \pi x}{L})}\exp{(-\frac{n^2\pi^2Dt}{L^2})} $$ where $C_n$ is obtained from the initial condition. Given the initial conditions $g(x)$ at $t=0$, we get, $$ g(x) = \sum_{n=1} C_n \sin{(\frac{n \pi x}{L})} $$ This can be seen as a Fourier series of the function $g(x)$. Thus, this can be rewritten as, $$ C_n = 2/L \int_0^L g(x) \sin{(\frac{n \pi x}{L})} dx $$

Numerical solution

n addition to the analytical solution, it is also possible to solve the equation numerically. Two methods are applied: Explicit and Implicit. The explicit method, also called the Euler method, is the easiest to implement. It discretises the equations in the following way, $$ \frac{u^{n+1}-u^{n}}{\Delta t} = D \frac{u^n_{i+1}-2u^n_{i} + u^n_{i-1}}{\Delta x} $$ where $n$ represents the time index and $i$ the spatial index. The method is easy as all terms are taken at the previous time step, making it possible to explicitly solve for the next time step. $$ u^{n+1} = u^{n} + F (u^n_{i+1}-2u^n_{i} + u^n_{i-1}) $$ with $F$ is, $$ F = \frac{D\Delta t}{\Delta x ^2} $$ This method is only stable when $F$ is smaller than 0.5. Alternatively, an implicit scheme can be applied. This scheme takes all the values at the new time step resulting in the follwing scheme, $$ \frac{u^{n+1}-u^{n}}{\Delta t} = D \frac{u^{n+1}_{i+1}-2u^{n+1}_{i} + u^{n+1}_{i-1}}{\Delta x} $$ This scheme is unconditionally stable and does not have limitations in terms of time step or space discretisation, but this comes with additional computational cost. A system of equations has to be solved. This system of equations is given by $$ -F \eta^{n+1}_{i+1} + (1+2F) \eta^{n+1}_{i} - F \eta^{n+1}_{i-1} = \eta^{n}_{i} $$ or in matrix form, $$ \begin{bmatrix} -F & (1+2F) & -F & ... &0 & 0 \\ 0 & -F & (1+2F) & -F & 0 & 0 \\ ... & ... & ... & ... & ...& 0\\ 0 & ... & ... & -F & (1+2)F& -F\\ \end{bmatrix} \begin{bmatrix} \eta^{n+1}_{1} \\ \eta^{n+1}_{2} \\ ... \\ ... \\ ... \\ \eta^{n+1}_{i} \end{bmatrix} = \begin{bmatrix} \eta^{n}_{1} \\ \eta^{n}_{2} \\ ... \\ ... \\ ... \\ \eta^{n}_{i} \end{bmatrix} $$

Numerical solution

let’s solve the diffusion equation with the two methods. notebook