Analysis of schemes for the diffusion equation¶

The numerical experiments in the sections diffu:pde1:FE:experiments and diffu:pde1:theta:experiments reveal that there are some numerical problems with the Forward Euler and Crank-Nicolson schemes: sawtooth-like noise is sometimes present in solutions that are, from a mathematical point of view, expected to be smooth. This section presents a mathematical analysis that explains the observed behavior and arrives at criteria for obtaining numerical solutions that reproduce the qualitative properties of the exact solutions. In short, we shall explain what is observed in Figures diffu:pde1:FE:fig:F=0.5-diffu:pde1:CN:fig:F=10.

Properties of the solution¶

A particular characteristic of diffusive processes, governed by an equation like

$$ \begin{equation} u_t = \dfc u_{xx}, \label{diffu:pde1:eq} \tag{1} \end{equation} $$

is that the initial shape $u(x,0)=I(x)$ spreads out in space with time, along with a decaying amplitude. Three different examples will illustrate the spreading of $u$ in space and the decay in time.

Similarity solution¶

The diffusion equation (1) admits solutions that depend on $\eta = (x-c)/\sqrt{4\dfc t}$ for a given value of $c$. One particular solution is

$$ \begin{equation} u(x,t) = a\,\mbox{erf}(\eta) + b, \label{diffu:pdf1:erf:sol} \tag{2} \end{equation} $$

where

$$ \begin{equation} \mbox{erf}(\eta) = \frac{2}{\sqrt{\pi}}\int_0^\eta e^{-\zeta^2}d\zeta, \label{diffu:analysis:erf:def} \tag{3} \end{equation} $$

is the error function, and $a$ and $b$ are arbitrary constants. The error function lies in $(-1,1)$, is odd around $\eta =0$, and goes relatively quickly to $\pm 1$:

$$ \begin{align*} \lim_{\eta\rightarrow -\infty}\mbox{erf}(\eta) &=-1,\\ \lim_{\eta\rightarrow \infty}\mbox{erf}(\eta) &=1,\\ \mbox{erf}(\eta) &= -\mbox{erf}(-\eta),\\ \mbox{erf}(0) &=0,\\ \mbox{erf}(2) &=0.99532227,\\ \mbox{erf}(3) &=0.99997791 \thinspace . \end{align*} $$

As $t\rightarrow 0$, the error function approaches a step function centered at $x=c$. For a diffusion problem posed on the unit interval $[0,1]$, we may choose the step at $x=1/2$ (meaning $c=1/2$), $a=-1/2$, $b=1/2$. Then

$$ \begin{equation} u(x,t) = \frac{1}{2}\left(1 - \mbox{erf}\left(\frac{x-\frac{1}{2}}{\sqrt{4\dfc t}}\right)\right) = \frac{1}{2}\mbox{erfc}\left(\frac{x-\frac{1}{2}}{\sqrt{4\dfc t}}\right), \label{diffu:analysis:pde1:step:erf:sol} \tag{4} \end{equation} $$

where we have introduced the complementary error function $\mbox{erfc}(\eta) = 1-\mbox{erf}(\eta)$. The solution (4) implies the boundary conditions

$$ \begin{equation} u(0,t) = \frac{1}{2}\left(1 - \mbox{erf}\left(\frac{-1/2}{\sqrt{4\dfc t}}\right)\right), \label{diffu:analysis:pde1:p1:erf:uL} \tag{5} \end{equation} $$

$$ \begin{equation} u(1,t) = \frac{1}{2}\left(1 - \mbox{erf}\left(\frac{1/2}{\sqrt{4\dfc t}}\right)\right) \label{diffu:analysis:pde1:p1:erf:uR} \tag{6} \thinspace . \end{equation} $$

For small enough $t$, $u(0,t)\approx 1$ and $u(1,t)\approx 0$, but as $t\rightarrow\infty$, $u(x,t)\rightarrow 1/2$ on $[0,1]$.

Solution for a Gaussian pulse¶

The standard diffusion equation $u_t = \dfc u_{xx}$ admits a Gaussian function as solution:

$$ \begin{equation} u(x,t) = \frac{1}{\sqrt{4\pi\dfc t}} \exp{\left({-\frac{(x-c)^2}{4\dfc t}}\right)} \label{diffu:pde1:sol:Gaussian} \tag{7} \thinspace . \end{equation} $$

At $t=0$ this is a Dirac delta function, so for computational purposes one must start to view the solution at some time $t=t_\epsilon>0$. Replacing $t$ by $t_\epsilon +t$ in (7) makes it easy to operate with a (new) $t$ that starts at $t=0$ with an initial condition with a finite width. The important feature of (7) is that the standard deviation $\sigma$ of a sharp initial Gaussian pulse increases in time according to $\sigma = \sqrt{2\dfc t}$, making the pulse diffuse and flatten out.

