The famous diffusion equation, also known as the heat equation, reads

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

where $u(x,t)$ is the unknown function to be solved for, $x$ is a coordinate in space, and $t$ is time. The coefficient $\dfc$ is the diffusion coefficient and determines how fast $u$ changes in time. A quick short form for the diffusion equation is $u_t = \dfc u_{xx}$.

Compared to the wave equation, $u_{tt}=c^2u_{xx}$, which looks very similar, the diffusion equation features solutions that are very different from those of the wave equation. Also, the diffusion equation makes quite different demands to the numerical methods.

Typical diffusion problems may experience rapid change in the very beginning, but then the evolution of $u$ becomes slower and slower. The solution is usually very smooth, and after some time, one cannot recognize the initial shape of $u$. This is in sharp contrast to solutions of the wave equation where the initial shape is preserved in homogeneous media - the solution is then basically a moving initial condition. The standard wave equation $u_{tt}=c^2u_{xx}$ has solutions that propagate with speed $c$ forever, without changing shape, while the diffusion equation converges to a stationary solution $\bar u(x)$ as $t\rightarrow\infty$. In this limit, $u_t=0$, and $\bar u$ is governed by $\bar u''(x)=0$. This stationary limit of the diffusion equation is called the Laplace equation and arises in a very wide range of applications throughout the sciences.

mathcal{I}_t is possible to solve for $u(x,t)$ using an explicit scheme, as we do in the section An explicit method for the 1D diffusion equation, but the time step restrictions soon become much less favorable than for an explicit scheme applied to the wave equation. And of more importance, since the solution $u$ of the diffusion equation is very smooth and changes slowly, small time steps are not convenient and not required by accuracy as the diffusion process converges to a stationary state. Therefore, implicit schemes (as described in the section Implicit methods for the 1D diffusion equation) are popular, but these require solutions of systems of algebraic equations. We shall use ready-made software for this purpose, but also program some simple iterative methods. The exposition is, as usual in this book, very basic and focuses on the basic ideas and how to implement. More comprehensive mathematical treatments and classical analysis of the methods are found in lots of textbooks. A favorite of ours in this respect is the one by LeVeque [LeVeque_2007]. The books by Strikwerda [Strikwerda_2007] and by Lapidus and Pinder [Lapidus_Pinder_1982] are also highly recommended as additional material on the topic.

An explicit method for the 1D diffusion equation¶

Explicit finite difference methods for the wave equation $u_{tt}=c^2u_{xx}$ can be used, with small modifications, for solving $u_t = \dfc u_{xx}$ as well. % if BOOK == "book": The exposition below assumes that the reader is familiar with the basic ideas of discretization and implementation of wave equations from the chapter ch:wave. Readers not familiar with the Forward Euler, Backward Euler, and Crank-Nicolson (or centered or midpoint) discretization methods in time should consult, e.g., Section 1.1 in [Langtangen_decay]. % endif

The initial-boundary value problem for 1D diffusion¶

To obtain a unique solution of the diffusion equation, or equivalently, to apply numerical methods, we need initial and boundary conditions. The diffusion equation goes with one initial condition $u(x,0)=I(x)$, where $I$ is a prescribed function. One boundary condition is required at each point on the boundary, which in 1D means that $u$ must be known, $u_x$ must be known, or some combination of them.

We shall start with the simplest boundary condition: $u=0$. The complete initial-boundary value diffusion problem in one space dimension can then be specified as

$$ \begin{equation} \frac{\partial u}{\partial t} = \dfc \frac{\partial^2 u}{\partial x^2} + f, \quad x\in (0,L),\ t\in (0,T] \label{diffu:pde1} \tag{1} \end{equation} $$

$$ \begin{equation} u(x,0) = I(x), \quad x\in [0,L] \label{diffu:pde1:ic:u} \tag{2} \end{equation} $$

$$ \begin{equation} u(0,t) = 0, \quad t>0, \label{diffu:pde1:bc:0} \tag{3} \end{equation} $$

$$ \begin{equation} u(L,t) = 0, \quad t>0\thinspace . \label{diffu:pde1:bc:L} \tag{4} \end{equation} $$

With only a first-order derivative in time, only one initial condition is needed, while the second-order derivative in space leads to a demand for two boundary conditions. We have added a source term $f=f(x,t)$, which is convenient when testing implementations.

Diffusion equations like (1) have a wide range of applications throughout physical, biological, and financial sciences. One of the most common applications is propagation of heat, where $u(x,t)$ represents the temperature of some substance at point $x$ and time $t$. Other applications are listed in the section diffu:app.

