$$ \newcommand{\half}{\frac{1}{2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} $$

Study guide: Truncation error analysis

Hans Petter Langtangen [1, 2]

[1] Center for Biomedical Computing, Simula Research Laboratory [2] Department of Informatics, University of Oslo

Oct 17, 2015

Overview of what truncation errors are
Abstract problem setting
Various error measures
Truncation errors in finite difference formulas
Example: The backward difference for $ u'(t) $
Taylor series
Taylor series inserted in the backward difference approximation
The forward difference for $ u'(t) $
The central difference for $ u'(t) $ (1)
The central difference for $ u'(t) $ (1)
Leading-order error terms in finite differences (1)
Leading-order error terms in finite differences (2)
Leading-order error terms in mean values (1)
Leading-order error terms in mean values (2)
Software for computing truncation errors
Symbolic computing with difference operators
Truncation errors in exponential decay ODE
Truncation error of the Forward Euler scheme
Truncation error of the Crank-Nicolson scheme
Test the understanding!
Truncation error of the $ \theta $-rule
Using symbolic software
Empirical verification of the truncation error (1)
Empirical verification of the truncation error (2)
Empirical verification of the truncation error in the Forward Euler scheme
Empirical verification of the truncation error in the Forward Euler scheme
Increasing the accuracy by adding correction terms
Lowering the order of the derivative in the correction term
With a correction term Forward Euler becomes Crank-Nicolson
Correction terms in the Crank-Nicolson scheme (1)
Correction terms in the Crank-Nicolson scheme (2)
Extension to variable coefficients
Exact solutions of the finite difference equations
Computing truncation errors in nonlinear problems (1)
Computing truncation errors in nonlinear problems (2)
Truncation errors in vibration ODEs
Linear model without damping
Truncation errors in the initial condition
Computing correction terms
Model with damping and nonlinearity
Carrying out the truncation error analysis
Extension to quadratic damping
The truncation error for quadratic damping (1)
The truncation error for quadratic damping (2)
The general model formulated as first-order ODEs
The forward-backward scheme
Truncation error analysis
A centered scheme on a staggered mesh
Truncation error analysis (1)
Truncation error analysis (2)

Overview of what truncation errors are

Definition: The truncation error is the discrepancy that arises from performing a finite number of steps to approximate a process with infinitely many steps.

Widely used: truncation of infinite series, finite precision arithmetic, finite differences, and differential equations.

Why? The truncation error is an error measure that is easy to compute.

Abstract problem setting

Consider an abstract differential equation $$ \mathcal{L}(u)=0$$ Example: $ \mathcal{L}(u)=u'(t)+a(t)u(t)-b(t) $.

The corresponding discrete equation: $$ \mathcal{L}_{\Delta}(u) =0$$ Let now

$ u $ be the numerical solution of the discrete equations, computed at mesh points: $ u^n $, $ n=0,\ldots,N_t $

$ \uex $ the exact solution of the differential equation

$$ \begin{align*} \mathcal{L}(\uex)&=0\\ \mathcal{L}_\Delta(u)&=0 \end{align*} $$ $ u $ is computed at mesh points

Various error measures

Dream: the true error $ e = \uex-u $, but usually impossible

Must find other error measures that are easier to calculate

Derive formulas for $ u $ in (very) special, simplified cases

Compute empirical convergence rates for special choices of $ \uex $ (usually non-physical $ \uex $)

To what extent does $ \uex $ fulfill $ \mathcal{L}_\Delta(\uex)=0 $?

It does not fit, but we can measure the error $ \mathcal{L}_\Delta(\uex)=R $

$ R $ is the truncation error and it is easy to compute in general, without considering special cases

Truncation errors in finite difference formulas

Example: The backward difference for $ u'(t) $

Backward difference approximation to $ u' $: $$ \lbrack D_t^- u\rbrack^n = \frac{u^{n} - u^{n-1}}{\Delta t} \approx u'(t_n) $$

Define the truncation error of this approximation as $$ R^n = [D^-_tu]^n - u'(t_n) $$

The common way of calculating $ R^n $ is to

expand $ u(t) $ in a Taylor series around the point where the derivative is evaluated, here $ t_n $,

insert this Taylor series in \eqref{trunc:fd:bw3}, and

collect terms that cancel and simplify the expression.

Taylor series

General Taylor series expansion from calculus: $$ f(x+h) = \sum_{i=0}^\infty \frac{1}{i!}\frac{d^if}{dx^i}(x)h^i$$

Here: expand $ u^{n-1} $ around $ t_n $: $$ \begin{align*} u(t_{n-1}) = u(t-\Delta t) &= \sum_{i=0}^\infty \frac{1}{i!}\frac{d^iu}{dt^i}(t_n)(-\Delta t)^i\\ & = u(t_n) - u'(t_n)\Delta t + {\half}u''(t_n)\Delta t^2 + \Oof{\Delta t^3} \end{align*} $$

$ \Oof{\Delta t^3} $: power-series in $ \Delta t $ where the lowest power is $ \Delta t^3 $

Small $ \Delta t $: $ \Delta t \gg \Delta t^3 \gg \Delta t^4 $

Taylor series inserted in the backward difference approximation

$$ \begin{align*} [D_t^-u]^n - u'(t_n) &= \frac{u(t_n) - u(t_{n-1})}{\Delta t} - u'(t_n)\\ &= \frac{u(t_n) - (u(t_n) - u'(t_n)\Delta t + {\half}u''(t_n)\Delta t^2 + \Oof{\Delta t^3} )}{\Delta t}\\ &\quad -u'(t_n)\\ &= -{\half}u''(t_n)\Delta t + \Oof{\Delta t^2} ) \end{align*} $$

Result: $$ R^n = -{\half}u''(t_n)\Delta t + \Oof{\Delta t^2} ) $$ The difference approximation is of first order in $ \Delta t $. It is exact for linear $ \uex $.

The forward difference for $ u'(t) $

Now consider a forward difference: $$ u'(t_n) \approx [D_t^+ u]^n = \frac{u^{n+1}-u^n}{\Delta t}$$ Define the truncation error: $$ R^n = [D_t^+ u]^n - u'(t_n)$$ Expand $ u^{n+1} $ in a Taylor series around $ t_n $, $$ u(t_{n+1}) = u(t_n) + u'(t_n)\Delta t + {\half}u''(t_n)\Delta t^2 + \Oof{\Delta t^3} $$ We get $$ R = {\half}u''(t_n)\Delta t + \Oof{\Delta t^2}$$

The central difference for $ u'(t) $ (1)

For the central difference approximation, $$ u'(t_n)\approx [ D_tu]^n, \quad [D_tu]^n = \frac{u^{n+\half} - u^{n-\half}}{\Delta t} $$ the truncation error is $$ R^n = [ D_tu]^n - u'(t_n)$$ Expand $ u(t_{n+\half}) $ and $ u(t_{n-1/2}) $ in Taylor series around the point $ t_n $ where the derivative is evaluated: $$ \begin{align*} u(t_{n+\half}) = &u(t_n) + u'(t_n)\half\Delta t + {\half}u''(t_n)(\half\Delta t)^2 + \\ & \frac{1}{6}u'''(t_n) (\half\Delta t)^3 + \frac{1}{24}u''''(t_n) (\half\Delta t)^4 + \Oof{\Delta t^5}\\ u(t_{n-1/2}) = &u(t_n) - u'(t_n)\half\Delta t + {\half}u''(t_n)(\half\Delta t)^2 - \\ & \frac{1}{6}u'''(t_n) (\half\Delta t)^3 + \frac{1}{24}u''''(t_n) (\half\Delta t)^4 + \Oof{\Delta t^5} \end{align*} $$

The central difference for $ u'(t) $ (1)

$$ u(t_{n+\half}) - u(t_{n-1/2}) = u'(t_n)\Delta t + \frac{1}{24}u'''(t_n) \Delta t^3 + \Oof{\Delta t^5} \tp $$ By collecting terms in $ [D_t u]^n - u(t_n) $ we find $ R^n $ to be $$ R^n = \frac{1}{24}u'''(t_n)\Delta t^2 + \Oof{\Delta t^4} $$ Note:

Second-order accuracy since the leading term is $ \Delta t^2 $

Only even powers of $ \Delta t $

Leading-order error terms in finite differences (1)

$$ \begin{align*} \lbrack D_tu \rbrack^n &= \frac{u^{n+\half} - u^{n-\half}}{\Delta t} = u'(t_n) + R^n\\ R^n &= \frac{1}{24}u'''(t_n)\Delta t^2 + \Oof{\Delta t^4}\\ \lbrack D_{2t}u \rbrack^n &= \frac{u^{n+1} - u^{n-1}}{2\Delta t} = u'(t_n) + R^n\\ R^n &= \frac{1}{6}u'''(t_n)\Delta t^2 + \Oof{\Delta t^4}\\ \lbrack D_t^-u \rbrack^n &= \frac{u^{n} - u^{n-1}}{\Delta t} = u'(t_n) + R^n\\ R^n &= -{\half}u''(t_n)\Delta t + \Oof{\Delta t^2}\\ \lbrack D_t^+u \rbrack^n &= \frac{u^{n+1} - u^{n}}{\Delta t} = u'(t_n) + R^n\\ R^n &= {\half}u''(t_n)\Delta t + \Oof{\Delta t^2} \end{align*} $$

Leading-order error terms in finite differences (2)

$$ \begin{align*} [\bar D_tu]^{n+\theta} &= \frac{u^{n+1} - u^{n}}{\Delta t} = u'(t_{n+\theta}) + R^{n+\theta}\\ R^{n+\theta} &= \half(1-2\theta)u''(t_{n+\theta})\Delta t - \frac{1}{6}((1 - \theta)^3 - \theta^3)u'''(t_{n+\theta})\Delta t^2 + \Oof{\Delta t^3}\\ \lbrack D_t^{2-}u \rbrack^n &= \frac{3u^{n} - 4u^{n-1} + u^{n-2}}{2\Delta t} = u'(t_n) + R^n\\ R^n &= -\frac{1}{3}u'''(t_n)\Delta t^2 + \Oof{\Delta t^3}\\ \lbrack D_tD_t u \rbrack^n &= \frac{u^{n+1} - 2u^{n} + u^{n-1}}{\Delta t^2} = u''(t_n) + R^n\\ R^n &= \frac{1}{12}u''''(t_n)\Delta t^2 + \Oof{\Delta t^4}\\ \end{align*} $$

Leading-order error terms in mean values (1)

Weighted arithmetic mean: $$ \begin{align*} [\overline{u}^{t,\theta}]^{n+\theta} & = \theta u^{n+1} + (1-\theta)u^n = u(t_{n+\theta}) + R^{n+\theta}\\ R^{n+\theta} &= {\half}u''(t_{n+\theta})\Delta t^2\theta (1-\theta) + \Oof{\Delta t^3} \end{align*} $$ Standard arithmetic mean: $$ \begin{align*} [\overline{u}^{t}]^{n} &= \half(u^{n-\half} + u^{n+\half}) = u(t_n) + R^{n}\\ R^{n} &= \frac{1}{8}u''(t_{n})\Delta t^2 + \frac{1}{384}u''''(t_n)\Delta t^4 + \Oof{\Delta t^6} \end{align*} $$

Leading-order error terms in mean values (2)

Geometric mean: $$ \begin{align*} [\overline{u^2}^{t,g}]^{n} &= u^{n-\half}u^{n+\half} = (u^n)^2 + R^n\\ R^n &= - \frac{1}{4}u'(t_n)^2\Delta t^2 + \frac{1}{4}u(t_n)u''(t_n)\Delta t^2 + \Oof{\Delta t^4} \end{align*} $$

Harmonic mean: $$ \begin{align*} [\overline{u}^{t,h}]^{n} &= u^n = \frac{2}{\frac{1}{u^{n-\half}} + \frac{1}{u^{n+\half}}} + R^{n+\half}\\ R^n &= - \frac{u'(t_n)^2}{4u(t_n)}\Delta t^2 + \frac{1}{8}u''(t_n)\Delta t^2 \end{align*} $$

Software for computing truncation errors

Can use sympy to automate calculations with Taylor series.

Tool: course module truncation_errors

>>> from truncation_errors import TaylorSeries
>>> from sympy import *
>>> u, dt = symbols('u dt')
>>> u_Taylor = TaylorSeries(u, 4)
>>> u_Taylor(dt)
D1u*dt + D2u*dt**2/2 + D3u*dt**3/6 + D4u*dt**4/24 + u
>>> FE = (u_Taylor(dt) - u)/dt
>>> FE
(D1u*dt + D2u*dt**2/2 + D3u*dt**3/6 + D4u*dt**4/24)/dt
>>> simplify(FE)
D1u + D2u*dt/2 + D3u*dt**2/6 + D4u*dt**3/24

Notation: D1u for $ u' $, D2u for $ u'' $, etc.

See trunc/truncation_errors.py.

Symbolic computing with difference operators

A class DiffOp represents many common difference operators:

>>> from truncation_errors import DiffOp
>>> from sympy import *
>>> u = Symbol('u')
>>> diffop = DiffOp(u, independent_variable='t')
>>> diffop['geometric_mean']
-D1u**2*dt**2/4 - D1u*D3u*dt**4/48 + D2u**2*dt**4/64 + ...
>>> diffop['Dtm']
D1u + D2u*dt/2 + D3u*dt**2/6 + D4u*dt**3/24
>>> diffop.operator_names()
['geometric_mean', 'harmonic_mean', 'Dtm', 'D2t', 'DtDt',
 'weighted_arithmetic_mean', 'Dtp', 'Dt']

Names in diffop: Dtp for $ D^+_t $, Dtm for $ D_t^- $, Dt for $ D_t $, D2t for $ D_{2t} $, DtDt for $ D_tD_t $.

Truncation errors in exponential decay ODE

$$ u'(t)=-au(t) $$

Truncation error of the Forward Euler scheme

The Forward Euler scheme: $$ \lbrack D_t^+ u = -au \rbrack^n \tp $$ Definition of the truncation error $ R^n $: $$ \lbrack D_t^+ \uex + a\uex = R \rbrack^n \tp $$ From \eqref{trunc:table:fd1:fw:eq}-\eqref{trunc:table:fd1:fw}: $$ [D_t^+ \uex]^n = \uex'(t_n) + \half\uex''(t_n)\Delta t + \Oof{\Delta t^2}\tp$$ Inserted in \eqref{trunc:decay:FE:uex}: $$ \uex'(t_n) + \half\uex''(t_n)\Delta t + \Oof{\Delta t^2} + a\uex(t_n) = R^n \tp $$ Note: $ \uex'(t_n) + a\uex^n = 0 $ since $ \uex $ solves the ODE. Then $$ R^n = \half\uex''(t_n)\Delta t + \Oof{\Delta t^2} \tp $$

Truncation error of the Crank-Nicolson scheme

Crank-Nicolson: $$ [D_t u = -au]^{n+\half} $$ Truncation error: $$ [D_t \uex + a\overline{\uex}^{t} = R]^{n+\half} $$

From \eqref{trunc:table:fd1:center:eq}-\eqref{trunc:table:fd1:center} and \eqref{trunc:table:avg:arith:eq}-\eqref{trunc:table:avg:arith}: $$ \begin{align*} \lbrack D_t\uex \rbrack^{n+\half} &= u'(t_{n+\half}) + \frac{1}{24}\uex'''(t_{n+\half})\Delta t^2 + \Oof{\Delta t^4}\\ [a\overline{\uex}^{t}]^{n+\half} &= u(t_{n+\half}) + \frac{1}{8}u''(t_{n})\Delta t^2 + + \Oof{\Delta t^4} \end{align*} $$ Inserted in the scheme we get $$ R^{n+\half} = \left( \frac{1}{24}\uex'''(t_{n+\half}) + \frac{1}{8}u''(t_{n}) \right)\Delta t^2 + \Oof{\Delta t^4} $$ $ R^n = \Oof{\Delta t^2} $ (second-order scheme)

Test the understanding!

Analyze the the truncation error of the Backward Euler scheme and show that it is $ \Oof{\Delta t} $ (first order scheme).

Truncation error of the $ \theta $-rule

The $ \theta $-rule: $$ [\bar D_t u = -a\overline{u}^{t,\theta}]^{n+\theta} \tp $$

Truncation error: $$ [\bar D_t \uex + a\overline{\uex}^{t,\theta} = R]^{n+\theta} \tp $$ Use \eqref{trunc:table:fd1:theta:eq}-\eqref{trunc:table:fd1:theta} and \eqref{trunc:table:avg:theta:eq}-\eqref{trunc:table:avg:theta} along with $ \uex'(t_{n+\theta}) + a\uex(t_{n+\theta})=0 $ to show

$$ \begin{align} R^{n+\theta} = &({\half}-\theta)\uex''(t_{n+\theta})\Delta t + \half\theta (1-\theta)\uex''(t_{n+\theta})\Delta t^2 + \nonumber\\ & \half(\theta^2 -\theta + 3)\uex'''(t_{n+\theta})\Delta t^2 + \Oof{\Delta t^3} \label{_auto1} \end{align} $$ Note: 2nd-order scheme if and only if $ \theta =1/2 $.

Using symbolic software

Can use sympy and the tools in truncation_errors.py:

def decay():
    u, a = sm.symbols('u a')
    diffop = DiffOp(u, independent_variable='t',
                    num_terms_Taylor_series=3)
    D1u = diffop.D(1)   # symbol for du/dt
    ODE = D1u + a*u     # define ODE

    # Define schemes
    FE = diffop['Dtp'] + a*u
    CN = diffop['Dt' ] + a*u
    BE = diffop['Dtm'] + a*u
    # Residuals (truncation errors)
    R = {'FE': FE-ODE, 'BE': BE-ODE, 'CN': CN-ODE}
    return R

The returned dictionary becomes

decay: {
 'BE': D2u*dt/2 + D3u*dt**2/6,
 'FE': -D2u*dt/2 + D3u*dt**2/6,
 'CN': D3u*dt**2/24,
}

$ \theta $-rule: see truncation_errors.py (long expression, very advantageous to automate the math!)

Empirical verification of the truncation error (1)

Ideas:

Compute $ R^n $ numerically

Run a sequence of meshes

Estimate the convergence rate of $ R^n $

For the Forward Euler scheme: $$ R^n = \lbrack D_t^+\uex + a\uex \rbrack^n \tp $$ Insert correct $ \uex(t)=Ie^{-at} $ (or use method of manufactured solution in more general cases).

Empirical verification of the truncation error (2)

Assume $ R^n = C\Delta t^r $

$ C $ and $ r $ will vary with $ n $ - must estimate $ r $ for each mesh point

Use a sequence of meshes with $ N_t = 2^{-k}N_0 $ intervals, $ k=1,2,\ldots $

Transform $ R^n $ data to the coarsest mesh and estimate $ r $ for each coarse mesh point

See the text for more details and an implementation.

Empirical verification of the truncation error in the Forward Euler scheme

Figure 1: Estimated truncation error at mesh points for different meshes.

Empirical verification of the truncation error in the Forward Euler scheme

Figure 2: Difference between theoretical and estimated truncation error at mesh points for different meshes.

Increasing the accuracy by adding correction terms

Question.

Can we add terms in the differential equation that can help increase the order of the truncation error?

To be precise for the Forward Euler scheme, can we find $ C $ to make $ R $ $\Oof{\Delta t^2}$? $$ \lbrack D_t^+ \uex + a\uex = C + R \rbrack^n\tp $$

$$ \half\uex''(t_n)\Delta t - \frac{1}{6}\uex'''(t_n)\Delta t^2 + \Oof{\Delta t^3} = C^n + R^n\tp$$ Choosing $$ C^n = \half\uex''(t_n)\Delta t$$ makes $$ R^n = \frac{1}{6}\uex'''(t_n)\Delta t^2 + \Oof{\Delta t^3}\tp$$

Lowering the order of the derivative in the correction term

$ C^n $ contains $ u'' $

Can discretize $ u'' $ (requires $ u^{n+1} $, $ u^n $, and $ u^{n-1} $)

Can also express $ u'' $ in terms of $ u' $ or $ u $

$$ u' = -au,\quad\Rightarrow\quad u''=-au'=a^2u\tp$$ Result for $ u''=a^2u $: apply Forward Euler to a perturbed ODE, $$ u' = -\hat au ,\quad \hat a = a(1 - {\half}a\Delta t) $$ to make a second-order scheme!

With a correction term Forward Euler becomes Crank-Nicolson

Use the other alternative $ u''=-au' $: $$ u'=-au - {\half}a\Delta t u'\quad\Rightarrow\quad \left( 1 + {\half}a\Delta t\right) u' = -au\tp$$ Apply Forward Euler: $$ \left( 1 + {\half}a\Delta t\right)\frac{u^{n+1}-u^n}{\Delta t} = -au^n$$ which after some algebra can be written as $$ u^{n+1} = \frac{1 - {\half}a\Delta t}{1+{\half}a\Delta t}u^n\tp$$ This is a Crank-Nicolson scheme (of second order)!

Correction terms in the Crank-Nicolson scheme (1)

$$ [D_t u = -a\overline{u}^t]^{n+\half}$$ Definition of the truncation error $ R $ and correction terms $ C $: $$ [D_t \uex + a\overline{\uex}^{t} = C + R]^{n+\half}\tp$$

Must Taylor expand

the derivative

the arithmetic mean

$$ C^{n+\half} + R^{n+\half} = \frac{1}{24}\uex'''(t_{n+\half})\Delta t^2 + \frac{a}{8}\uex''(t_{n+\half})\Delta t^2 + \Oof{\Delta t^4}\tp$$ Let $ C^{n+\half} $ cancel the $ \Delta t^2 $ terms: $$ C^{n+\half} = \frac{1}{24}\uex'''(t_{n+\half})\Delta t^2 + \frac{a}{8}\uex''(t_{n})\Delta t^2\tp$$

Correction terms in the Crank-Nicolson scheme (2)

Must replace $ u''' $ and $ u'' $ in correction term

Using $ u'=-au $: $ u''=a^2u $ and $ u'''=-a^3u $

Result: solve the perturbed ODE by a Crank-Nicolson method, $$ u' = -\hat a u,\quad \hat a = a(1 - \frac{1}{12}a^2\Delta t^2)\tp$$ and experience an error $ \Oof{\Delta t^4} $.

Extension to variable coefficients

$$ u'(t) = -a(t)u(t) + b(t)$$

Forward Euler: $$ [D_t^+ u = -au + b]^n \tp $$ The truncation error is found from $$ [D_t^+ \uex + a\uex - b = R]^n \tp $$ Using \eqref{trunc:table:fd1:fw:eq}-\eqref{trunc:table:fd1:fw}: $$ \uex'(t_n) - \half\uex''(t_n)\Delta t + \Oof{\Delta t^2} + a(t_n)\uex(t_n) - b(t_n) = R^n \tp $$ Because of the ODE, $ \uex'(t_n) + a(t_n)\uex(t_n) - b(t_n) =0 $, and $$ R^n = -\half\uex''(t_n)\Delta t + \Oof{\Delta t^2} \tp $$ No problems with variable coefficients!

Exact solutions of the finite difference equations

How does the truncation error depend on $ \uex $ in finite differences?

One-sided differences: $ \uex''\Delta t $ (lowest order)

Centered differences: $ \uex'''\Delta t^2 $ (lowest order)

Only harmonic and geometric mean involve $ \uex' $ or $ \uex $

Consequence:

$ \uex(t)=ct+d $ will very often give exact solution of the discrete equations ($ R=0 $)!

Ideal for verification

Centered schemes allow quadratic $ \uex $

Problem: harmonic and geometric mean (error depends on $ \uex' $ and $ \uex $)

Computing truncation errors in nonlinear problems (1)

$$ u'=f(u,t) $$ Crank-Nicolson scheme: $$ [D_t u =\overline{f}^{t}]^{n+\half} $$ Truncation error: $$ [D_t \uex - \overline{f}^{t}= R]^{n+\half}\tp $$ Using \eqref{trunc:table:avg:arith:eq}-\eqref{trunc:table:avg:arith} for the arithmetic mean: $$ \begin{align*} \lbrack\overline{f}^{t}\rbrack^{n+\half} &= \half(f(\uex^n,t_n) + f(\uex^{n+1},t_{n+1}))\\ &= f(\uex^{n+\half},t_{n+\half}) + \frac{1}{8}\uex''(t_{n+\half})\Delta t^2 + \Oof{\Delta t^4}\tp \end{align*} $$

Computing truncation errors in nonlinear problems (2)

With \eqref{trunc:table:fd1:center:eq}-\eqref{trunc:table:fd1:center}, \eqref{trunc:decay:gen:ode:CN} leads to $ R^{n+\half} $ equal to $$ \uex'(t_{n+\half}) + \frac{1}{24}\uex'''(t_{n+\half})\Delta t^2 - f(\uex^{n+\half},t_{n+\half}) - \frac{1}{8}\uex''(t_{n+\half})\Delta t^2 + \Oof{\Delta t^4} \tp $$ Since $ \uex'(t_{n+\half}) - f(\uex^{n+\half},t_{n+\half})=0 $, the truncation error becomes $$ R^{n+\half} = (\frac{1}{24}\uex'''(t_{n+\half}) - \frac{1}{8}\uex''(t_{n+\half})) \Delta t^2\tp $$ The computational techniques worked well even for this nonlinear ODE!

Truncation errors in vibration ODEs

Linear model without damping

$$ u''(t) + \omega^2 u(t) = 0,\quad u(0)=I,\ u'(0)=0\tp $$ Centered difference approximation: $$ [D_tD_t u + \omega^2u=0]^n $$ Truncation error: $$ [D_tD_t \uex + \omega^2\uex =R]^n \tp $$ Use \eqref{trunc:table:fd2:center:eq}-\eqref{trunc:table:fd2:center} to expand $ [D_t D_t\uex]^n $: $$ [D_tD_t \uex]^n = \uex''(t_n) + \frac{1}{12}\uex''''(t_n)\Delta t^2$$ Collect terms: $ \uex''(t) + \omega^2\uex(t)=0 $. Then, $$ R^n = \frac{1}{12}\uex''''(t_n)\Delta t^2 + \Oof{\Delta t^4} \tp $$

Truncation errors in the initial condition

Initial conditions: $ u(0)=I $, $ u'(0)=V $

Need discretization of $ u'(0) $

Standard, centered difference: $ [D_{2t} u=V]^0 $, $ R^0=\Oof{\Delta t^2} $

Simpler, forward difference: $ [D_t^+u=V]^0 $, $ R^0=\Oof{\Delta t} $

Does the lower order of the forward scheme impact the order of the whole simulation?

Answer: run experiments!

Computing correction terms

Can we add terms to the ODE such that the truncation error is improved?

$$ [D_tD_t \uex + \omega^2\uex =C + R]^n$$

Idea: choose $ C^n $ such that it absorbs the $ \Delta t^2 $ term in $ R^n $,

$$ C^n = \frac{1}{12}\uex''''(t_n)\Delta t^2\tp$$

Downside: got a $ u'''' $ term

Remedy: use the ODE $ u''=-\omega^2u $ to see that $ u'''' = \omega^4 u $.

Just apply the standard scheme to a modified ODE:

$$ [D_tD_t u + \omega^2(1 - \frac{1}{12}\omega^2\Delta t^2)u=0]^n$$

Accuracy is $ \Oof{\Delta t^4} $.

Model with damping and nonlinearity

Linear damping $ \beta u' $, nonlinear spring force $ s(u) $, and excitation $ F $: $$ mu'' + \beta u' + s(u) =F(t) $$ Central difference discretization: $$ [mD_tD_t u + \beta D_{2t} u + s(u)=F]^n \tp $$ Truncation error is defined by $$ [mD_tD_t \uex + \beta D_{2t} \uex + s(\uex)=F + R]^n \tp $$

Carrying out the truncation error analysis

Using \eqref{trunc:table:fd2:center:eq}-\eqref{trunc:table:fd2:center} and \eqref{trunc:table:fd1:center2:eq}-\eqref{trunc:table:fd1:center2} we get $$ \begin{align*} \lbrack mD_tD_t \uex + \beta D_{2t} \uex\rbrack^n &= m\uex''(t_n) + \beta\uex'(t_n) + \\ &\quad \left(\frac{m}{12}\uex''''(t_n) + \frac{\beta}{6}\uex'''(t_n)\right)\Delta t^2 + \Oof{\Delta t^4} \end{align*} $$ The terms $$ m\uex''(t_n) + \beta\uex'(t_n) + \omega^2\uex(t_n) + s(\uex(t_n)) - F^n$$ correspond to the ODE (= zero).

Result: accuracy of $ \Oof{\Delta t^2} $ since $$ R^n = \left(\frac{m}{12}\uex''''(t_n) + \frac{\beta}{6}\uex'''(t_n)\right)\Delta t^2 + \Oof{\Delta t^4} $$

Correction terms: complicated when the ODE has many terms...

Extension to quadratic damping

$$ mu'' + \beta |u'|u' + s(u) =F(t) $$ Centered scheme: $ |u'|u' $ gives rise to a nonlinearity.

Linearization trick: use a geometric mean, $$ [|u'|u']^n \approx |[u']^{n-\half}|[u']^{n+\half}\tp$$

Scheme: $$ [mD_t D_t u]^n + \beta |[D_{t} u]^{n-\half}|[D_t u]^{n+\half} + s(u^n)=F^n\tp $$

The truncation error for quadratic damping (1)

Definition of $ R^n $: $$ [mD_t D_t \uex]^n + \beta |[D_{t} \uex]^{n-\half}|[D_t \uex]^{n+\half} + s(\uex^n)-F^n = R^n\tp $$

Truncation error of the geometric mean, see \eqref{trunc:table:avg:geom:eq}-\eqref{trunc:table:avg:geom}, $$ \begin{align*} |[D_{t} \uex]^{n-\half}|[D_t \uex]^{n+\half} &= [|D_t\uex|D_t\uex]^n - \frac{1}{4}u'(t_n)^2\Delta t^2 + \\ &\quad \frac{1}{4}u(t_n)u''(t_n)\Delta t^2 + \Oof{\Delta t^4}\tp \end{align*} $$ Using \eqref{trunc:table:fd1:center:eq}-\eqref{trunc:table:fd1:center} for the $ D_t\uex $ factors results in $$ \begin{align*} [|D_t\uex|D_t\uex]^n &= |\uex' + \frac{1}{24}\uex'''(t_n)\Delta t^2 + \Oof{\Delta t^4}|\times\\ &\quad (\uex' + \frac{1}{24}\uex'''(t_n)\Delta t^2 + \Oof{\Delta t^4}) \end{align*} $$

