Diffusion equations

Linear diffusion equation in 1D

The standard, linear, 1D diffusion equation takes the form

$$\tag{738} \frac{\partial u}{\partial t} = {\alpha}\frac{\partial^2 u}{\partial x^2} + f(x,t),\quad x\in (0, L),\ t\in (0,T],$$

where ${\alpha} > 0$ is the constant diffusion coefficient. A more compact form of the diffusion equation is $u_t = {\alpha} u_{xx}+f$.

The spatial derivative in the diffusion equation, ${\alpha} u_{xx}$, is commonly discretized as $[D_x D_xu]^n_i$. The time-derivative, however, can be treated by a variety of methods.

The Forward Euler scheme in time

Let us start with the simple Forward Euler scheme:

$$[D_t^+ u = {\alpha} D_xD_x u + f]^n{\thinspace .}$$

The truncation error arises as the residual $R$ when inserting the exact solution ${u_{\small\mbox{e}}}$ in the discrete equations:

$$[D_t^+ {u_{\small\mbox{e}}} = {\alpha} D_xD_x {u_{\small\mbox{e}}} + f + R]^n_i{\thinspace .}$$

Now, using (664)-(665) and (670)-(671), we can transform the difference operators to derivatives:

$$\begin{split}\begin{align*} {u_{\small\mbox{e}, t}}(x_i,t_n) &+ \frac{1}{2}{u_{\small\mbox{e}, tt}}(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} = {\alpha}{u_{\small\mbox{e}, xx}}(x_i,t_n) + \\ &\frac{{\alpha}}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)} + f(x_i,t_n) + R^n_i{\thinspace .} \end{align*}\end{split}$$

The terms ${u_{\small\mbox{e}, t}}(x_i,t_n) - {\alpha}{u_{\small\mbox{e}, xx}}(x_i,t_n) - f(x_i,t_n)$ vanish because ${u_{\small\mbox{e}}}$ solves the PDE. The truncation error then becomes

$$R^n_i = \frac{1}{2}{u_{\small\mbox{e}, tt}}(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} - \frac{{\alpha}}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)}{\thinspace .}$$

The Crank-Nicolson scheme in time

The Crank-Nicolson method consists of using a centered difference for $u_t$ and an arithmetic average of the $u_{xx}$ term:

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

The equation for the truncation error is

$$[D_t {u_{\small\mbox{e}}}]^{n+\frac{1}{2}}_i = {\alpha}\frac{1}{2}([D_xD_x {u_{\small\mbox{e}}}]^n_i + [D_xD_x {u_{\small\mbox{e}}}]^{n+1}_i) + f^{n+\frac{1}{2}}_i + R^{n+\frac{1}{2}}_i{\thinspace .}$$

To find the truncation error, we start by expressing the arithmetic average in terms of values at time $t_{n+\frac{1}{2}}$. According to (674)-(675),

$$\frac{1}{2}([D_xD_x {u_{\small\mbox{e}}}]^n_i + [D_xD_x {u_{\small\mbox{e}}}]^{n+1}_i) = [D_xD_x{u_{\small\mbox{e}}}]^{n+\frac{1}{2}}_i + \frac{1}{8}[D_xD_x{u_{\small\mbox{e}, tt}}]_i^{n+\frac{1}{2}}\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .}$$

With (670)-(671) we can express the difference operator $D_xD_xu$ in terms of a derivative:

$$[D_xD_x{u_{\small\mbox{e}}}]^{n+\frac{1}{2}}_i = {u_{\small\mbox{e}, xx}}(x_i, t_{n+\frac{1}{2}}) + \frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i, t_{n+\frac{1}{2}})\Delta x^2 + {\mathcal{O}(\Delta x^4)}{\thinspace .}$$

The error term from the arithmetic mean is similarly expanded,

$$\frac{1}{8}[D_xD_x{u_{\small\mbox{e}, tt}}]_i^{n+\frac{1}{2}}\Delta t^2 = \frac{1}{8}{u_{\small\mbox{e}, ttxx}}(x_i, t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^2\Delta x^2)}$$

The time derivative is analyzed using (658)-(659):