Solution for a sine component¶

Also, (1) admits a solution of the form

$$ \begin{equation} u(x,t) = Qe^{-at}\sin\left( kx\right) \label{diffu:pde1:sol1} \tag{8} \thinspace . \end{equation} $$

The parameters $Q$ and $k$ can be freely chosen, while inserting (8) in (1) gives the constraint

$$ a = -\dfc k^2 \thinspace . $$

A very important feature is that the initial shape $I(x)=Q\sin\left( kx\right)$ undergoes a damping $\exp{(-\dfc k^2t)}$, meaning that rapid oscillations in space, corresponding to large $k$, are very much faster dampened than slow oscillations in space, corresponding to small $k$. This feature leads to a smoothing of the initial condition with time. (In fact, one can use a few steps of the diffusion equation as a method for removing noise in signal processing.) To judge how good a numerical method is, we may look at its ability to smoothen or dampen the solution in the same way as the PDE does.

The following example illustrates the damping properties of (8). We consider the specific problem

$$ \begin{align*} u_t &= u_{xx},\quad x\in (0,1),\ t\in (0,T],\\ u(0,t) &= u(1,t) = 0,\quad t\in (0,T],\\ u(x,0) & = \sin (\pi x) + 0.1\sin(100\pi x) \thinspace . \end{align*} $$

The initial condition has been chosen such that adding two solutions like (8) constructs an analytical solution to the problem:

$$ \begin{equation} u(x,t) = e^{-\pi^2 t}\sin (\pi x) + 0.1e^{-\pi^2 10^4 t}\sin (100\pi x) \label{diffu:pde1:sol2} \tag{9} \thinspace . \end{equation} $$

Figure illustrates the rapid damping of rapid oscillations $\sin (100\pi x)$ and the very much slower damping of the slowly varying $\sin (\pi x)$ term. After about $t=0.5\cdot10^{-4}$ the rapid oscillations do not have a visible amplitude, while we have to wait until $t\sim 0.5$ before the amplitude of the long wave $\sin (\pi x)$ becomes very small.

Evolution of the solution of a diffusion problem: initial condition (upper left), 1/100 reduction of the small waves (upper right), 1/10 reduction of the long wave (lower left), and 1/100 reduction of the long wave (lower right).

Analysis of discrete equations¶

A counterpart to (8) is the complex representation of the same function:

$$ u(x,t) = Qe^{-at}e^{ikx}, $$

where $i=\sqrt{-1}$ is the imaginary unit. We can add such functions, often referred to as wave components, to make a Fourier representation of a general solution of the diffusion equation:

$$ \begin{equation} u(x,t) \approx \sum_{k\in K} b_k e^{-\dfc k^2t}e^{ikx}, \label{diffu:pde1:u:Fourier} \tag{10} \end{equation} $$

where $K$ is a set of an infinite number of $k$ values needed to construct the solution. In practice, however, the series is truncated and $K$ is a finite set of $k$ values needed to build a good approximate solution. Note that (9) is a special case of (10) where $K=\{\pi, 100\pi\}$, $b_{\pi}=1$, and $b_{100\pi}=0.1$.

The amplitudes $b_k$ of the individual Fourier waves must be determined from the initial condition. At $t=0$ we have $u\approx\sum_kb_k\exp{(ikx)}$ and find $K$ and $b_k$ such that

$$ \begin{equation} I(x) \approx \sum_{k\in K} b_k e^{ikx}\thinspace . \label{_auto1} \tag{11} \end{equation} $$

(The relevant formulas for $b_k$ come from Fourier analysis, or equivalently, a least-squares method for approximating $I(x)$ in a function space with basis $\exp{(ikx)}$.)

Much insight about the behavior of numerical methods can be obtained by investigating how a wave component $\exp{(-\dfc k^2 t)}\exp{(ikx)}$ is treated by the numerical scheme. mathcal{I}_t appears that such wave components are also solutions of the schemes, but the damping factor $\exp{(-\dfc k^2 t)}$ varies among the schemes. To ease the forthcoming algebra, we write the damping factor as $A^n$. The exact amplification factor corresponding to $A$ is $\Aex = \exp{(-\dfc k^2\Delta t)}$.

Analysis of the finite difference schemes¶

We have seen that a general solution of the diffusion equation can be built as a linear combination of basic components

$$ e^{-\dfc k^2t}e^{ikx} \thinspace . $$