Forward Euler scheme¶

The first step in the discretization procedure is to replace the domain $[0,L]\times [0,T]$ by a set of mesh points. Here we apply equally spaced mesh points

$$ x_i=i\Delta x,\quad i=0,\ldots,N_x, $$

and

$$ t_n=n\Delta t,\quad n=0,\ldots,N_t \thinspace . $$

Moreover, $u^n_i$ denotes the mesh function that approximates $u(x_i,t_n)$ for $i=0,\ldots,N_x$ and $n=0,\ldots,N_t$. Requiring the PDE (1) to be fulfilled at a mesh point $(x_i,t_n)$ leads to the equation

$$ \begin{equation} \frac{\partial}{\partial t} u(x_i, t_n) = \dfc\frac{\partial^2}{\partial x^2} u(x_i, t_n) + f(x_i,t_n), \label{diffu:pde1:step2} \tag{5} \end{equation} $$

The next step is to replace the derivatives by finite difference approximations. The computationally simplest method arises from using a forward difference in time and a central difference in space:

$$ \begin{equation} [D_t^+ u = \dfc D_xD_x u + f]^n_i \thinspace . \label{diffu:pde1:step3a} \tag{6} \end{equation} $$

Written out,

$$ \begin{equation} \frac{u^{n+1}_i-u^n_i}{\Delta t} = \dfc \frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\Delta x^2} + f_i^n\thinspace . \label{diffu:pde1:step3b} \tag{7} \end{equation} $$

We have turned the PDE into algebraic equations, also often called discrete equations. The key property of the equations is that they are algebraic, which makes them easy to solve. As usual, we anticipate that $u^n_i$ is already computed such that $u^{n+1}_i$ is the only unknown in (7). Solving with respect to this unknown is easy:

$$ \begin{equation} u^{n+1}_i = u^n_i + F\left( u^{n}_{i+1} - 2u^n_i + u^n_{i-1}\right) + \Delta t f_i^n, \label{diffu:pde1:step4} \tag{8} \end{equation} $$

where we have introduced the mesh Fourier number:

$$ \begin{equation} F = \dfc\frac{\Delta t}{\Delta x^2}\thinspace . \label{_auto1} \tag{9} \end{equation} $$

$F$ is the key parameter in the discrete diffusion equation.

Note that $F$ is a dimensionless number that lumps the key physical parameter in the problem, $\dfc$, and the discretization parameters $\Delta x$ and $\Delta t$ into a single parameter. Properties of the numerical method are critically dependent upon the value of $F$ (see the section diffu:pde1:analysis for details).

The computational algorithm then becomes

1. compute $u^0_i=I(x_i)$ for $i=0,\ldots,N_x$

2. for $n=0,1,\ldots,N_t$:

a. apply (8) for all the internal spatial points $i=1,\ldots,N_x-1$

b. set the boundary values $u^{n+1}_i=0$ for $i=0$ and $i=N_x$

The algorithm is compactly and fully specified in Python:

Note that we use a for $\dfc$ in the code, motivated by easy visual mapping between the variable name and the mathematical symbol in formulas.

We need to state already now that the shown algorithm does not produce meaningful results unless $F\leq 1/2$. Why is explained in the section diffu:pde1:analysis.

Implementation¶

The file diffu1D_u0.py contains a complete function solver_FE_simple for solving the 1D diffusion equation with $u=0$ on the boundary as specified in the algorithm above:

A faster alternative is available in the function solver_FE, which adds the possibility of solving the finite difference scheme by vectorization. The vectorized version replaces the explicit loop

by arithmetics on displaced slices of the u array:

For example, the vectorized version runs 70 times faster than the scalar version in a case with 100 time steps and a spatial mesh of $10^5$ cells.

The solver_FE function also features a callback function such that the user can process the solution at each time level. The callback function looks like user_action(u, x, t, n), where u is the array containing the solution at time level n, x holds all the spatial mesh points, while t holds all the temporal mesh points. The solver_FE function is very similar to solver_FE_simple above:

Verification¶

Exact solution of discrete equations¶

Before thinking about running the functions in the previous section, we need to construct a suitable test example for verification. mathcal{I}_t appears that a manufactured solution that is linear in time and at most quadratic in space fulfills the Forward Euler scheme exactly. With the restriction that $u=0$ for $x=0,L$, we can try the solution

$$ u(x,t) = 5tx(L-x)\thinspace . $$

Inserted in the PDE, it requires a source term

