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)
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
sympyto 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.
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}
\tag{1}
\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
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$$
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} \)!