$$[D_t u]^{n+\frac{1}{2}}_i = {u_{\small\mbox{e}, t}}(x_i,t_{n+\frac{1}{2}}) + \frac{1}{24}{u_{\small\mbox{e}, ttt}}(x_i,t_{n+\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .}$$

Summing up all the contributions and notifying that

$${u_{\small\mbox{e}, t}}(x_i,t_{n+\frac{1}{2}}) = {\alpha}{u_{\small\mbox{e}, xx}}(x_i, t_{n+\frac{1}{2}}) + f(x_i,t_{n+\frac{1}{2}}),$$

the truncation error is given by

$$\begin{split}\begin{align*} R^{n+\frac{1}{2}}_i & = \frac{1}{8}{u_{\small\mbox{e}, xx}}(x_i,t_{n+\frac{1}{2}})\Delta t^2 + \frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i, t_{n+\frac{1}{2}})\Delta x^2 +\\ &\quad \frac{1}{24}{u_{\small\mbox{e}, ttt}}(x_i,t_{n+\frac{1}{2}})\Delta t^2 + + {\mathcal{O}(\Delta x^4)} + {\mathcal{O}(\Delta t^4)} + {\mathcal{O}(\Delta t^2\Delta x^2)} \end{align*}\end{split}$$

Nonlinear diffusion equation in 1D

We address the PDE

$$\frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\left( {\alpha}(u)\frac{\partial u}{\partial x}\right) + f(u),$$

with two potentially nonlinear coefficients $q(u)$ and ${\alpha}(u)$. We use a Backward Euler scheme with arithmetic mean for ${\alpha}(u)$,

$$[D^-u = D_x\overline{{\alpha}(u)}^{x}D_x u + f(u)]_i^n{\thinspace .}$$

Inserting ${u_{\small\mbox{e}}}$ defines the truncation error $R$:

$$[D^-{u_{\small\mbox{e}}} = D_x\overline{{\alpha}({u_{\small\mbox{e}}})}^{x}D_x {u_{\small\mbox{e}}} + f({u_{\small\mbox{e}}})]_i^n{\thinspace .}$$

The most computationally challenging part is the variable coefficient with ${\alpha}(u)$, but we can use the same setup as in the section Extension to variable coefficients and arrive at a truncation error ${\mathcal{O}(\Delta x^2)}$ for the $x$-derivative term. The nonlinear term $[f({u_{\small\mbox{e}}})]=^n_{i} = f({u_{\small\mbox{e}}}(x_i, t_n))$ matches $x$ and $t$ derivatives of ${u_{\small\mbox{e}}}$ in the PDE. We end up with

$$R^n_i = -{\frac{1}{2}}\frac{\partial^2}{\partial t^2}{u_{\small\mbox{e}}}(x_i,t_n)\Delta t + {\mathcal{O}(\Delta x^2)}{\thinspace .}$$

Exercises

Exercise B.1: Truncation error of a weighted mean

Derive the truncation error of the weighted mean in (672)-(673).

Hint. Expand ${u_{\small\mbox{e}}}^{n+1}$ and ${u_{\small\mbox{e}}}^n$ around $t_{n+\theta}$.

Filename: trunc_weighted_mean.

Exercise B.2: Simulate the error of a weighted mean

We consider the weighted mean

$${u_{\small\mbox{e}}}(t_n) \approx \theta {u_{\small\mbox{e}}}^{n+1} + (1-\theta){u_{\small\mbox{e}}}^n{\thinspace .}$$

Choose some specific function for ${u_{\small\mbox{e}}}(t)$ and compute the error in this approximation for a sequence of decreasing $\Delta t = t_{n+1}-t_n$ and for $\theta = 0, 0.25, 0.5, 0.75, 1$. Assuming that the error equals $C\Delta t^r$, for some constants $C$ and $r$, compute $r$ for the two smallest $\Delta t$ values for each choice of $\theta$ and compare with the truncation error (672)-(673). Filename: trunc_theta_avg.

Exercise B.3: Verify a truncation error formula

Set up a numerical experiment as explained in the section Empirical verification of the truncation error for verifying the formulas (668)-(669). Filename: trunc_backward_2level.

Problem B.4: Truncation error of the Backward Euler scheme

Derive the truncation error of the Backward Euler scheme for the decay ODE $u'=-au$ with constant $a$. Extend the analysis to cover the variable-coefficient case $u'=-a(t)u + b(t)$. Filename: trunc_decay_BE.

Exercise B.5: Empirical estimation of truncation errors

Use the ideas and tools from the section Empirical verification of the truncation error to estimate the rate of the truncation error of the Backward Euler and Crank-Nicolson schemes applied to the exponential decay model $u'=-au$, $u(0)=I$.

Hint. In the Backward Euler scheme, the truncation error can be estimated at mesh points $n=1,\ldots,N$, while the truncation error must be estimated at midpoints $t_{n+\frac{1}{2}}$, $n=0,\ldots,N-1$ for the Crank-Nicolson scheme. The truncation_error(dt, N) function to be supplied to the estimate function needs to carefully implement these details and return the right t array such that t[i] is the time point corresponding to the quantities R[i] and R_a[i].