$$ f(x,t) = 10\dfc t + 5x(L-x)\thinspace . $$

% if BOOK == 'book': With the formulas from sec:form:fdtn we can easily check % else: Let us check % endif that the manufactured u fulfills the scheme:

$$ \begin{align*} \lbrack D_t^+ u = \dfc D_x D_x u + f\rbrack^n_i &= \lbrack 5x(L-x)D_t^+ t = 5 t\dfc D_x D_x (xL-x^2) +\\ &\quad\quad 10\dfc t + 5x(L-x)\rbrack^n_i\\ &= \lbrack 5x(L-x) = 5 t\dfc (-2) + 10\dfc t + 5x(L-x) \rbrack^n_i, \end{align*} $$

which is a 0=0 expression. The computation of the source term, given any $u$, is easily automated with sympy:

Now we can choose any expression for u and automatically get the suitable source term f. However, the manufactured solution u will in general not be exactly reproduced by the scheme: only constant and linear functions are differentiated correctly by a forward difference, while only constant, linear, and quadratic functions are differentiated exactly by a $[D_xD_x u]^n_i$ difference.

The numerical code will need to access the u and f above as Python functions. The exact solution is wanted as a Python function u_exact(x, t), while the source term is wanted as f(x, t). The parameters a and L in u and f above are symbols and must be replaced by float objects in a Python function. This can be done by redefining a and L as float objects and performing substitutions of symbols by numbers in u and f. The appropriate code looks like this:

Here we also make a function I for the initial condition.

The idea now is that our manufactured solution should be exactly reproduced by the code (to machine precision). For this purpose we make a test function for comparing the exact and numerical solutions at the end of the time interval:

The critical value $F=0.5$.

We emphasize that the value F=0.5 is critical: the tests above will fail if F has a larger value. This is because the Forward Euler scheme is unstable for $F>1/2$.

The reader may wonder if $F=1/2$ is safe or if $F<1/2$ should be required. Experiments show that $F=1/2$ works fine for $u_t=\dfc u_{xx}$, so there is no accumulation of rounding errors in this case and hence no need to introduce any safety factor to keep $F$ away from the limiting value 0.5.

Checking convergence rates¶

If our chosen exact solution does not satisfy the discrete equations exactly, we are left with checking the convergence rates, just as we did previously for the wave equation. However, with the Euler scheme here, we have different accuracies in time and space, since we use a second order approximation to the spatial derivative and a first order approximation to the time derivative. Thus, we must expect different convergence rates in time and space. For the numerical error,

$$ E = C_t\Delta t^r + C_x\Delta x^p, $$

we should get convergence rates $r=1$ and $p=2$ ($C_t$ and $C_x$ are unknown constants). As previously, in the section wave:pde2:fd:MMS, we simplify matters by introducing a single discretization parameter $h$:

$$ h = \Delta t,\quad \Delta x = Kh^{r/p}, $$

where $K$ is any constant. This allows us to factor out only one discretization parameter $h$ from the formula:

$$ E = C_t h + C_x (Kh^{r/p})^p = \tilde C h^r,\quad \tilde C = C_t + C_sK^r\thinspace . $$

The computed rate $r$ should approach 1 with increasing resolution.

mathcal{I}_t is tempting, for simplicity, to choose $K=1$, which gives $\Delta x = h^{r/p}$, expected to be $\sqrt{\Delta t}$. However, we have to control the stability requirement: $F\leq\frac{1}{2}$, which means

$$ \frac{\dfc\Delta t}{\Delta x^2}\leq\frac{1}{2}\quad\Rightarrow \quad \Delta x \geq \sqrt{2\dfc}h^{1/2} , $$

implying that $K=\sqrt{2\dfc}$ is our choice in experiments where we lie on the stability limit $F=1/2$.

Numerical experiments¶

When a test function like the one above runs silently without errors, we have some evidence for a correct implementation of the numerical method. The next step is to do some experiments with more interesting solutions.

We target a scaled diffusion problem where $x/L$ is a new spatial coordinate and $\dfc t/L^2$ is a new time coordinate. The source term $f$ is omitted, and $u$ is scaled by $\max_{x\in [0,L]}|I(x)|$ (see Section 3.2 in [Langtangen_scaling] for details). The governing PDE is then

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

in the spatial domain $[0,L]$, with boundary conditions $u(0)=u(1)=0$. Two initial conditions will be tested: a discontinuous plug,

$$ I(x) = \left\lbrace\begin{array}{ll} 0, & |x-L/2| > 0.1\\ 1, & \hbox{otherwise} \end{array}\right. $$