A fundamental question is whether such components are also solutions of the finite difference schemes. This is indeed the case, but the amplitude $\exp{(-\dfc k^2t)}$ might be modified (which also happens when solving the ODE counterpart $u'=-\dfc u$). We therefore look for numerical solutions of the form

$$ \begin{equation} u^n_q = A^n e^{ikq\Delta x} = A^ne^{ikx}, \label{diffu:pde1:analysis:uni} \tag{12} \end{equation} $$

where the amplification factor $A$ must be determined by inserting the component into an actual scheme. Note that $A^n$ means $A$ raised to the power of $n$, $n$ being the index in the time mesh, while the superscript $n$ in $u^n_q$ just denotes $u$ at time $t_n$.

Stability¶

The exact amplification factor is $\Aex=\exp{(-\dfc^2 k^2\Delta t)}$. We should therefore require $|A| < 1$ to have a decaying numerical solution as well. If $-1\leq A<0$, $A^n$ will change sign from time level to time level, and we get stable, non-physical oscillations in the numerical solutions that are not present in the exact solution.

Accuracy¶

To determine how accurately a finite difference scheme treats one wave component (12), we see that the basic deviation from the exact solution is reflected in how well $A^n$ approximates $\Aex^n$, or how well $A$ approximates $\Aex$. We can plot $\Aex$ and the various expressions for $A$, and we can make Taylor expansions of $A/\Aex$ to see the error more analytically.

Truncation error¶

As an alternative to examining the accuracy of the damping of a wave component, we can perform a general truncation error analysis as explained in "Truncation error analysis": "" [Langtangen_deqbook_trunc]. Such results are more general, but less detailed than what we get from the wave component analysis. The truncation error can almost always be computed and represents the error in the numerical model when the exact solution is substituted into the equations. In particular, the truncation error analysis tells the order of the scheme, which is of fundamental importance when verifying codes based on empirical estimation of convergence rates.

Analysis of the Forward Euler scheme¶

The Forward Euler finite difference scheme for $u_t = \dfc u_{xx}$ can be written as

$$ [D_t^+ u = \dfc D_xD_x u]^n_q\thinspace . $$

Inserting a wave component (12) in the scheme demands calculating the terms

$$ e^{ikq\Delta x}[D_t^+ A]^n = e^{ikq\Delta x}A^n\frac{A-1}{\Delta t}, $$

and

$$ A^nD_xD_x [e^{ikx}]_q = A^n\left( - e^{ikq\Delta x}\frac{4}{\Delta x^2} \sin^2\left(\frac{k\Delta x}{2}\right)\right) \thinspace . $$

Inserting these terms in the discrete equation and dividing by $A^n e^{ikq\Delta x}$ leads to

$$ \frac{A-1}{\Delta t} = -\dfc \frac{4}{\Delta x^2}\sin^2\left( \frac{k\Delta x}{2}\right), $$

and consequently

$$ \begin{equation} A = 1 -4F\sin^2 p \label{_auto2} \tag{13} \end{equation} $$

where

$$ \begin{equation} F = \frac{\dfc\Delta t}{\Delta x^2} \label{_auto3} \tag{14} \end{equation} $$

is the numerical Fourier number, and $p=k\Delta x/2$. The complete numerical solution is then

$$ \begin{equation} u^n_q = \left(1 -4F\sin^2 p\right)^ne^{ikq\Delta x} \thinspace . \label{_auto4} \tag{15} \end{equation} $$

Stability¶

We easily see that $A\leq 1$. However, the $A$ can be less than $-1$, which will lead to growth of a numerical wave component. The criterion $A\geq -1$ implies

$$ 4F\sin^2 (p/2)\leq 2 \thinspace . $$

The worst case is when $\sin^2 (p/2)=1$, so a sufficient criterion for stability is

$$ \begin{equation} F\leq {\frac{1}{2}}, \label{_auto5} \tag{16} \end{equation} $$

or expressed as a condition on $\Delta t$:

$$ \begin{equation} \Delta t\leq \frac{\Delta x^2}{2\dfc}\thinspace . \label{_auto6} \tag{17} \end{equation} $$

Note that halving the spatial mesh size, $\Delta x \rightarrow {\frac{1}{2}} \Delta x$, requires $\Delta t$ to be reduced by a factor of $1/4$. The method hence becomes very expensive for fine spatial meshes.

Accuracy¶

Since $A$ is expressed in terms of $F$ and the parameter we now call $p=k\Delta x/2$, we should also express $\Aex$ by $F$ and $p$. The exponent in $\Aex$ is $-\dfc k^2\Delta t$, which equals $-F k^2\Delta x^2=-F4p^2$. Consequently,

$$ \Aex = \exp{(-\dfc k^2\Delta t)} = \exp{(-4Fp^2)} \thinspace . $$

All our $A$ expressions as well as $\Aex$ are now functions of the two dimensionless parameters $F$ and $p$.

Computing the Taylor series expansion of $A/\Aex$ in terms of $F$ can easily be done with aid of sympy:

The result is

$$ \frac{A}{\Aex} = 1 - 4 F \sin^{2}p + 2F p^{2} - 16F^{2} p^{2} \sin^{2}p + 8 F^{2} p^{4} + \cdots $$

Recalling that $F=\dfc\Delta t/\Delta x^2$, $p=k\Delta x/2$, and that $\sin^2p\leq 1$, we realize that the dominating terms in $A/\Aex$ are at most

$$ 1 - 4\dfc \frac{\Delta t}{\Delta x^2} + \dfc\Delta t - 4\dfc^2\Delta t^2 + \dfc^2 \Delta t^2\Delta x^2 + \cdots \thinspace . $$

Truncation error¶

We follow the theory explained in "Truncation error analysis": "" [Langtangen_deqbook_trunc]. The recipe is to set up the scheme in operator notation and use formulas from "Overview of leading-order error terms in finite difference formulas": "" [Langtangen_deqbook_trunc] to derive an expression for the residual. The details are documented in "Linear diffusion equation in 1D": "" [Langtangen_deqbook_trunc]. We end up with a truncation error

$$ R^n_i = \Oof{\Delta t} + \Oof{\Delta x^2}\thinspace . $$

Although this is not the true error $\uex(x_i,t_n) - u^n_i$, it indicates that the true error is of the form

$$ E = C_t\Delta t + C_x\Delta x^2 $$

for two unknown constants $C_t$ and $C_x$.

Analysis of the Backward Euler scheme¶

Discretizing $u_t = \dfc u_{xx}$ by a Backward Euler scheme,

$$ [D_t^- u = \dfc D_xD_x u]^n_q, $$

and inserting a wave component (12), leads to calculations similar to those arising from the Forward Euler scheme, but since

$$ e^{ikq\Delta x}[D_t^- A]^n = A^ne^{ikq\Delta x}\frac{1 - A^{-1}}{\Delta t}, $$

we get

$$ \frac{1-A^{-1}}{\Delta t} = -\dfc \frac{4}{\Delta x^2}\sin^2\left( \frac{k\Delta x}{2}\right), $$

and then

$$ \begin{equation} A = \left(1 + 4F\sin^2p\right)^{-1} \label{diffu:pde1:analysis:BE:A} \tag{18} \thinspace . \end{equation} $$

The complete numerical solution can be written

$$ \begin{equation} u^n_q = \left(1 + 4F\sin^2 p\right)^{-n} e^{ikq\Delta x} \thinspace . \label{_auto7} \tag{19} \end{equation} $$

Stability¶

We see from (18) that $0<A<1$, which means that all numerical wave components are stable and non-oscillatory for any $\Delta t >0$.

Truncation error¶

The derivation of the truncation error for the Backward Euler scheme is almost identical to that for the Forward Euler scheme. We end up with

$$ R^n_i = \Oof{\Delta t} + \Oof{\Delta x^2}\thinspace . $$

Analysis of the Crank-Nicolson scheme¶

The Crank-Nicolson scheme can be written as

$$ [D_t u = \dfc D_xD_x \overline{u}^x]^{n+\frac{1}{2}}_q, $$

or

$$ [D_t u]^{n+\frac{1}{2}}_q = \frac{1}{2}\dfc\left( [D_xD_x u]^{n}_q + [D_xD_x u]^{n+1}_q\right) \thinspace . $$

Inserting (12) in the time derivative approximation leads to

$$ [D_t A^n e^{ikq\Delta x}]^{n+\frac{1}{2}} = A^{n+\frac{1}{2}} e^{ikq\Delta x}\frac{A^{\frac{1}{2}}-A^{-\frac{1}{2}}}{\Delta t} = A^ne^{ikq\Delta x}\frac{A-1}{\Delta t} \thinspace . $$

Inserting (12) in the other terms and dividing by $A^ne^{ikq\Delta x}$ gives the relation

$$ \frac{A-1}{\Delta t} = -\frac{1}{2}\dfc\frac{4}{\Delta x^2} \sin^2\left(\frac{k\Delta x}{2}\right) (1 + A), $$

and after some more algebra,

$$ \begin{equation} A = \frac{ 1 - 2F\sin^2p}{1 + 2F\sin^2p} \thinspace . \label{_auto8} \tag{20} \end{equation} $$

The exact numerical solution is hence

$$ \begin{equation} u^n_q = \left(\frac{ 1 - 2F\sin^2p}{1 + 2F\sin^2p}\right)^ne^{ikq\Delta x} \thinspace . \label{_auto9} \tag{21} \end{equation} $$

Stability¶

The criteria $A>-1$ and $A<1$ are fulfilled for any $\Delta t >0$. Therefore, the solution cannot grow, but it will oscillate if $1-2F\sin^p < 0$. To avoid such non-physical oscillations, we must demand $F\leq\frac{1}{2}$.

Truncation error¶

The truncation error is derived in "Linear diffusion equation in 1D": "" [Langtangen_deqbook_trunc]:

$$ R^{n+\frac{1}{2}}_i = \Oof{\Delta x^2} + \Oof{\Delta t^2}\thinspace . $$

Analysis of the Leapfrog scheme¶

An attractive feature of the Forward Euler scheme is the explicit time stepping and no need for solving linear systems. However, the accuracy in time is only $\Oof{\Delta t}$. We can get an explicit second-order scheme in time by using the Leapfrog method:

$$ [D_{2t} u = \dfc D_xDx u + f]^n_q\thinspace . $$

Written out,

$$ u_q^{n+1} = u_q^{n-1} + \frac{2\dfc\Delta t}{\Delta x^2} (u^{n}_{q+1} - 2u^n_q + u^n_{q-1}) + f(x_q,t_n)\thinspace . $$

We need some formula for the first step, $u^1_q$, but for that we can use a Forward Euler step.

Unfortunately, the Leapfrog scheme is always unstable for the diffusion equation. To see this, we insert a wave component $A^ne^{ikx}$ and get

$$ \frac{A - A^{-1}}{\Delta t} = -\dfc \frac{4}{\Delta x^2}\sin^2 p, $$

or

$$ A^2 + 4F \sin^2 p\, A - 1 = 0, $$

which has roots

$$ A = -2F\sin^2 p \pm \sqrt{4F^2\sin^4 p + 1}\thinspace . $$

Both roots have $|A|>1$ so the amplitude always grows, which is not in accordance with the physics of the problem. However, for a PDE with a first-order derivative in space, instead of a second-order one, the Leapfrog scheme performs very well.

Summary of accuracy of amplification factors¶

We can plot the various amplification factors against $p=k\Delta x/2$ for different choices of the $F$ parameter. Figures diffu:pde1:fig:A:err:C20, diffu:pde1:fig:A:err:C0.5, and diffu:pde1:fig:A:err:C0.1 show how long and small waves are damped by the various schemes compared to the exact damping. As long as all schemes are stable, the amplification factor is positive, except for Crank-Nicolson when $F>0.5$.

Amplification factors for large time steps.

Amplification factors for time steps around the Forward Euler stability limit.

Amplification factors for small time steps.

The effect of negative amplification factors is that $A^n$ changes sign from one time level to the next, thereby giving rise to oscillations in time in an animation of the solution. We see from Figure that for $F=20$, waves with $p\geq \pi/4$ undergo a damping close to $-1$, which means that the amplitude does not decay and that the wave component jumps up and down (flips amplitude) in time. For $F=2$ we have a damping of a factor of 0.5 from one time level to the next, which is very much smaller than the exact damping. Short waves will therefore fail to be effectively dampened. These waves will manifest themselves as high frequency oscillatory noise in the solution.

A value $p=\pi/4$ corresponds to four mesh points per wave length of $e^{ikx}$, while $p=\pi/2$ implies only two points per wave length, which is the smallest number of points we can have to represent the wave on the mesh.

To demonstrate the oscillatory behavior of the Crank-Nicolson scheme, we choose an initial condition that leads to short waves with significant amplitude. A discontinuous $I(x)$ will in particular serve this purpose: Figures diffu:pde1:CN:fig:F=3 and diffu:pde1:CN:fig:F=10 correspond to $F=3$ and $F=10$, respectively, and we see how short waves pollute the overall solution.

Analysis of the 2D diffusion equation¶

Diffusion in several dimensions is treated later, but it is appropriate to include the analysis here. We first consider the 2D diffusion equation

$$ u_{t} = \dfc(u_{xx} + u_{yy}), $$

which has Fourier component solutions of the form

$$ u(x,y,t) = Ae^{-\dfc k^2t}e^{i(k_x x + k_yy)}, $$

and the schemes have discrete versions of this Fourier component:

$$ u^{n}_{q,r} = A\xi^{n}e^{i(k_x q\Delta x + k_y r\Delta y)}\thinspace . $$

The Forward Euler scheme¶

For the Forward Euler discretization,

$$ [D_t^+u = \dfc(D_xD_x u + D_yD_y u)]_{q,r}^n, $$

we get

$$ \frac{\xi - 1}{\Delta t} = -\dfc\frac{4}{\Delta x^2}\sin^2\left(\frac{k_x\Delta x}{2}\right) - \dfc\frac{4}{\Delta y^2}\sin^2\left(\frac{k_y\Delta y}{2}\right)\thinspace . $$

Introducing

$$ p_x = \frac{k_x\Delta x}{2},\quad p_y = \frac{k_y\Delta y}{2}, $$

we can write the equation for $\xi$ more compactly as

$$ \frac{\xi - 1}{\Delta t} = -\dfc\frac{4}{\Delta x^2}\sin^2 p_x - \dfc\frac{4}{\Delta y^2}\sin^2 p_y, $$

and solve for $\xi$:

$$ \begin{equation} \xi = 1 - 4F_x\sin^2 p_x - 4F_y\sin^2 p_y\thinspace . \label{diffu:2D:analysis:xi} \tag{22} \end{equation} $$

The complete numerical solution for a wave component is

$$ \begin{equation} u^{n}_{q,r} = A(1 - 4F_x\sin^2 p_x - 4F_y\sin^2 p_y)^n e^{i(k_xq\Delta x + k_yr\Delta y)}\thinspace . \label{diffu:2D:analysis:FE:numexact} \tag{23} \end{equation} $$

For stability we demand $-1\leq\xi\leq 1$, and $-1\leq\xi$ is the critical limit, since clearly $\xi \leq 1$, and the worst case happens when the sines are at their maximum. The stability criterion becomes

$$ \begin{equation} F_x + F_y \leq \frac{1}{2}\thinspace . \label{diffu:2D:analysis:FE:stab} \tag{24} \end{equation} $$

For the special, yet common, case $\Delta x=\Delta y=h$, the stability criterion can be written as

$$ \Delta t \leq \frac{h^2}{2d\dfc}, $$

where $d$ is the number of space dimensions: $d=1,2,3$.

The Backward Euler scheme¶

The Backward Euler method,

$$ [D_t^-u = \dfc(D_xD_x u + D_yD_y u)]_{q,r}^n, $$

results in

$$ 1 - \xi^{-1} = - 4F_x \sin^2 p_x - 4F_y \sin^2 p_y, $$

and

$$ \xi = (1 + 4F_x \sin^2 p_x + 4F_y \sin^2 p_y)^{-1}, $$

which is always in $(0,1]$. The solution for a wave component becomes

$$ \begin{equation} u^{n}_{q,r} = A(1 + 4F_x\sin^2 p_x + 4F_y\sin^2 p_y)^{-n} e^{i(k_xq\Delta x + k_yr\Delta y)}\thinspace . \label{diffu:2D:analysis:BN:numexact} \tag{25} \end{equation} $$

The Crank-Nicolson scheme¶

With a Crank-Nicolson discretization,

$$ [D_tu]^{n+\frac{1}{2}}_{q,r} = \frac{1}{2} [\dfc(D_xD_x u + D_yD_y u)]_{q,r}^{n+1} + \frac{1}{2} [\dfc(D_xD_x u + D_yD_y u)]_{q,r}^n, $$

we have, after some algebra,

$$ \xi = \frac{1 - 2(F_x\sin^2 p_x + F_x\sin^2p_y)}{1 + 2(F_x\sin^2 p_x + F_x\sin^2p_y)}\thinspace . $$

The fraction on the right-hand side is always less than 1, so stability in the sense of non-growing wave components is guaranteed for all physical and numerical parameters. However, the fraction can become negative and result in non-physical oscillations. This phenomenon happens when

$$ F_x\sin^2 p_x + F_x\sin^2p_y > \frac{1}{2}\thinspace . $$

A criterion against non-physical oscillations is therefore

$$ F_x + F_y \leq \frac{1}{2}, $$

which is the same limit as the stability criterion for the Forward Euler scheme.

The exact discrete solution is

$$ \begin{equation} u^{n}_{q,r} = A \left( \frac{1 - 2(F_x\sin^2 p_x + F_x\sin^2p_y)}{1 + 2(F_x\sin^2 p_x + F_x\sin^2p_y)} \right)^n e^{i(k_xq\Delta x + k_yr\Delta y)}\thinspace . \label{diffu:2D:analysis:CN:numexact} \tag{26} \end{equation} $$

Explanation of numerical artifacts¶

The behavior of the solution generated by Forward Euler discretization in time (and centered differences in space) is summarized at the end of the section diffu:pde1:FE:experiments. Can we, from the analysis above, explain the behavior?

We may start by looking at Figure where $F=0.51$. The figure shows that the solution is unstable and grows in time. The stability limit for such growth is $F=0.5$ and since the $F$ in this simulation is slightly larger, growth is unavoidable.

Figure has unexpected features: we would expect the solution of the diffusion equation to be smooth, but the graphs in Figure contain non-smooth noise. Turning to Figure, which has a quite similar initial condition, we see that the curves are indeed smooth. The problem with the results in Figure is that the initial condition is discontinuous. To represent it, we need a significant amplitude on the shortest waves in the mesh. However, for $F=0.5$, the shortest wave ($p=\pi/2$) gives the amplitude in the numerical solution as $(1-4F)^n$, which oscillates between negative and positive values at subsequent time levels for $F>\frac{1}{4}$. Since the shortest waves have visible amplitudes in the solution profile, the oscillations becomes visible. The smooth initial condition in Figure, on the other hand, leads to very small amplitudes of the shortest waves. That these waves then oscillate in a non-physical way for $F=0.5$ is not a visible effect. The oscillations in time in the amplitude $(1-4F)^n$ disappear for $F\leq\frac{1}{4}$, and that is why also the discontinuous initial condition always leads to smooth solutions in Figure, where $F=\frac{1}{4}$.

Turning the attention to the Backward Euler scheme and the experiments in Figure, we see that even the discontinuous initial condition gives smooth solutions for $F=0.5$ (and in fact all other $F$ values). From the exact expression of the numerical amplitude, $(1 + 4F\sin^2p)^{-1}$, we realize that this factor can never flip between positive and negative values, and no instabilities can occur. The conclusion is that the Backward Euler scheme always produces smooth solutions. Also, the Backward Euler scheme guarantees that the solution cannot grow in time (unless we add a source term to the PDE, but that is meant to represent a physically relevant growth).

Finally, we have some small, strange artifacts when simulating the development of the initial plug profile with the Crank-Nicolson scheme, see Figure, where $F=3$. The Crank-Nicolson scheme cannot give growing amplitudes, but it may give oscillating amplitudes in time. The critical factor is $1 - 2F\sin^2p$, which for the shortest waves ($p=\pi/2$) indicates a stability limit $F=0.5$. With the discontinuous initial condition, we have enough amplitude on the shortest waves so their wrong behavior is visible, and this is what we see as small instabilities in Figure. The only remedy is to lower the $F$ value.

Exercises¶

Exercise 1: Explore symmetry in a 1D problem¶

This exercise simulates the exact solution (7). Suppose for simplicity that $c=0$.

a) Formulate an initial-boundary value problem that has (7) as solution in the domain $[-L,L]$. Use the exact solution (7) as Dirichlet condition at the boundaries. Simulate the diffusion of the Gaussian peak. Observe that the solution is symmetric around $x=0$.