The truncation error for quadratic damping (2)

For simplicity, remove the absolute value. The product becomes $$ [D_t\uex D_t\uex]^n = (\uex'(t_n))^2 + \frac{1}{12}\uex(t_n)\uex'''(t_n)\Delta t^2 + \Oof{\Delta t^4}\tp$$

With $$ m[D_t D_t\uex]^n = m\uex''(t_n) + \frac{m}{12}\uex''''(t_n)\Delta t^2 +\Oof{\Delta t^4}$$ and using $ mu'' + \beta (u')^2 + s(u)=F $, we end up with $$ R^n = (\frac{m}{12}\uex''''(t_n) + \frac{\beta}{12}\uex(t_n)\uex'''(t_n)) \Delta t^2 + \Oof{\Delta t^4}\tp$$ Second-order accuracy! Thanks to

difference approximation with error $ \Oof{\Delta t^2} $

geometric mean approximation with error $ \Oof{\Delta t^2} $

The general model formulated as first-order ODEs

$$ mu'' + \beta |u'|u' + s(u) =F(t) $$ Rewritten as first-order system: $$ \begin{align*} u' &= v\\ v' &= \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right) \end{align*} $$

To solution methods:

Forward-backward scheme

Centered scheme on a staggered mesh

The forward-backward scheme