and a smooth Gaussian function,

$$ I(x) = e^{-\frac{1}{2\sigma^2}(x-L/2)^2}\thinspace . $$

The functions plug and gaussian in diffu1D_u0.py run the two cases, respectively:

These functions make use of the function viz for running the solver and visualizing the solution using a callback function with plotting:

Notice that this viz function stores all the solutions in a list solutions in the callback function. Modern computers have hardly any problem with storing a lot of such solutions for moderate values of $N_x$ in 1D problems, but for 2D and 3D problems, this technique cannot be used and solutions must be stored in files.

[hpl 1: Better to show the scalable file solution here?]

Our experiments employ a time step $\Delta t = 0.0002$ and simulate for $t\in [0,0.1]$. First we try the highest value of $F$: $F=0.5$. This resolution corresponds to $N_x=50$. A possible terminal command is

    Terminal> python -c 'from diffu1D_u0 import gaussian
              gaussian("solver_FE", F=0.5, dt=0.0002)'

The $u(x,t)$ curve as a function of $x$ is shown in Figure at four time levels.

We see that the curves have saw-tooth waves in the beginning of the simulation. This non-physical noise is smoothed out with time, but solutions of the diffusion equations are known to be smooth, and this numerical solution is definitely not smooth. Lowering $F$ helps: $F\leq 0.25$ gives a smooth solution, see % if FORMAT == "pdflatex": Figure (and a movie). % else: Figure.

% endif

Increasing $F$ slightly beyond the limit 0.5, to $F=0.51$, gives growing, non-physical instabilities, as seen in Figure.

Forward Euler scheme for $F=0.5$.

Forward Euler scheme for $F=0.25$.

Forward Euler scheme for $F=0.51$.

Instead of a discontinuous initial condition we now try the smooth Gaussian function for $I(x)$. A simulation for $F=0.5$ is shown in Figure. Now the numerical solution is smooth for all times, and this is true for any $F\leq 0.5$.

% if FORMAT != "pdflatex":

% endif

Forward Euler scheme for $F=0.5$.

Experiments with these two choices of $I(x)$ reveal some important observations:

• The Forward Euler scheme leads to growing solutions if $F>\frac{1}{2}$.

• $I(x)$ as a discontinuous plug leads to a saw tooth-like noise for $F=\frac{1}{2}$, which is absent for $F\leq\frac{1}{4}$.

• The smooth Gaussian initial function leads to a smooth solution for all relevant $F$ values ($F\leq \frac{1}{2}$).

Implicit methods for the 1D diffusion equation¶

Simulations with the Forward Euler scheme show that the time step restriction, $F\leq\frac{1}{2}$, which means $\Delta t \leq \Delta x^2/(2\dfc)$, may be relevant in the beginning of the diffusion process, when the solution changes quite fast, but as time increases, the process slows down, and a small $\Delta t$ may be inconvenient. With implicit schemes, which lead to coupled systems of linear equations to be solved at each time level, any size of $\Delta t$ is possible (but the accuracy decreases with increasing $\Delta t$). The Backward Euler scheme, derived and implemented below, is the simplest implicit scheme for the diffusion equation.

Backward Euler scheme¶

In (5), we now apply a backward difference in time, but the same central difference in space:

$$ \begin{equation} [D_t^- u = D_xD_x u + f]^n_i, \label{diffu:pde1:step3aBE} \tag{10} \end{equation} $$

which written out reads

$$ \begin{equation} \frac{u^{n}_i-u^{n-1}_i}{\Delta t} = \dfc\frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\Delta x^2} + f_i^n\thinspace . \label{diffu:pde1:step3bBE} \tag{11} \end{equation} $$

Now we assume $u^{n-1}_i$ is already computed, but that all quantities at the "new" time level $n$ are unknown. This time it is not possible to solve with respect to $u_i^{n}$ because this value couples to its neighbors in space, $u^n_{i-1}$ and $u^n_{i+1}$, which are also unknown. Let us examine this fact for the case when $N_x=3$. Equation (11) written for $i=1,\ldots,Nx-1= 1,2$ becomes

$$ \begin{equation} \frac{u^{n}_1-u^{n-1}_1}{\Delta t} = \dfc\frac{u^{n}_{2} - 2u^n_1 + u^n_{0}}{\Delta x^2} + f_1^n \label{_auto2} \tag{12} \end{equation} $$

