Solving the Diffusion Equation Explicitly

This post is part of a series of Finite Difference Method Articles. Other posts in the series concentrate on Derivative Approximation, the Crank-Nicolson Implicit Method and the Tridiagonal Matrix Solver/Thomas Algorithm:

• Derivative Approximation via Finite Difference Methods
• Solving the Diffusion Equation Explicitly
• Crank-Nicolson Implicit Scheme
• Tridiagonal Matrix Solver via Thomas Algorithm

In Part 1 of the series on Finite Difference Methods, it was shown that continuous derivatives could be approximated and applied to a discrete domain. The next step is to apply these derivatives to a parabolic PDE. The heat equation is the canonical example of a parabolic PDE, and this will now be discretised.

The heat equation with initial condition g is given below:

$$\frac{\partial f}{\partial t}=\frac{{\partial }^{2}f}{\partial x^2},f(x,0)=g(x)$$

This is discretised by applying a forward difference to the time derivative and a centered second difference for the diffusion term to give:

$$\frac{f_i^{n+1}-f_i^n}{\mathrm{\Delta }t}=\frac{f_{i+1}^n-2f_i^n+f_{i-1}^n}{(\mathrm{\Delta }x{)}^2}$$

This equation can be rearranged for $f_i^{n+1}$ to give:

$$f_i^{n+1}=f_i^n+\frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x{)}^2}(f_{i+1}^n-2f_i^n+f_{i-1}^n)$$

Computation of Solution

For $n=1$ all of the approximations to the solution f are known on the right hand side of the equation. This makes the equation explicit. Thus the only remaining task is to determine $\mathrm{\Delta }t$, which is set to $(\mathrm{\Delta }x{)}^2$. Thus the equation reduces to:

$$f_i^{n+1}=f_{i+1}^n-f_i^n+f_{i-1}^n$$

This formula will allow calculation of $f_i^2$ for all $i\in 1,...,I$ as an approximation to g provides the discretised initial condition. Let $g_i$ be such an approximation, with three 'heat bumps' and take $I=13$, then the values of $g_i$ are given by:

$$\begin{array}{ccccccccccccc}0& 0& 1& 0& 0& 1& 2& 1& 0& 0& 1& 0& 0\end{array}$$

This will allow a solution to be generated by progressing from the initial condition in a temporal manner. This is referred to as a time marched solution:

$$\begin{array}{cccccccccccccc}n=1& 0& 0& 1& 0& 0& 1& 2& 1& 0& 0& 1& 0& 0\\ n=2& 0& 1& -1& 1& 1& 1& 0& 1& 1& 1& -1& 1& 0\\ n=3& 1& -2& 3& -1& 1& 0& 2& 0& 1& -1& 3& -2& 1\\ n=4& -3& 6& -6& 5& -2& 3& -2& 3& -2& 5& -6& 6& -3\end{array}$$

For $n=4$ the discretised solution is wildly oscillating and bears no resemblance to the continuous solution, which has the physical interpretation of diffusing the heat bumps. Clearly one of the assumptions made about the formulae is incorrect.

Stability

It turns out that the issue is with the choice of time step size $\mathrm{\Delta }t$ relative to the spatial step size $\mathrm{\Delta }x$. The correct stability criterion is given by:

$$\frac{\mathrm{\Delta }t}{(\mathrm{\Delta }x{)}^2}\le \frac{1}{2}$$

It is beyond the scope of this article to prove this result, but note that it puts a significant restriction on the time step - a problem which will be remedied in the next article.

Applying this new stability criterion, with $\mathrm{\Delta }t=0.5\times (\mathrm{\Delta }x{)}^2$, to the above formulae produces:

$$f_i^{n+1}=\frac{1}{2}{f}_{i+1}^n+\frac{1}{2}{f}_{i-1}^n$$

The calculation for this formulae provides:

$$\begin{array}{cccccccccccccc}n=1& 0& 0& 1& 0& 0& 1& 2& 1& 0& 0& 1& 0& 0\\ n=2& 0& 0.5& 0& 0.5& 0.5& 1& 1& 1& 0.5& 0.5& 0& 0.5& 0\\ n=3& 0.25& 0& 0.5& 0.25& 0.75& 0.75& 1& 0.75& 0.75& 0.25& 0.5& 0& 0.25\\ n=4& 0& 0.375& 0.125& 0.625& 0.5& 0.875& 0.75& 0.875& 0.5& 0.625& 0.125& 0.375& 0\end{array}$$

This is an accurate approximation to the true solution. There are no negative values and the physical interpretation of the heat diffusing through a 1D bar fits with the solution.

Conclusion

The choice of time step is very restrictive. In practice, with a refined spatial mesh, the time step would be proihibitively small. The next article will introduce another method to overcome the small timestep issue.

Reference: Solving the Diffusion Equation Explicitly

From https://www.quantstart.com/articles/Solving-the-Diffusion-Equation-Explicitly/