Truncation errors in diffusion equations¶
Linear diffusion equation in 1D¶
The standard, linear, 1D diffusion equation takes the form
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:
The truncation error arises as the residual \(R\) when inserting the exact solution \({u_{\small\mbox{e}}}\) in the discrete equations:
Now, using (11)-(12) and (17)-(18), we can transform the difference operators to derivatives:
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
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:
The equation for the truncation error is
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 (21)-(22),
With (17)-(18) we can express the difference operator \(D_xD_xu\) in terms of a derivative:
The error term from the arithmetic mean is similarly expanded,
The time derivative is analyzed using (5)-(6):
Summing up all the contributions and notifying that
the truncation error is given by
Linear diffusion equation in 2D/3D¶
A nonlinear diffusion equation in 2D¶
Exercises¶
Exercise 1: Truncation error of a weighted mean¶
Derive the truncation error of the weighted mean in (19)-(20).
Hint. Expand \({u_{\small\mbox{e}}}^{n+1}\) and \({u_{\small\mbox{e}}}^n\) around \(t_{n+\theta}\).
Filename: trunc_weighted_mean.
Exercise 2: Simulate the error of a weighted mean¶
We consider the weighted mean
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
(19)-(20).
Filename: trunc_theta_avg.
Exercise 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
(15)-(16).
Filename: trunc_backward_2level.
Exercise 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 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 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.
Exercise 7: Verify the effect of correction terms¶
The program decay_convrate.py
solves \(u'=-au\), \(u(0)=I\), by the \(\theta\)-rule and computes
convergence rates. Copy this file and
adjust \(a\) in the solver function 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
\({\mathcal{O}(\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.
Exercise 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}}\):
Compute the truncation error in both cases.
Filename: trunc_decay_CN_vc.
Exercise 9: Truncation error of \(u'=f(u,t)\)¶
Consider the general nonlinear first-order scalar ODE
Show that the truncation error in the Forward Euler scheme,
and in the Backward Euler scheme,
both are of first order, regardless of what \(f\) is.
Showing the order of the truncation error in the Crank-Nicolson scheme,
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
Check that the derived truncation error is consistent with previous
results for the case \(f(u,t)=-au\).
Filename: trunc_nonlinear_ODE.
Exercise 10: Truncation error of \([D_t D_tu]^n\)¶
Derive the truncation error of the finite difference approximation
(17)-(18) to
the second-order derivative.
Filename: trunc_d2u.
Exercise 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.
Exercise 12: Investigate the accuracy of a simplified scheme¶
Consider the ODE
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
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.