$$ \begin{equation} \frac{u^{n}_2-u^{n-1}_2}{\Delta t} = \dfc\frac{u^{n}_{3} - 2u^n_2 + u^n_{1}}{\Delta x^2} + f_2^n \label{_auto3} \tag{13} \end{equation} $$

The boundary values $u^n_0$ and $u^n_3$ are known as zero. Collecting the unknown new values $u^n_1$ and $u^n_2$ on the left-hand side and multiplying by $\Delta t$ gives

$$ \begin{equation} \left(1+ 2F\right) u^{n}_1 - F u^{n}_{2} = u^{n-1}_1 + \Delta t f_1^n, \label{_auto4} \tag{14} \end{equation} $$

$$ \begin{equation} - F u^{n}_{1} + \left(1+ 2F\right) u^{n}_2 = u^{n-1}_2 + \Delta t f_2^n\thinspace . \label{_auto5} \tag{15} \end{equation} $$

This is a coupled $2\times 2$ system of algebraic equations for the unknowns $u^n_1$ and $u^n_2$. The equivalent matrix form is

$$ \left(\begin{array}{cc} 1+ 2F & - F\\ - F & 1+ 2F \end{array}\right) \left(\begin{array}{c} u^{n}_1\\ u^{n}_2 \end{array}\right) = \left(\begin{array}{c} u^{n-1}_1 + \Delta t f_1^n\\ u^{n-1}_2 + \Delta t f_2^n \end{array}\right) $$

Terminology: implicit vs. explicit methods.

Discretization methods that lead to a coupled system of equations for the unknown function at a new time level are said to be implicit methods. The counterpart, explicit methods, refers to discretization methods where there is a simple explicit formula for the values of the unknown function at each of the spatial mesh points at the new time level. From an implementational point of view, implicit methods are more comprehensive to code since they require the solution of coupled equations, i.e., a matrix system, at each time level. With explicit methods we have a closed-form formula for the value of the unknown at each mesh point.

Very often explicit schemes have a restriction on the size of the time step that can be relaxed by using implicit schemes. In fact, implicit schemes are frequently unconditionally stable, so the size of the time step is governed by accuracy and not by stability. This is the great advantage of implicit schemes.

In the general case, (11) gives rise to a coupled $(N_x-1)\times (N_x-1)$ system of algebraic equations for all the unknown $u^n_i$ at the interior spatial points $i=1,\ldots,N_x-1$. Collecting the unknowns on the left-hand side, (11) can be written

$$ \begin{equation} - F u^n_{i-1} + \left(1+ 2F \right) u^{n}_i - F u^n_{i+1} = u_{i-1}^{n-1}, \label{diffu:pde1:step4BE} \tag{16} \end{equation} $$

for $i=1,\ldots,N_x-1$. One can either view these equations as a system where the $u^{n}_i$ values at the internal mesh points, $i=1,\ldots,N_x-1$, are unknown, or we may append the boundary values $u^n_0$ and $u^n_{N_x}$ to the system. In the latter case, all $u^n_i$ for $i=0,\ldots,N_x$ are considered unknown, and we must add the boundary equations to the $N_x-1$ equations in (16):

$$ \begin{equation} u_0^n = 0,\label{diffu:pde1:step4BE:BC:0} \tag{17} \end{equation} $$

$$ \begin{equation} u_{N_x}^n = 0\thinspace . \label{diffu:pde1:step4BE:BC:L} \tag{18} \end{equation} $$

A coupled system of algebraic equations can be written on matrix form, and this is important if we want to call up ready-made software for solving the system. The equations (16) and (17)--(18) correspond to the matrix equation

$$ AU = b $$

where $U=(u^n_0,\ldots,u^n_{N_x})$, and the matrix $A$ has the following structure:

$$ \begin{equation} A = \left( \begin{array}{cccccccccc} A_{0,0} & A_{0,1} & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ A_{1,0} & A_{1,1} & A_{1,2} & \ddots & & & & & \vdots \\ 0 & A_{2,1} & A_{2,2} & A_{2,3} & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & A_{i,i-1} & A_{i,i} & A_{i,i+1} & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & A_{N_x-1,N_x} \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & A_{N_x,N_x-1} & A_{N_x,N_x} \end{array} \right) \label{diffu:pde1:matrix:sparsity} \tag{19} \end{equation} $$

The nonzero elements are given by

$$ \begin{equation} A_{i,i-1} = -F \label{_auto6} \tag{20} \end{equation} $$

$$ \begin{equation} A_{i,i} = 1+ 2F \label{_auto7} \tag{21} \end{equation} $$