Forward step for $ u $, backward step for $ v $: $$ \begin{align*} [D_t^+ u &= v]^n\\ [D_t^-v &= \frac{1}{m}( F(t) - \beta |v|v - s(u))]^{n+1} \end{align*} $$

Note:

step $ u $ forward with known $ v $ in \eqref{trunc:vib:gen:2x2model:ode:u:fw}

step $ v $ forward with known $ u $ in \eqref{trunc:vib:gen:2x2model:ode:v:bw}

Problem: $ |v|v $ gives nonlinearity $ |v^{n+1}|v^{n+1} $.

Remedy: linearized as $ |v^{n}|v^{n+1} $

$$ \begin{align*} [D_t^+ u &= v]^n\\ [D_t^-v]^{n+1} &= \frac{1}{m}( F(t_{n+1}) - \beta |v^n|v^{n+1} - s(u^{n+1})) \end{align*} $$

Truncation error analysis

Aim (as always): turn difference operators into derivatives + truncation error terms

One-sided forward/backward differences: error $ \Oof{\Delta t} $

Linearization of $ |v^{n+1}|v^{n+1} $ to $ |v^n|v^{n+1} $: error $ \Oof{\Delta t} $

All errors are $ \Oof{\Delta t} $

First-order scheme? No!

"Symmetric" use of the $ \Oof{\Delta t} $ building blocks yields in fact a $ \Oof{\Delta t^2} $ scheme (!)