b) Show from (7) that $u_x(c,t)=0$. Since the solution is symmetric around $x=c=0$, we can solve the numerical problem in frac{1}{2} of the domain, using a symmetry boundary condition $u_x=0$ at $x=0$. Set up the initial-boundary value problem in this case. Simulate the diffusion problem in $[0,L]$ and compare with the solution in a).

Solution.

$$ \begin{align*} u_t &= \dfc u_xx,\\ u_x(0,t) &= 0,\\ u(L,t)& =\frac{1}{\sqrt{4\pi\dfc t}} \exp{\left({-\frac{x^2}{4\dfc t}}\right)}\thinspace . \end{align*} $$

Filename: diffu_symmetric_gaussian.

Exercise 2: Investigate approximation errors from a $u_x=0$ boundary condition¶

We consider the problem solved in Exercise 1: Explore symmetry in a 1D problem part b). The boundary condition $u_x(0,t)=0$ can be implemented in two ways: 1) by a standard symmetric finite difference $[D_{2x}u]_i^n=0$, or 2) by a one-sided difference $[D^+u=0]^n_i=0$. Investigate the effect of these two conditions on the convergence rate in space.

Hint. If you use a Forward Euler scheme, choose a discretization parameter $h=\Delta t = \Delta x^2$ and assume the error goes like $E\sim h^r$. The error in the scheme is $\Oof{\Delta t,\Delta x^2}$ so one should expect that the estimated $r$ approaches 1. The question is if a one-sided difference approximation to $u_x(0,t)=0$ destroys this convergence rate.

