The standard, linear, 1D diffusion equation takes the form $$ \begin{equation} \frac{\partial u}{\partial t} = \dfc\frac{\partial^2 u}{\partial x^2} + f(x,t),\quad x\in (0, L),\ t\in (0,T], \tag{730} \end{equation} $$ where \( \dfc > 0 \) is the constant diffusion coefficient. A more compact form of the diffusion equation is \( u_t = \dfc u_{xx}+f \).
The spatial derivative in the diffusion equation, \( \dfc u_{xx} \), is commonly discretized as \( [D_x D_xu]^n_i \). The time-derivative, however, can be treated by a variety of methods.
Let us start with the simple Forward Euler scheme: $$ [D_t^+ u = \dfc D_xD_x u + f]^n\tp$$ The truncation error arises as the residual \( R \) when inserting the exact solution \( \uex \) in the discrete equations: $$ [D_t^+ \uex = \dfc D_xD_x \uex + f + R]^n_i\tp$$ Now, using (656)-(657) and (662)-(663), we can transform the difference operators to derivatives: $$ \begin{align*} \uexd{t}(x_i,t_n) &+ \half\uexd{tt}(t_n)\Delta t + \Oof{\Delta t^2} = \dfc\uexd{xx}(x_i,t_n) + \\ &\frac{\dfc}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2 + \Oof{\Delta x^4} + f(x_i,t_n) + R^n_i\tp \end{align*} $$ The terms \( \uexd{t}(x_i,t_n) - \dfc\uexd{xx}(x_i,t_n) - f(x_i,t_n) \) vanish because \( \uex \) solves the PDE. The truncation error then becomes $$ R^n_i = \half\uexd{tt}(t_n)\Delta t + \Oof{\Delta t^2} - \frac{\dfc}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2 + \Oof{\Delta x^4}\tp $$
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+\half}_i = \dfc\half([D_xD_x u]^n_i + [D_xD_x u]^{n+1}_i + f^{n+\half}_i\tp$$ The equation for the truncation error is $$ [D_t \uex]^{n+\half}_i = \dfc\half([D_xD_x \uex]^n_i + [D_xD_x \uex]^{n+1}_i) + f^{n+\half}_i + R^{n+\half}_i\tp$$ To find the truncation error, we start by expressing the arithmetic average in terms of values at time \( t_{n+\half} \). According to (666)-(667), $$ \half([D_xD_x \uex]^n_i + [D_xD_x \uex]^{n+1}_i) = [D_xD_x\uex]^{n+\half}_i + \frac{1}{8}[D_xD_x\uexd{tt}]_i^{n+\half}\Delta t^2 + \Oof{\Delta t^4}\tp $$ With (662)-(663) we can express the difference operator \( D_xD_xu \) in terms of a derivative: $$ [D_xD_x\uex]^{n+\half}_i = \uexd{xx}(x_i, t_{n+\half}) + \frac{1}{12}\uexd{xxxx}(x_i, t_{n+\half})\Delta x^2 + \Oof{\Delta x^4}\tp $$ The error term from the arithmetic mean is similarly expanded, $$ \frac{1}{8}[D_xD_x\uexd{tt}]_i^{n+\half}\Delta t^2 = \frac{1}{8}\uexd{ttxx}(x_i, t_{n+\half})\Delta t^2 + \Oof{\Delta t^2\Delta x^2} $$
The time derivative is analyzed using (650)-(651): $$ [D_t u]^{n+\half}_i = \uexd{t}(x_i,t_{n+\half}) + \frac{1}{24}\uexd{ttt}(x_i,t_{n+\half})\Delta t^2 + \Oof{\Delta t^4}\tp $$
Summing up all the contributions and notifying that $$ \uexd{t}(x_i,t_{n+\half}) = \dfc\uexd{xx}(x_i, t_{n+\half}) + f(x_i,t_{n+\half}),$$ the truncation error is given by $$ \begin{align*} R^{n+\half}_i & = \frac{1}{8}\uexd{xx}(x_i,t_{n+\half})\Delta t^2 + \frac{1}{12}\uexd{xxxx}(x_i, t_{n+\half})\Delta x^2 +\\ &\quad \frac{1}{24}\uexd{ttt}(x_i,t_{n+\half})\Delta t^2 + + \Oof{\Delta x^4} + \Oof{\Delta t^4} + \Oof{\Delta t^2\Delta x^2} \end{align*} $$
We address the PDE $$ \frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\left( \dfc(u)\frac{\partial u}{\partial x}\right) + f(u),$$ with two potentially nonlinear coefficients \( q(u) \) and \( \dfc(u) \). We use a Backward Euler scheme with arithmetic mean for \( \dfc(u) \), $$ [D^-u = D_x\overline{\dfc(u)}^{x}D_x u + f(u)]_i^n\tp$$ Inserting \( \uex \) defines the truncation error \( R \): $$ [D^-\uex = D_x\overline{\dfc(\uex)}^{x}D_x \uex + f(\uex)]_i^n\tp$$
The most computationally challenging part is the variable coefficient with \( \dfc(u) \), but we can use the same setup as in the section Extension to variable coefficients and arrive at a truncation error \( \Oof{\Delta x^2} \) for the \( x \)-derivative term. The nonlinear term \( [f(\uex)]=^n_{i} = f(\uex(x_i, t_n)) \) matches \( x \) and \( t \) derivatives of \( \uex \) in the PDE. We end up with $$ R^n_i = -{\half}\frac{\partial^2}{\partial t^2}\uex(x_i,t_n)\Delta t + \Oof{\Delta x^2}\tp$$
Derive the truncation error of the weighted mean in (664)-(665).
Hint. Expand \( \uex^{n+1} \) and \( \uex^n \) around \( t_{n+\theta} \).
Filename: trunc_weighted_mean.
We consider the weighted mean
$$ \uex(t_n) \approx \theta \uex^{n+1} + (1-\theta)\uex^n\tp $$
Choose some specific function for \( \uex(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
(664)-(665).
Filename: trunc_theta_avg.
Set up a numerical experiment as explained in
the section Empirical verification of the truncation error for verifying the formulas
(660)-(661).
Filename: trunc_backward_2level.
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.
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+\half} \), \( 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.
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.
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 \( \Oof{\Delta t^4} \).
Filename: trunc_decay_corr_verify.
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+\half} \):
$$
\begin{align}
\lbrack D_t u &= -a\overline{u}^t + b \rbrack^{n+\half},
\tag{731}\\
\lbrack D_t u &= \overline{-au+b}^t \rbrack^{n+\half}
\tp
\tag{732}
\end{align}
$$
Compute the truncation error in both cases.
Filename: trunc_decay_CN_vc.
Consider the general nonlinear first-order scalar ODE $$ u'(t) = f(u(t), t) \tp $$ 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+\half}, $$
is somewhat more involved: Taylor expand \( \uex^n \), \( \uex^{n+1} \),
\( f(\uex^n, t_n) \), and \( f(\uex^{n+1}, t_{n+1}) \) around \( t_{n+\half} \),
and use that
$$ \frac{df}{dt} = \frac{\partial f}{\partial u}u' + \frac{\partial f}{\partial t}
\tp $$
Check that the derived truncation error is consistent with previous
results for the case \( f(u,t)=-au \).
Filename: trunc_nonlinear_ODE.
Derive the truncation error of the finite difference approximation
(662)-(663) to
the second-order derivative.
Filename: trunc_d2u.
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.
Consider the ODE
$$ mu'' + \beta |u'|u' + s(u) = F(t)\tp$$
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\tp$$
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.