Truncation errors in vibration ODEs

Linear model without damping

The next example on computing the truncation error involves the following ODE for vibration problems:

(1)\[ u''(t) + \omega^2 u(t) = 0{\thinspace .}\]

Here, \(\omega\) is a given constant.

The truncation error of a centered finite difference scheme

Using a standard, second-ordered, central difference for the second-order derivative time, we have the scheme

(2)\[ [D_tD_t u + \omega^2u=0]^n\]\[ {\thinspace .}\]

Inserting the exact solution \({u_{\small\mbox{e}}}\) in this equation and adding a residual \(R\) so that \({u_{\small\mbox{e}}}\) can fulfill the equation results in

\[[D_tD_t {u_{\small\mbox{e}}} + \omega^2{u_{\small\mbox{e}}} =R]^n {\thinspace .}\]

To calculate the truncation error \(R^n\), we use (1.15)-(1.16), i.e.,

\[[D_tD_t {u_{\small\mbox{e}}}]^n = {u_{\small\mbox{e}}}''(t_n) + \frac{1}{12}{u_{\small\mbox{e}}}''''(t_n)\Delta t^2,\]

and the fact that \({u_{\small\mbox{e}}}''(t) + \omega^2{u_{\small\mbox{e}}}(t)=0\). The result is

\[R^n = \frac{1}{12}{u_{\small\mbox{e}}}''''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)} {\thinspace .}\]

The truncation error of approximating \(u'(0)\)

The initial conditions for (1) are \(u(0)=I\) and \(u'(0)=V\). The latter involves a finite difference approximation. The standard choice

\[[D_{2t}u=V]^0,\]

where \(u^{-1}\) is eliminated with the aid of the discretized ODE for \(n=0\), involves a centered difference with an \({\mathcal{O}(\Delta t^2)}\) truncation error given by (1.5)-(1.6). The simpler choice

\[[D_t^+u = V]^0,\]

is based on a forward difference with a truncation error \({\mathcal{O}(\Delta t)}\). A central question is if this initial error will impact the order of the scheme throughout the simulation. Exercise 11: Investigate the impact of approximating asks you to quickly perform an experiment to investigate this question.

Truncation error of the equation for the first step

We have shown that the truncation error of the difference used to approximate the initial condition \(u'(0)=0\) is \({\mathcal{O}(\Delta t^2)}\), but can also investigate the difference equation used for the first step. In a truncation error setting, the right way to view this equation is not to use the initial condition \([D_{2t}u=V]^0\) to express \(u^{-1}=u^1 - 2\Delta t V\) in order to eliminate \(u^{-1}\) from the discretized differential equation, but the other way around: the fundamental equation is the discretized initial condition \([D_{2t}u=V]^0\) and we use the discretized ODE \([D_tD_t + \omega^2 u=0]^0\) to eliminate \(u^{-1}\) in the disretized initial condition. From \([D_tD_t + \omega^2 u=0]^0\) we have

\[u^{-1} = 2u^0 - u^1 - \Delta t^2\omega^2 u^0,\]

which inserted in \([D_{2t}u = V]^0\) gives

(3)\[ \frac{u^1 - u^0}{\Delta t} + \frac{1}{2}\omega^2\Delta t u^0 = V{\thinspace .}\]

The first term can be recognized as a forward difference such that the equation can be written in operator notation as

\[[D_t^+ u + \frac{1}{2}\omega^2\Delta t u = V]^0{\thinspace .}\]

The truncation error is defined as

\[[D_t^+ {u_{\small\mbox{e}}} + \frac{1}{2}\omega^2\Delta t {u_{\small\mbox{e}}} - V = R]^0{\thinspace .}\]

Using (1.9)-(1.10) with one more term in the Taylor series, we get that

\[{u_{\small\mbox{e}}}'(0) + \frac{1}{2}{u_{\small\mbox{e}}}''(0)\Delta t + \frac{1}{6}{u_{\small\mbox{e}}}'''(0)\Delta t^2 + {\mathcal{O}(\Delta t^3)} + \frac{1}{2}\omega^2\Delta t {u_{\small\mbox{e}}}(0) - V = R^n{\thinspace .}\]

Now, \({u_{\small\mbox{e}}}'(0)=V\) and \({u_{\small\mbox{e}}}''(0)=-\omega^2 {u_{\small\mbox{e}}}(0)\) so we get

\[R^n = \frac{1}{6}{u_{\small\mbox{e}}}'''(0)\Delta t^2 + {\mathcal{O}(\Delta t^3)}{\thinspace .}\]

There is another way of analyzing the discrete initial condition, because eliminating \(u^{-1}\) via the discretized ODE can be expressed as

(4)\[ [ D_{2t} u + \Delta t(D_tD_t u - \omega^2 u) = V]^0{\thinspace .}\]

Writing out (4) shows that the equation is equivalent to (3). The truncation error is defined by

\[[ D_{2t} {u_{\small\mbox{e}}} + \Delta t(D_tD_t {u_{\small\mbox{e}}} - \omega^2 {u_{\small\mbox{e}}}) = V + R]^0{\thinspace .}\]

Replacing the difference via (1.5)-(1.6) and (1.15)-(1.16), as well as using \({u_{\small\mbox{e}}}'(0)=V\) and \({u_{\small\mbox{e}}}''(0) = -\omega^2{u_{\small\mbox{e}}}(0)\), gives

\[R^n = \frac{1}{6}u'''(0)\Delta t^2 + {\mathcal{O}(\Delta t^3)}{\thinspace .}\]

Computing correction terms

The idea of using correction terms to increase the order of \(R^n\) can be applied as described in the section Increasing the accuracy by adding correction terms. We look at

\[[D_tD_t {u_{\small\mbox{e}}} + \omega^2{u_{\small\mbox{e}}} =C + R]^n,\]

and observe that \(C^n\) must be chosen to cancel the \(\Delta t^2\) term in \(R^n\). That is,

\[C^n = \frac{1}{12}{u_{\small\mbox{e}}}''''(t_n)\Delta t^2{\thinspace .}\]

To get rid of the 4th-order derivative we can use the differential equation: \(u''=-\omega^u\), which implies \(u'''' = \omega^4 u\). Adding the correction term to the ODE results in

(5)\[ u'' + \omega^2(1 - \frac{1}{12}\omega^2\Delta t^2)u = 0{\thinspace .}\]

Solving this equation by the standard scheme

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

will result in a scheme with trunction error \({\mathcal{O}(\Delta t^4)}\).

We can use another set of arguments to justify that (5) leads to a higher-order method. Mathematical analysis of the scheme (2) reveals that the numerical frequency \(\tilde\omega\) is (approximately as \(\Delta t\rightarrow 0\))

\[\tilde\omega = \omega (1+\frac{1}{24}\omega^2\Delta t^2){\thinspace .}\]

One can therefore attempt to replace \(\omega\) in the ODE by a slightly smaller \(\omega\) since the numerics will make it larger:

\[[ u'' + (\omega(1 - \frac{1}{24}\omega^2\Delta t^2))^2 u = 0{\thinspace .}\]

Expanding the squared term and omitting the higher-order term \(\Delta t^4\) gives exactly the ODE (5). Experiments show that \(u^n\) is computed to 4th order in \(\Delta t\).

Model with damping and nonlinearity

The model (1) can be extended to include damping \(\beta u'\), a nonlinear restoring (spring) force \(s(u)\), and some known excitation force \(F(t)\):

(6)\[ mu'' + \beta u' + s(u) =F(t){\thinspace .}\]

The coefficient \(m\) usually represents the mass of the system. This governing equation can by discretized by centered differences:

\[[mD_tD_t u + \beta D_{2t} u + s(u)=F]^n {\thinspace .}\]

The exact solution \({u_{\small\mbox{e}}}\) fulfills the discrete equations with a residual term:

\[[mD_tD_t {u_{\small\mbox{e}}} + \beta D_{2t} {u_{\small\mbox{e}}} + s({u_{\small\mbox{e}}})=F + R]^n {\thinspace .}\]

Using (1.15)-(1.16) and (1.5)-(1.6) we get

\[\begin{split}\lbrack mD_tD_t {u_{\small\mbox{e}}} + \beta D_{2t} {u_{\small\mbox{e}}}\rbrack^n &= m{u_{\small\mbox{e}}}''(t_n) + \beta{u_{\small\mbox{e}}}'(t_n) + \\ &\quad \left(\frac{m}{12}{u_{\small\mbox{e}}}''''(t_n) + \frac{\beta}{6}{u_{\small\mbox{e}}}'''(t_n)\right)\Delta t^2 + {\mathcal{O}(\Delta t^4)}\end{split}\]

Combining this with the previous equation, we can collect the terms

\[m{u_{\small\mbox{e}}}''(t_n) + \beta{u_{\small\mbox{e}}}'(t_n) + \omega^2{u_{\small\mbox{e}}}(t_n) + s({u_{\small\mbox{e}}}(t_n)) - F^n,\]

and set this sum to zero because \({u_{\small\mbox{e}}}\) solves the differential equation. We are left with the truncation error

(7)\[ R^n = \left(\frac{m}{12}{u_{\small\mbox{e}}}''''(t_n) + \frac{\beta}{6}{u_{\small\mbox{e}}}'''(t_n)\right)\Delta t^2 + {\mathcal{O}(\Delta t^4)},\]

so the scheme is of second order.

According to (7), we can add correction terms

\[C^n = \left(\frac{m}{12}{u_{\small\mbox{e}}}''''(t_n) + \frac{\beta}{6}{u_{\small\mbox{e}}}'''(t_n)\right)\Delta t^2,\]

to the right-hand side of the ODE to obtain a fourth-order scheme. However, expressing \(u''''\) and \(u'''\) in terms of lower-order derivatives is now harder because the differential equation is more complicated:

\[\begin{split}u''' &= \frac{1}{m}(F' - \beta u'' - s'(u)u'),\\ u'''' &= \frac{1}{m}(F'' - \beta u''' - s''(u)(u')^2 - s'(u)u''),\\ &= \frac{1}{m}(F'' - \beta \frac{1}{m}(F' - \beta u'' - s'(u)u') - s''(u)(u')^2 - s'(u)u''){\thinspace .}\end{split}\]

It is not impossible to discretize the resulting modified ODE, but it is up to debate whether correction terms are feasible and the way to go. Computing with a smaller \(\Delta t\) is usually always possible in these problems to achieve the desired accuracy.

Extension to quadratic damping

Instead of the linear damping term \(\beta u'\) in (6) we now consider quadratic damping \(\beta |u'|u'\):

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

A centered difference for \(u'\) gives rise to a nonlinearity, which can be linearized using a geometric mean: \([|u'|u']^n \approx |[u']^{n-\frac{1}{2}}|[u']^{n+\frac{1}{2}}\). The resulting scheme becomes

\[[mD_t D_t u]^n + \beta |[D_{t} u]^{n-\frac{1}{2}}|[D_t u]^{n+\frac{1}{2}} + s(u^n)=F^n{\thinspace .}\]

The truncation error is defined through

\[[mD_t D_t {u_{\small\mbox{e}}}]^n + \beta |[D_{t} {u_{\small\mbox{e}}}]^{n-\frac{1}{2}}|[D_t {u_{\small\mbox{e}}}]^{n+\frac{1}{2}} + s({u_{\small\mbox{e}}}^n)-F^n = R^n{\thinspace .}\]

We start with expressing the truncation error of the geometric mean. According to (1.21)-(1.22),

\[|[D_{t} {u_{\small\mbox{e}}}]^{n-\frac{1}{2}}|[D_t {u_{\small\mbox{e}}}]^{n+\frac{1}{2}} = [|D_t{u_{\small\mbox{e}}}|D_t{u_{\small\mbox{e}}}]^n - \frac{1}{4}u'(t_n)^2\Delta t^2 + \frac{1}{4}u(t_n)u''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .}\]

Using (1.3)-(1.4) for the \(D_t{u_{\small\mbox{e}}}\) factors results in

\[[|D_t{u_{\small\mbox{e}}}|D_t{u_{\small\mbox{e}}}]^n = |{u_{\small\mbox{e}}}' + \frac{1}{24}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)}|({u_{\small\mbox{e}}}' + \frac{1}{24}{u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)})\]

We can remove the absolute value since it essentially gives a factor 1 or -1 only. Calculating the product, we have the leading-order terms

\[[D_t{u_{\small\mbox{e}}} D_t{u_{\small\mbox{e}}}]^n = ({u_{\small\mbox{e}}}'(t_n))^2 + \frac{1}{12}{u_{\small\mbox{e}}}(t_n){u_{\small\mbox{e}}}'''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .}\]

With

\[m[D_t D_t{u_{\small\mbox{e}}}]^n = m{u_{\small\mbox{e}}}''(t_n) + \frac{m}{12}{u_{\small\mbox{e}}}''''(t_n)\Delta t^2 +{\mathcal{O}(\Delta t^4)},\]

and using the differential equation on the form \(mu'' + \beta (u')^2 + s(u)=F\), we end up with

\[R^n = (\frac{m}{12}{u_{\small\mbox{e}}}''''(t_n) + \frac{\beta}{12}{u_{\small\mbox{e}}}(t_n){u_{\small\mbox{e}}}'''(t_n)) \Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .}\]

This result demonstrates that we have second-order accuracy also with quadratic damping. The key elements that lead to the second-order accuracy is that the difference approximations are \({\mathcal{O}(\Delta t^2)}\) and the geometric mean approximation is also of \({\mathcal{O}(\Delta t^2)}\).

The general model formulated as first-order ODEs

The second-order model (8) can be formulated as a first-order system,

(9)\[ u' = v,\]

(10)\[ v' = \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right){\thinspace .}\]

The system (9)-(9) can be solved either by a forward-backward scheme or a centered scheme on a staggered mesh.

The forward-backward scheme

The discretization is based on the idea of stepping (9) forward in time and then using a backward difference in (10) with the recently computed (and therefore known) \(u\):

(11)\[ \lbrack D_t^+ u = v \rbrack^n,\]

(12)\[ \lbrack D_t^-v = \frac{1}{m}( F(t) - \beta |v|v - s(u)) \rbrack^{n+1}{\thinspace .}\]

The term \(|v|v\) gives rise to a nonlinearity \(|v^{n+1}|v^{n+1}\), which can be linearized as \(|v^{n}|v^{n+1}\):

(13)\[ \lbrack D_t^+ u = v \rbrack^n,\]

(14)\[ \lbrack D_t^-v \rbrack^{n+1} = \frac{1}{m}( F(t_{n+1}) - \beta |v^n|v^{n+1} - s(u^{n+1})){\thinspace .}\]

Each ODE will have a truncation error when inserting the exact solutions \({u_{\small\mbox{e}}}\) and \({v_{\small\mbox{e}}}\) in (11)-(12):

(15)\[ \lbrack D_t^+ {u_{\small\mbox{e}}} = {v_{\small\mbox{e}}} + R_u \rbrack^n,\]

(16)\[ \lbrack D_t^-{v_{\small\mbox{e}}} \rbrack^{n+1} = \frac{1}{m}( F(t_{n+1}) - \beta |{v_{\small\mbox{e}}}(t_n)|{v_{\small\mbox{e}}}(t_{n+1}) - s({u_{\small\mbox{e}}}(t_{n+1}))) + R_v^{n+1}{\thinspace .}\]

Application of (1.9)-(1.10) and (1.7)-(1.8) in (15) and (16), respectively, gives

(17)\[ {u_{\small\mbox{e}}}'(t_n) + \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)} = {v_{\small\mbox{e}}}(t_n) + R_u^n,\]

\[{v_{\small\mbox{e}}}'(t_{n+1}) - \frac{1}{2}{v_{\small\mbox{e}}}''(t_{n+1})\Delta t + {\mathcal{O}(\Delta t^2)} = \frac{1}{m}(F(t_{n+1}) - \beta|{v_{\small\mbox{e}}}(t_n)|{v_{\small\mbox{e}}}(t_{n+1}) +\nonumber\]

(18)\[ \quad s({u_{\small\mbox{e}}}(t_{n+1}))+ R_v^n{\thinspace .}\]

Since \({u_{\small\mbox{e}}} ' = {v_{\small\mbox{e}}}\), (17) gives

\[R_u^n = \frac{1}{2}{u_{\small\mbox{e}}}''(t_n)\Delta t + {\mathcal{O}(\Delta t^2)}{\thinspace .}\]

In (18) we can collect the terms that constitute the ODE, but the damping term has the wrong form. Let us drop the absolute value in the damping term for simplicity. Adding a substracting the right form \(v^{n+1}v^{n+1}\) helps:

\[\begin{split}{v_{\small\mbox{e}}}'(t_{n+1}) &- \frac{1}{m}(F(t_{n+1}) - \beta {v_{\small\mbox{e}}}(t_{n+1}){v_{\small\mbox{e}}}(t_{n+1}) + s({u_{\small\mbox{e}}}(t_{n+1})) + \\ & (\beta {v_{\small\mbox{e}}}(t_n){v_{\small\mbox{e}}}(t_{n+1}) - \beta {v_{\small\mbox{e}}}(t_{n+1}){v_{\small\mbox{e}}}(t_{n+1}))),\end{split}\]

which reduces to

\[\begin{split}\frac{\beta}{m}{v_{\small\mbox{e}}}(t_{n+1}({v_{\small\mbox{e}}}(t_n) - {v_{\small\mbox{e}}}(t_{n+1})) &= \frac{\beta}{m}{v_{\small\mbox{e}}}(t_{n+1}[D_t^-{v_{\small\mbox{e}}}]^{n+1}\Delta t\\ &= \frac{\beta}{m}{v_{\small\mbox{e}}}(t_{n+1}({v_{\small\mbox{e}}}'(t_{n+1})\Delta t + -\frac{1}{2}{v_{\small\mbox{e}}}'''(t_{n+1})\Delta t^ + {\mathcal{O}(\Delta t^3)}){\thinspace .}\end{split}\]

We end with \(R_u^n\) and \(R_v^{n+1}\) as \({\mathcal{O}(\Delta t)}\), simply because all the building blocks in the schemes (the forward and backward differences and the linearization trick) are only first-order accurate. However, this analysis is misleading: the building blocks play together in a way that makes the scheme second-order accurate. This is shown by considering an alternative, yet equivalent, formulation of the above scheme.

A centered scheme on a staggered mesh

We now introduce a staggered mesh where we seek \(u\) at mesh points \(t_n\) and \(v\) at points \(t_{n+\frac{1}{2}}\) in between the \(u\) points. The staggered mesh makes it easy to formulate centered differences in the system (9)-(9):

(19)\[ \lbrack D_t u = v \rbrack^{n-\frac{1}{2}},\]

(20)\[ \lbrack D_t v = \frac{1}{m}( F(t) - \beta |v|v - s(u)) \rbrack^{n}{\thinspace .}\]

The term \(|v^n|v^n\) causes trouble since \(v^n\) is not computed, only \(v^{n-\frac{1}{2}}\) and \(v^{n+\frac{1}{2}}\). Using geometric mean, we can express \(|v^n|v^n\) in terms of known quantities: \(|v^n|v^n \approx |v^{n-\frac{1}{2}}|v^{n+\frac{1}{2}}\). We then have

(21)\[ \lbrack D_t u \rbrack^{n-\frac{1}{2}} = v^{n-\frac{1}{2}},\]

(22)\[ \lbrack D_t v \rbrack^n = \frac{1}{m}( F(t_n) - \beta |v^{n-\frac{1}{2}}|v^{n+\frac{1}{2}} - s(u^n)){\thinspace .}\]

The truncation error in each equation fulfills

\[\begin{split}\lbrack D_t {u_{\small\mbox{e}}} \rbrack^{n-\frac{1}{2}} &= {v_{\small\mbox{e}}}(t_{n-\frac{1}{2}}) + R_u^{n-\frac{1}{2}},\\ \lbrack D_t {v_{\small\mbox{e}}} \rbrack^n &= \frac{1}{m}( F(t_n) - \beta |{v_{\small\mbox{e}}}(t_{n-\frac{1}{2}})|{v_{\small\mbox{e}}}(t_{n+\frac{1}{2}}) - s(u^n)) + R_v^n{\thinspace .}\end{split}\]

The truncation error of the centered differences is given by (1.3)-(1.4), and the geometric mean approximation analysis can be taken from (1.21)-(1.22). These results lead to

\[{u_{\small\mbox{e}}}'(t_{n-\frac{1}{2}}) + \frac{1}{24}{u_{\small\mbox{e}}}'''(t_{n-\frac{1}{2}})\Delta t^2 + {\mathcal{O}(\Delta t^4)} = {v_{\small\mbox{e}}}(t_{n-\frac{1}{2}}) + R_u^{n-\frac{1}{2}},\]

and

\[{v_{\small\mbox{e}}}'(t_n) = \frac{1}{m}( F(t_n) - \beta |{v_{\small\mbox{e}}}(t_n)|{v_{\small\mbox{e}}}(t_n) + {\mathcal{O}(\Delta t^2)} - s(u^n)) + R_v^n{\thinspace .}\]

The ODEs fulfilled by \({u_{\small\mbox{e}}}\) and \({v_{\small\mbox{e}}}\) are evident in these equations, and we achieve second-order accuracy for the truncation error in both equations:

\[R_u^{n-\frac{1}{2}}= {\mathcal{O}(\Delta t^2)}, \quad R_v^n = {\mathcal{O}(\Delta t^2)}{\thinspace .}\]

Comparing (21)-(22) with (13)-(14), we can hopefully realize that these schemes are equivalent (which becomes clear when we implement both). The obvious advantage with the staggered mesh approach is that we can all the way use second-order accurate building blocks and in this way concince ourselves that the resulting scheme has an error of \({\mathcal{O}(\Delta t^2)}\).