Filename: diffu_onesided_fd.

Exercise 3: Experiment with open boundary conditions in 1D¶

We address diffusion of a Gaussian function as in Exercise 1: Explore symmetry in a 1D problem, in the domain $[0,L]$, but now we shall explore different types of boundary conditions on $x=L$. In real-life problems we do not know the exact solution on $x=L$ and must use something simpler.

a) Imagine that we want to solve the problem numerically on $[0,L]$, with a symmetry boundary condition $u_x=0$ at $x=0$, but we do not know the exact solution and cannot of that reason assign a correct Dirichlet condition at $x=L$. One idea is to simply set $u(L,t)=0$ since this will be an accurate approximation before the diffused pulse reaches $x=L$ and even thereafter it might be a satisfactory condition if the exact $u$ has a small value. Let $\uex$ be the exact solution and let $u$ be the solution of $u_t=\dfc u_{xx}$ with an initial Gaussian pulse and the boundary conditions $u_x(0,t)=u(L,t)=0$. Derive a diffusion problem for the error $e=\uex - u$. Solve this problem numerically using an exact Dirichlet condition at $x=L$. Animate the evolution of the error and make a curve plot of the error measure

$$ E(t)=\sqrt{\frac{\int_0^L e^2dx}{\int_0^L udx}}\thinspace . $$