Filename: trunc_decay_BNCN.

Exercise B.6: Correction term for a Backward Euler scheme

Consider the model $u'=-au$, $u(0)=I$. Use the ideas of the section Increasing the accuracy by adding correction terms to add a correction term to the ODE such that the Backward Euler scheme applied to the perturbed ODE problem is of second order in $\Delta t$. Find the amplification factor. Filename: trunc_decay_BE_corr.

Problem B.7: Verify the effect of correction terms

Make a program that solves $u'=-au$, $u(0)=I$, by the $\theta$-rule and computes convergence rates. Adjust $a$ such that it incorporates correction terms. Run the program to verify that the error from the Forward and Backward Euler schemes with perturbed $a$ is $\Oof{\Delta t^2}$, while the error arising from the Crank-Nicolson scheme with perturbed $a$ is ${\mathcal{O}(\Delta t^4)}$. Filename: trunc_decay_corr_verify.

Problem B.8: Truncation error of the Crank-Nicolson scheme

The variable-coefficient ODE $u'=-a(t)u+b(t)$ can be discretized in two different ways by the Crank-Nicolson scheme, depending on whether we use averages for $a$ and $b$ or compute them at the midpoint $t_{n+\frac{1}{2}}$:

$$\tag{739} \lbrack D_t u = -a\overline{u}^t + b \rbrack^{n+\frac{1}{2}},$$

$$\tag{740} \lbrack D_t u = \overline{-au+b}^t \rbrack^{n+\frac{1}{2}} {\thinspace .}$$

Compute the truncation error in both cases. Filename: trunc_decay_CN_vc.

Problem B.9: Truncation error of $u'=f(u,t)$

Consider the general nonlinear first-order scalar ODE

$$u'(t) = f(u(t), t) {\thinspace .}$$

Show that the truncation error in the Forward Euler scheme,

$$[D_t^+ u = f(u,t)]^n,$$

and in the Backward Euler scheme,

$$[D_t^- u = f(u,t)]^n,$$

both are of first order, regardless of what $f$ is.

Showing the order of the truncation error in the Crank-Nicolson scheme,

$$[D_t u = f(u,t)]^{n+\frac{1}{2}},$$

is somewhat more involved: Taylor expand ${u_{\small\mbox{e}}}^n$, ${u_{\small\mbox{e}}}^{n+1}$, $f({u_{\small\mbox{e}}}^n, t_n)$, and $f({u_{\small\mbox{e}}}^{n+1}, t_{n+1})$ around $t_{n+\frac{1}{2}}$, and use that

$$\frac{df}{dt} = \frac{\partial f}{\partial u}u' + \frac{\partial f}{\partial t} {\thinspace .}$$

Check that the derived truncation error is consistent with previous results for the case $f(u,t)=-au$. Filename: trunc_nonlinear_ODE.

Exercise B.10: Truncation error of $[D_t D_tu]^n$

Derive the truncation error of the finite difference approximation (670)-(671) to the second-order derivative. Filename: trunc_d2u.

Exercise B.11: Investigate the impact of approximating $u'(0)$

The section Linear model without damping describes two ways of discretizing the initial condition $u'(0)=V$ for a vibration model $u''+\omega^2u=0$: a centered difference $[D_{2t}u=V]^0$ or a forward difference $[D_t^+u=V]^0$. The program vib_undamped.py solves $u''+\omega^2u=0$ with $[D_{2t}u=0]^0$ and features a function convergence_rates for computing the order of the error in the numerical solution. Modify this program such that it applies the forward difference $[D_t^+u=0]^0$ and report how this simpler and more convenient approximation impacts the overall convergence rate of the scheme. Filename: trunc_vib_ic_fw.

Problem B.12: Investigate the accuracy of a simplified scheme

Consider the ODE

$$mu'' + \beta |u'|u' + s(u) = F(t){\thinspace .}$$

The term $|u'|u'$ quickly gives rise to nonlinearities and complicates the scheme. Why not simply apply a backward difference to this term such that it only involves known values? That is, we propose to solve

$$[mD_tD_tu + \beta |D_t^-u|D_t^-u + s(u) = F]^n{\thinspace .}$$

Drop the absolute value for simplicity and find the truncation error of the scheme. Perform numerical experiments with the scheme and compared with the one based on centered differences. Can you illustrate the accuracy loss visually in real computations, or is the asymptotic analysis here mainly of theoretical interest? Filename: trunc_vib_bw_damping.