The next example on computing the truncation error involves the following ODE for vibration problems: $$ \begin{equation} u''(t) + \omega^2 u(t) = 0\tp \tag{7.48} \end{equation} $$ Here, \( \omega \) is a given constant.
Using a standard, second-ordered, central difference for the second-order derivative in time, we have the scheme $$ \begin{equation} [D_tD_t u + \omega^2u=0]^n \tag{7.49} \tp \end{equation} $$
Inserting the exact solution \( \uex \) in this equation and adding a residual \( R \) so that \( \uex \) can fulfill the equation results in $$ \begin{equation} [D_tD_t \uex + \omega^2\uex =R]^n \tp \tag{7.50} \end{equation} $$ To calculate the truncation error \( R^n \), we use (7.17)-(7.18), i.e., $$ [D_tD_t \uex]^n = \uex''(t_n) + \frac{1}{12}\uex''''(t_n)\Delta t^2 + \Oof{\Delta t^4},$$ and the fact that \( \uex''(t) + \omega^2\uex(t)=0 \). The result is $$ \begin{equation} R^n = \frac{1}{12}\uex''''(t_n)\Delta t^2 + \Oof{\Delta t^4} \tp \tag{7.51} \end{equation} $$
The initial conditions for (7.48) 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 \( \Oof{\Delta t^2} \) truncation error given by (7.7)-(7.8). The simpler choice $$ [D_t^+u = V]^0,$$ is based on a forward difference with a truncation error \( \Oof{\Delta t} \). A central question is if this initial error will impact the order of the scheme throughout the simulation. Exercise B.11: Investigate the impact of approximating \( u'(0) \) asks you to perform an experiment to investigate this question.
We have shown that the truncation error of the difference used to approximate the initial condition \( u'(0)=0 \) is \( \Oof{\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 discretized 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 $$ \begin{equation} \frac{u^1 - u^0}{\Delta t} + \half\omega^2\Delta t u^0 = V\tp \tag{7.52} \end{equation} $$ The first term can be recognized as a forward difference such that the equation can be written in operator notation as $$ [D_t^+ u + \half\omega^2\Delta t u = V]^0\tp$$ The truncation error is defined as $$ [D_t^+ \uex + \half\omega^2\Delta t \uex - V = R]^0\tp$$ Using (7.11)-(7.12) with one more term in the Taylor series, we get that $$ \uex'(0) + \half\uex''(0)\Delta t + \frac{1}{6}\uex'''(0)\Delta t^2 + \Oof{\Delta t^3} + \half\omega^2\Delta t \uex(0) - V = R^n\tp$$ Now, \( \uex'(0)=V \) and \( \uex''(0)=-\omega^2 \uex(0) \) so we get $$ R^n = \frac{1}{6}\uex'''(0)\Delta t^2 + \Oof{\Delta t^3}\tp$$
There is another way of analyzing the discrete initial condition, because eliminating \( u^{-1} \) via the discretized ODE can be expressed as $$ \begin{equation} [ D_{2t} u + \Delta t(D_tD_t u - \omega^2 u) = V]^0\tp \tag{7.53} \end{equation} $$ Writing out (7.53) shows that the equation is equivalent to (7.52). The truncation error is defined by $$ [ D_{2t} \uex + \Delta t(D_tD_t \uex - \omega^2 \uex) = V + R]^0\tp$$ Replacing the difference via (7.7)-(7.8) and (7.17)-(7.18), as well as using \( \uex'(0)=V \) and \( \uex''(0) = -\omega^2\uex(0) \), gives $$ R^n = \frac{1}{6}\uex'''(0)\Delta t^2 + \Oof{\Delta t^3}\tp$$
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 \uex + \omega^2\uex =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}\uex''''(t_n)\Delta t^2\tp$$ To get rid of the 4th-order derivative we can use the differential equation: \( u''=-\omega^2u \), which implies \( u'''' = \omega^4 u \). Adding the correction term to the ODE results in $$ \begin{equation} u'' + \omega^2(1 - \frac{1}{12}\omega^2\Delta t^2)u = 0\tp \tag{7.54} \end{equation} $$ 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 truncation error \( \Oof{\Delta t^4} \).
We can use another set of arguments to justify that (7.54) leads to a higher-order method.
Mathematical analysis of the scheme (7.49)
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)\tp$$
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\tp$$
Expanding the squared term and omitting the higher-order term \( \Delta t^4 \)
gives exactly the ODE (7.54). Experiments
show that \( u^n \) is computed to 4th order in \( \Delta t \). You can confirm
this by running a little program in the vib directory:
from vib_undamped import convergence_rates, solver_adjust_w
r = convergence_rates(
m=5, solver_function=solver_adjust_w, num_periods=8)
One will see that the rates r lie around 4.
The model (7.48) can be extended to include damping \( \beta u' \), a nonlinear restoring (spring) force \( s(u) \), and some known excitation force \( F(t) \): $$ \begin{equation} mu'' + \beta u' + s(u) =F(t)\tp \tag{7.55} \end{equation} $$ The coefficient \( m \) usually represents the mass of the system. This governing equation can by discretized by centered differences: $$ \begin{equation} [mD_tD_t u + \beta D_{2t} u + s(u)=F]^n \tp \tag{7.56} \end{equation} $$ The exact solution \( \uex \) fulfills the discrete equations with a residual term: $$ \begin{equation} [mD_tD_t \uex + \beta D_{2t} \uex + s(\uex)=F + R]^n \tp \tag{7.57} \end{equation} $$ Using (7.17)-(7.18) and (7.7)-(7.8) 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*} $$ Combining this with the previous equation, we can collect the terms $$ m\uex''(t_n) + \beta\uex'(t_n) + \omega^2\uex(t_n) + s(\uex(t_n)) - F^n,$$ and set this sum to zero because \( \uex \) solves the differential equation. We are left with the truncation error $$ \begin{equation} R^n = \left(\frac{m}{12}\uex''''(t_n) + \frac{\beta}{6}\uex'''(t_n)\right)\Delta t^2 + \Oof{\Delta t^4}, \tag{7.58} \end{equation} $$ so the scheme is of second order.
According to (7.58), we can add correction terms $$ C^n = \left(\frac{m}{12}\uex''''(t_n) + \frac{\beta}{6}\uex'''(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{align*} 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'')\tp \end{align*} $$ 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.
Instead of the linear damping term \( \beta u' \) in (7.55) we now consider quadratic damping \( \beta |u'|u' \): $$ \begin{equation} mu'' + \beta |u'|u' + s(u) =F(t)\tp \tag{7.59} \end{equation} $$ A centered difference for \( u' \) gives rise to a nonlinearity, which can be linearized using a geometric mean: \( [|u'|u']^n \approx |[u']^{n-\half}|[u']^{n+\half} \). The resulting scheme becomes $$ \begin{equation} [mD_t D_t u]^n + \beta |[D_{t} u]^{n-\half}|[D_t u]^{n+\half} + s(u^n)=F^n\tp \tag{7.60} \end{equation} $$ The truncation error is defined through $$ \begin{equation} [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 \tag{7.61} \end{equation} $$
We start with expressing the truncation error of the geometric mean. According to (7.23)-(7.24), $$ \begin{align*} |[D_{t} \uex]^{n-\half}|[D_t \uex]^{n+\half} &= [|D_t\uex|D_t\uex]^n - \frac{1}{4}\uex'(t_n)^2\Delta t^2 +\\ &\quad \frac{1}{4}\uex(t_n)\uex''(t_n)\Delta t^2 + \Oof{\Delta t^4}\tp \end{align*} $$ Using (7.5)-(7.6) for the \( D_t\uex \) factors results in $$ [|D_t\uex|D_t\uex]^n = |\uex' + \frac{1}{24}\uex'''(t_n)\Delta t^2 + \Oof{\Delta t^4}|(\uex' + \frac{1}{24}\uex'''(t_n)\Delta t^2 + \Oof{\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\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 the differential equation on the form \( 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$$ 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 \( \Oof{\Delta t^2} \) and the geometric mean approximation is also \( \Oof{\Delta t^2} \).
The second-order model (7.59) can be formulated as a first-order system, $$ \begin{align} v' &= \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right), \tag{7.62}\\ u' &= v\tp \tag{7.63} \end{align} $$ The system (7.63)-(7.63) can be solved either by a forward-backward scheme (the Euler-Cromer method) or 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+\half} \) in between the \( u \) points. The staggered mesh makes it easy to formulate centered differences in the system (7.63)-(7.63): $$ \begin{align} \lbrack D_t u &= v \rbrack^{n-\half}, \tag{7.64} \\ \lbrack D_t v &= \frac{1}{m}( F(t) - \beta |v|v - s(u)) \rbrack^{n}\tp \tag{7.65} \end{align} $$ The term \( |v^n|v^n \) causes trouble since \( v^n \) is not computed, only \( v^{n-\half} \) and \( v^{n+\half} \). Using geometric mean, we can express \( |v^n|v^n \) in terms of known quantities: \( |v^n|v^n \approx |v^{n-\half}|v^{n+\half} \). We then have $$ \begin{align} \lbrack D_t u \rbrack^{n-\half} &= v^{n-\half}, \tag{7.66} \\ \lbrack D_t v \rbrack^n &= \frac{1}{m}( F(t_n) - \beta |v^{n-\half}|v^{n+\half} - s(u^n))\tp \tag{7.67} \end{align} $$ The truncation error in each equation fulfills $$ \begin{align*} \lbrack D_t \uex \rbrack^{n-\half} &= \vex(t_{n-\half}) + R_u^{n-\half},\\ \lbrack D_t \vex \rbrack^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*} $$ The truncation error of the centered differences is given by (7.5)-(7.6), and the geometric mean approximation analysis can be taken from (7.23)-(7.24). These results lead to $$ \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 $$ The ODEs fulfilled by \( \uex \) and \( \vex \) are evident in these equations, and we achieve second-order accuracy for the truncation error in both equations: $$ R_u^{n-\half}= \Oof{\Delta t^2}, \quad R_v^n = \Oof{\Delta t^2}\tp$$