Is this a suitable error measure for the present problem?

b) Instead of using $u(L,t)=0$ as approximate boundary condition for letting the diffused Gaussian pulse move out of our finite domain, one may try $u_x(L,t)=0$ since the solution for large $t$ is quite flat. Argue that this condition gives a completely wrong asymptotic solution as $t\rightarrow 0$. To do this, integrate the diffusion equation from $0$ to $L$, integrate $u_{xx}$ by parts (or use Gauss' divergence theorem in 1D) to arrive at the important property

$$ \frac{d}{dt}\int_{0}^L u(x,t)dx = 0, $$

implying that $\int_0^Ludx$ must be constant in time, and therefore

$$ \int_{0}^L u(x,t)dx = \int_{0}^LI(x)dx\thinspace . $$

The integral of the initial pulse is 1.

c) Another idea for an artificial boundary condition at $x=L$ is to use a cooling law

$$ \begin{equation} -\dfc u_x = q(u - u_S), \label{diffu:pde1:Gaussian:xL:cooling} \tag{27} \end{equation} $$

where $q$ is an unknown heat transfer coefficient and $u_S$ is the surrounding temperature in the medium outside of $[0,L]$. (Note that arguing that $u_S$ is approximately $u(L,t)$ gives the $u_x=0$ condition from the previous subexercise that is qualitatively wrong for large $t$.) Develop a diffusion problem for the error in the solution using (27) as boundary condition. Assume one can take $u_S=0$ "outside the domain" since $\uex\rightarrow 0$ as $x\rightarrow\infty$. Find a function $q=q(t)$ such that the exact solution obeys the condition (27). Test some constant values of $q$ and animate how the corresponding error function behaves. Also compute $E(t)$ curves as defined above.