$$ \begin{equation} A_{i,i+1} = -F \label{_auto8} \tag{22} \end{equation} $$

in the equations for internal points, $i=1,\ldots,N_x-1$. The first and last equation correspond to the boundary condition, where we know the solution, and therefore we must have

$$ \begin{equation} A_{0,0} = 1, \label{_auto9} \tag{23} \end{equation} $$

$$ \begin{equation} A_{0,1} = 0, \label{_auto10} \tag{24} \end{equation} $$

$$ \begin{equation} A_{N_x,N_x-1} = 0, \label{_auto11} \tag{25} \end{equation} $$

$$ \begin{equation} A_{N_x,N_x} = 1\thinspace . \label{_auto12} \tag{26} \end{equation} $$

The right-hand side $b$ is written as

$$ \begin{equation} b = \left(\begin{array}{c} b_0\\ b_1\\ \vdots\\ b_i\\ \vdots\\ b_{N_x} \end{array}\right) \label{_auto13} \tag{27} \end{equation} $$

with

$$ \begin{equation} b_0 = 0, \label{_auto14} \tag{28} \end{equation} $$

$$ \begin{equation} b_i = u^{n-1}_i,\quad i=1,\ldots,N_x-1, \label{_auto15} \tag{29} \end{equation} $$

$$ \begin{equation} b_{N_x} = 0 \thinspace . \label{_auto16} \tag{30} \end{equation} $$

We observe that the matrix $A$ contains quantities that do not change in time. Therefore, $A$ can be formed once and for all before we enter the recursive formulas for the time evolution. The right-hand side $b$, however, must be updated at each time step. This leads to the following computational algorithm, here sketched with Python code:

Regarding verification, the same considerations apply as for the Forward Euler method (the section Verification).

Sparse matrix implementation¶

We have seen from (19) that the matrix $A$ is tridiagonal. The code segment above used a full, dense matrix representation of $A$, which stores a lot of values we know are zero beforehand, and worse, the solution algorithm computes with all these zeros. With $N_x+1$ unknowns, the work by the solution algorithm is $\frac{1}{3} (N_x+1)^3$ and the storage requirements $(N_x+1)^2$. By utilizing the fact that $A$ is tridiagonal and employing corresponding software tools that work with the three diagonals, the work and storage demands can be proportional to $N_x$ only. This leads to a dramatic improvement: with $N_x=200$, which is a realistic resolution, the code runs about 40,000 times faster and reduces the storage to just 1.5%! mathcal{I}_t is no doubt that we should take advantage of the fact that $A$ is tridiagonal.

The key idea is to apply a data structure for a tridiagonal or sparse matrix. The scipy.sparse package has relevant utilities. For example, we can store only the nonzero diagonals of a matrix. The package also has linear system solvers that operate on sparse matrix data structures. The code below illustrates how we can store only the main diagonal and the upper and lower diagonals.

The scipy.sparse.linalg.spsolve function utilizes the sparse storage structure of A and performs, in this case, a very efficient Gaussian elimination solve.

The program diffu1D_u0.py contains a function solver_BE, which implements the Backward Euler scheme sketched above. As mentioned in the section Forward Euler scheme, the functions plug and gaussian run the case with $I(x)$ as a discontinuous plug or a smooth Gaussian function. All experiments point to two characteristic features of the Backward Euler scheme: 1) it is always stable, and 2) it always gives a smooth, decaying solution.

Crank-Nicolson scheme¶

The idea in the Crank-Nicolson scheme is to apply centered differences in space and time, combined with an average in time. We demand the PDE to be fulfilled at the spatial mesh points, but midway between the points in the time mesh:

$$ \frac{\partial}{\partial t} u(x_i, t_{n+\frac{1}{2}}) = \dfc\frac{\partial^2}{\partial x^2}u(x_i, t_{n+\frac{1}{2}}) + f(x_i,t_{n+\frac{1}{2}}), $$

for $i=1,\ldots,N_x-1$ and $n=0,\ldots, N_t-1$.

With centered differences in space and time, we get

$$ [D_t u = \dfc D_xD_x u + f]^{n+\frac{1}{2}}_i\thinspace . $$

On the right-hand side we get an expression

$$ \frac{1}{\Delta x^2}\left(u^{n+\frac{1}{2}}_{i-1} - 2u^{n+\frac{1}{2}}_i + u^{n+\frac{1}{2}}_{i+1}\right) + f_i^{n+\frac{1}{2}}\thinspace . $$