Why? See next slide...

A centered scheme on a staggered mesh

Staggered mesh:

$ u $ is computed at mesh points $ t_n $

$ v $ is computed at points $ t_{n+\half} $

Centered differences in \eqref{trunc:vib:gen:2x2model:ode:u}-\eqref{trunc:vib:gen:2x2model:ode:u}: $$ \begin{align*} [D_t u &= v]^{n-\half}\\ [D_t v &= \frac{1}{m}( F(t) - \beta |v|v - s(u))]^{n} \end{align*} $$

Problem: $ |v^n|v^n $ because $ v^n $ is not computed directly

Remedy: Geometric mean,

$$ |v^n|v^n \approx |v^{n-\half}|v^{n+\half}\tp$$

Truncation error analysis (1)

Resulting scheme: $$ \begin{align*} [D_t u]^{n-\half} &= v^{n-\half}\\ [D_t v]^n &= \frac{1}{m}( F(t_n) - \beta |v^{n-\half}|v^{n+\half} - s(u^n)) \end{align*} $$

The truncation error in each equation is found from $$ \begin{align*} [D_t \uex]^{n-\half} &= \vex(t_{n-\half}) + R_u^{n-\half}\\ [D_t \vex]^n &= \frac{1}{m}( F(t_n) - \beta |\vex(t_{n-\half})|\vex(t_{n+\half}) - s(u^n)) + R_v^n\tp \end{align*} $$ Using \eqref{trunc:table:fd1:center:eq}-\eqref{trunc:table:fd1:center} for derivatives and \eqref{trunc:table:avg:geom:eq}-\eqref{trunc:table:avg:geom} for the geometric mean: $$ \uex'(t_{n-\half}) + \frac{1}{24}\uex'''(t_{n-\half})\Delta t^2 + \Oof{\Delta t^4} = \vex(t_{n-\half}) + R_u^{n-\half}$$ and $$ \vex'(t_n) = \frac{1}{m}( F(t_n) - \beta |\vex(t_n)|\vex(t_n) + \Oof{\Delta t^2} - s(u^n)) + R_v^n\tp $$