Filename: diffu_open_BC.

Exercise 4: Simulate a diffused Gaussian peak in 2D/3D¶

a) Generalize (7) to multi dimensions by assuming that one-dimensional solutions can be multiplied to solve $u_t = \dfc\nabla^2 u$. Set $c=0$ such that the peak of the Gaussian is at the origin.

b) One can from the exact solution show that $u_x=0$ on $x=0$, $u_y=0$ on $y=0$, and $u_z=0$ on $z=0$. The approximately correct condition $u=0$ can be set on the remaining boundaries (say $x=L$, $y=L$, $z=L$), cf. Exercise 3: Experiment with open boundary conditions in 1D. Simulate a 2D case and make an animation of the diffused Gaussian peak.

c) The formulation in b) makes use of symmetry of the solution such that we can solve the problem in the first quadrant (2D) or octant (3D) only. To check that the symmetry assumption is correct, formulate the problem without symmetry in a domain $[-L,L]\times [L,L]$ in 2D. Use $u=0$ as approximately correct boundary condition. Simulate the same case as in b), but in a four times as large domain. Make an animation and compare it with the one in b).

Filename: diffu_symmetric_gaussian_2D.

Exercise 5: Examine stability of a diffusion model with a source term¶

Consider a diffusion equation with a linear $u$ term:

$$ u_t = \dfc u_{xx} + \beta u\thinspace . $$

a) Derive in detail the Forward Euler, Backward Euler, and Crank-Nicolson schemes for this type of diffusion model. Thereafter, formulate a $\theta$-rule to summarize the three schemes.

b) Assume a solution like (8) and find the relation between $a$, $k$, $\dfc$, and $\beta$.

Hint. Insert (8) in the PDE problem.

c) Calculate the stability of the Forward Euler scheme. Design numerical experiments to confirm the results.

Hint. Insert the discrete counterpart to (8) in the numerical scheme. Run experiments at the stability limit and slightly above.

d) Repeat c) for the Backward Euler scheme.

e) Repeat c) for the Crank-Nicolson scheme.

f) How does the extra term $bu$ impact the accuracy of the three schemes?

Hint. For analysis of the accuracy, compare the numerical and exact amplification factors, in graphs and/or by Taylor series expansion.

Filename: diffu_stability_uterm.