This expression is problematic since $u^{n+\frac{1}{2}}_i$ is not one of the unknowns we compute. A possibility is to replace $u^{n+\frac{1}{2}}_i$ by an arithmetic average:

$$ u^{n+\frac{1}{2}}_i\approx \frac{1}{2}\left(u^{n}_i +u^{n+1}_{i}\right)\thinspace . $$

In the compact notation, we can use the arithmetic average notation $\overline{u}^t$:

$$ [D_t u = \dfc D_xD_x \overline{u}^t + f]^{n+\frac{1}{2}}_i\thinspace . $$

We can also use an average for $f_i^{n+\frac{1}{2}}$:

$$ [D_t u = \dfc D_xD_x \overline{u}^t + \overline{f}^t]^{n+\frac{1}{2}}_i\thinspace . $$

After writing out the differences and average, multiplying by $\Delta t$, and collecting all unknown terms on the left-hand side, we get

$$ u^{n+1}_i - \frac{1}{2} F(u^{n+1}_{i-1} - 2u^{n+1}_i + u^{n+1}_{i+1}) = u^{n}_i + \frac{1}{2} F(u^{n}_{i-1} - 2u^{n}_i + u^{n}_{i+1})\nonumber $$

$$ \begin{equation} \qquad + \frac{1}{2} f_i^{n+1} + \frac{1}{2} f_i^n\thinspace . \label{_auto17} \tag{31} \end{equation} $$

Also here, as in the Backward Euler scheme, the new unknowns $u^{n+1}_{i-1}$, $u^{n+1}_{i}$, and $u^{n+1}_{i+1}$ are coupled in a linear system $AU=b$, where $A$ has the same structure as in (19), but with slightly different entries:

$$ \begin{equation} A_{i,i-1} = -\frac{1}{2} F \label{_auto18} \tag{32} \end{equation} $$

$$ \begin{equation} A_{i,i} = 1 + F \label{_auto19} \tag{33} \end{equation} $$

$$ \begin{equation} A_{i,i+1} = -\frac{1}{2} F \label{_auto20} \tag{34} \end{equation} $$

in the equations for internal points, $i=1,\ldots,N_x-1$. The equations for the boundary points correspond to

$$ \begin{equation} A_{0,0} = 1, \label{_auto21} \tag{35} \end{equation} $$

$$ \begin{equation} A_{0,1} = 0, \label{_auto22} \tag{36} \end{equation} $$

$$ \begin{equation} A_{N_x,N_x-1} = 0, \label{_auto23} \tag{37} \end{equation} $$

$$ \begin{equation} A_{N_x,N_x} = 1\thinspace . \label{_auto24} \tag{38} \end{equation} $$

The right-hand side $b$ has entries

$$ \begin{equation} b_0 = 0, \label{_auto25} \tag{39} \end{equation} $$

$$ \begin{equation} b_i = u^{n-1}_i + \frac{1}{2}(f_i^n + f_i^{n+1}),\quad i=1,\ldots,N_x-1, \label{_auto26} \tag{40} \end{equation} $$

$$ \begin{equation} b_{N_x} = 0 \thinspace . \label{_auto27} \tag{41} \end{equation} $$

When verifying some implementation of the Crank-Nicolson scheme by convergence rate testing, one should note that the scheme is second order accurate in both space and time. The numerical error then reads

$$ E = C_t\Delta t^r + C_x\Delta x^r, $$

where $r=2$ ($C_t$ and $C_x$ are unknown constants, as before). When introducing a single discretization parameter, we may now simply choose

$$ h = \Delta x = \Delta t, $$

which gives

$$ E = C_th^r + C_xh^r = (C_t + C_x)h^r, $$

where $r$ should approach 2 as resolution is increased in the convergence rate computations.

The unifying $\theta$ rule¶

For the equation

$$ \frac{\partial u}{\partial t} = G(u), $$

where $G(u)$ is some spatial differential operator, the $\theta$-rule looks like

$$ \frac{u^{n+1}_i - u^n_i}{\Delta t} = \theta G(u^{n+1}_i) + (1-\theta) G(u^{n}_i)\thinspace . $$

The important feature of this time discretization scheme is that we can implement one formula and then generate a family of well-known and widely used schemes:

• $\theta=0$ gives the Forward Euler scheme in time

• $\theta=1$ gives the Backward Euler scheme in time

• $\theta=\frac{1}{2}$ gives the Crank-Nicolson scheme in time

In the compact difference notation, we write the $\theta$ rule as

$$ [\overline{D}_t u = \dfc D_xD_x u]^{n+\theta}\thinspace . $$

