$$\newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\acc}{\boldsymbol{a}} \newcommand{\rpos}{\boldsymbol{r}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\q}{\boldsymbol{q}} \newcommand{\g}{\boldsymbol{g}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\stepzero}{*} \newcommand{\stephalf}{***} \newcommand{\stepone}{**} \newcommand{\baspsi}{\psi} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integer}{\mathbb{Z}}$$

# Vibration ODEs

## Linear model without damping

The next example on computing the truncation error involves the following ODE for vibration problems: $$$$u''(t) + \omega^2 u(t) = 0\tp \tag{7.48}$$$$ 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 in time, we have the scheme $$$$[D_tD_t u + \omega^2u=0]^n \tag{7.49} \tp$$$$

Inserting the exact solution $$\uex$$ in this equation and adding a residual $$R$$ so that $$\uex$$ can fulfill the equation results in $$$$[D_tD_t \uex + \omega^2\uex =R]^n \tp \tag{7.50}$$$$ 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 $$$$R^n = \frac{1}{12}\uex''''(t_n)\Delta t^2 + \Oof{\Delta t^4} \tp \tag{7.51}$$$$

### The truncation error of approximating $$u'(0)$$

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.

### 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 $$\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 $$$$\frac{u^1 - u^0}{\Delta t} + \half\omega^2\Delta t u^0 = V\tp \tag{7.52}$$$$ 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 $$$$[ D_{2t} u + \Delta t(D_tD_t u - \omega^2 u) = V]^0\tp \tag{7.53}$$$$ 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$$

### 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 \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 $$$$u'' + \omega^2(1 - \frac{1}{12}\omega^2\Delta t^2)u = 0\tp \tag{7.54}$$$$ 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(


One will see that the rates r lie around 4.

## Model with damping and nonlinearity

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)$$: $$$$mu'' + \beta u' + s(u) =F(t)\tp \tag{7.55}$$$$ 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 \tp \tag{7.56}$$$$ The exact solution $$\uex$$ fulfills the discrete equations with a residual term: $$$$[mD_tD_t \uex + \beta D_{2t} \uex + s(\uex)=F + R]^n \tp \tag{7.57}$$$$ 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 $$$$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}$$$$ 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'$$: $$$$mu'' + \beta |u'|u' + s(u) =F(t)\tp \tag{7.59}$$$$ 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 $$$$[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}$$$$ The truncation error is defined through $$$$[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}$$$$

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 general model formulated as first-order ODEs

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.

### 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$$