Truncation error analysis (2)

Resulting truncation error is $ \Oof{\Delta t^2} $: $$ R_u^{n-\half}= \Oof{\Delta t^2}, \quad R_v^n = \Oof{\Delta t^2}\tp$$

Observation.

Comparing The schemes \eqref{trunc:vib:gen:2x2model:ode:u:staggered2}-\eqref{trunc:vib:gen:2x2model:ode:v:staggered2} and \eqref{trunc:vib:gen:2x2model:ode:u:fw2}-\eqref{trunc:vib:gen:2x2model:ode:v:bw2} are equivalent. Therefore, the forward/backward scheme with ad hoc linearization is also $ \Oof{\Delta t^2} $!

Study guide: Truncation error analysis

Oct 17, 2015

Table of contents

Overview of what truncation errors are

Abstract problem setting

Various error measures

Truncation errors in finite difference formulas

Example: The backward difference for \( u'(t) \)

Taylor series

Taylor series inserted in the backward difference approximation

The forward difference for \( u'(t) \)

The central difference for \( u'(t) \) (1)

The central difference for \( u'(t) \) (1)

Leading-order error terms in finite differences (1)

Leading-order error terms in finite differences (2)

Leading-order error terms in mean values (1)

Leading-order error terms in mean values (2)

Software for computing truncation errors

Symbolic computing with difference operators

Truncation errors in exponential decay ODE

Truncation error of the Forward Euler scheme

Truncation error of the Crank-Nicolson scheme

Test the understanding!

Truncation error of the \( \theta \)-rule

Using symbolic software

Empirical verification of the truncation error (1)

Empirical verification of the truncation error (2)

Empirical verification of the truncation error in the Forward Euler scheme

Empirical verification of the truncation error in the Forward Euler scheme

Increasing the accuracy by adding correction terms

Lowering the order of the derivative in the correction term

With a correction term Forward Euler becomes Crank-Nicolson

Correction terms in the Crank-Nicolson scheme (1)

Correction terms in the Crank-Nicolson scheme (2)

Extension to variable coefficients

Exact solutions of the finite difference equations

Computing truncation errors in nonlinear problems (1)

Computing truncation errors in nonlinear problems (2)

Truncation errors in vibration ODEs

Linear model without damping

Truncation errors in the initial condition

Computing correction terms

Model with damping and nonlinearity

Carrying out the truncation error analysis

Extension to quadratic damping

The truncation error for quadratic damping (1)

The truncation error for quadratic damping (2)

The general model formulated as first-order ODEs

The forward-backward scheme

Truncation error analysis

A centered scheme on a staggered mesh

Truncation error analysis (1)

Truncation error analysis (2)