We have that $t_{n+\theta} = \theta t_{n+1} + (1-\theta)t_n$.

Applied to the 1D diffusion problem, the $\theta$-rule gives

$$ \begin{align*} \frac{u^{n+1}_i-u^n_i}{\Delta t} &= \dfc\left( \theta \frac{u^{n+1}_{i+1} - 2u^{n+1}_i + u^{n+1}_{i-1}}{\Delta x^2} + (1-\theta) \frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\Delta x^2}\right)\\ &\qquad + \theta f_i^{n+1} + (1-\theta)f_i^n \thinspace . \end{align*} $$

This scheme also leads to a matrix system with entries

$$ A_{i,i-1} = -F\theta,\quad A_{i,i} = 1+2F\theta\quad, A_{i,i+1} = -F\theta, $$

while right-hand side entry $b_i$ is

$$ b_i = u^n_{i} + F(1-\theta) \frac{u^{n}_{i+1} - 2u^n_i + u^n_{i-1}}{\Delta x^2} + \Delta t\theta f_i^{n+1} + \Delta t(1-\theta)f_i^n\thinspace . $$

The corresponding entries for the boundary points are as in the Backward Euler and Crank-Nicolson schemes listed earlier.

Note that convergence rate testing with implementations of the theta rule must adjust the error expression according to which of the underlying schemes is actually being run. That is, if $\theta=0$ (i.e., Forward Euler) or $\theta=1$ (i.e., Backward Euler), there should be first order convergence, whereas with $\theta=0.5$ (i.e., Crank-Nicolson), one should get second order convergence (as outlined in previous sections).

Experiments¶

We can repeat the experiments from the section Numerical experiments to see if the Backward Euler or Crank-Nicolson schemes have problems with sawtooth-like noise when starting with a discontinuous initial condition. We can also verify that we can have $F>\frac{1}{2}$, which allows larger time steps than in the Forward Euler method.

Backward Euler scheme for $F=0.5$.

The Backward Euler scheme always produces smooth solutions for any $F$. Figure shows one example. Note that the mathematical discontinuity at $t=0$ leads to a linear variation on a mesh, but the approximation to a jump becomes better as $N_x$ increases. In our simulation, we specify $\Delta t$ and $F$, and set $N_x$ to $L/\sqrt{\dfc\Delta t/F}$. Since $N_x\sim\sqrt{F}$, the discontinuity looks sharper in the Crank-Nicolson simulations with larger $F$.

The Crank-Nicolson method produces smooth solutions for small $F$, $F\leq\frac{1}{2}$, but small noise gets more and more evident as $F$ increases. Figures diffu:pde1:CN:fig:F=3 and diffu:pde1:CN:fig:F=10 demonstrate the effect for $F=3$ and $F=10$, respectively. The section diffu:pde1:analysis explains why such noise occur.

Crank-Nicolson scheme for $F=3$.

Crank-Nicolson scheme for $F=10$.

The Laplace and Poisson equation¶

The Laplace equation, $\nabla^2 u = 0$, and the Poisson equation, $-\nabla^2 u = f$, occur in numerous applications throughout science and engineering. In 1D these equations read $u''(x)=0$ and $-u''(x)=f(x)$, respectively. We can solve 1D variants of the Laplace equations with the listed software, because we can interpret $u_{xx}=0$ as the limiting solution of $u_t = \dfc u_{xx}$ when $u$ reaches a steady state limit where $u_t\rightarrow 0$. Similarly, Poisson's equation $-u_{xx}=f$ arises from solving $u_t = u_{xx} + f$ and letting $t\rightarrow\infty$ so $u_t\rightarrow 0$.

Technically in a program, we can simulate $t\rightarrow\infty$ by just taking one large time step: $\Delta t\rightarrow\infty$. In the limit, the Backward Euler scheme gives

$$ -\frac{u^{n+1}_{i+1} - 2u^{n+1}_i + u^{n+1}_{i-1}}{\Delta x^2} = f^{n+1}_i, $$

which is nothing but the discretization $[-D_xD_x u = f]^{n+1}_i=0$ of $-u_{xx}=f$.

The result above means that the Backward Euler scheme can solve the limit equation directly and hence produce a solution of the 1D Laplace equation. With the Forward Euler scheme we must do the time stepping since $\Delta t > \Delta x^2/\dfc$ is illegal and leads to instability. We may interpret this time stepping as solving the equation system from $-u_{xx}=f$ by iterating on a pseudo time variable.

[hpl 2: Better to say the last sentence when we treat iterative methods.]