Vibration ODEs¶
Linear model without damping¶
The next example on computing the truncation error involves the following ODE for vibration problems:
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
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
To calculate the truncation error \(R^n\), we use (670)-(671), i.e.,
and the fact that \({u_{\small\mbox{e}}}''(t) + \omega^2{u_{\small\mbox{e}}}(t)=0\). The result is
The truncation error of approximating \(u'(0)\)¶
The initial conditions for (701) are \(u(0)=I\) and \(u'(0)=V\). The latter involves a finite difference approximation. The standard choice
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 (660)-(661). The simpler choice
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 B.11: Investigate the impact of approximating 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 \({\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 discretized initial condition. From \([D_tD_t + \omega^2 u=0]^0\) we have
which inserted in \([D_{2t}u = V]^0\) gives
The first term can be recognized as a forward difference such that the equation can be written in operator notation as
The truncation error is defined as
Using (664)-(665) with one more term in the Taylor series, we get that
Now, \({u_{\small\mbox{e}}}'(0)=V\) and \({u_{\small\mbox{e}}}''(0)=-\omega^2 {u_{\small\mbox{e}}}(0)\) so we get
There is another way of analyzing the discrete initial condition, because eliminating \(u^{-1}\) via the discretized ODE can be expressed as
Writing out (706) shows that the equation is equivalent to (705). The truncation error is defined by
Replacing the difference via (660)-(661) and (670)-(671), as well as using \({u_{\small\mbox{e}}}'(0)=V\) and \({u_{\small\mbox{e}}}''(0) = -\omega^2{u_{\small\mbox{e}}}(0)\), gives
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
and observe that \(C^n\) must be chosen to cancel the \(\Delta t^2\) term in \(R^n\). That is,
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
Solving this equation by the standard scheme
will result in a scheme with truncation error \({\mathcal{O}(\Delta t^4)}\).
We can use another set of arguments to justify that (707) leads to a higher-order method. Mathematical analysis of the scheme (702) reveals that the numerical frequency \(\tilde\omega\) is (approximately as \(\Delta t\rightarrow 0\))
One can therefore attempt to replace \(\omega\) in the ODE by a slightly smaller \(\omega\) since the numerics will make it larger:
Expanding the squared term and omitting the higher-order term \(\Delta t^4\)
gives exactly the ODE (707). 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.
Model with damping and nonlinearity¶
The model (701) can be extended to include damping \(\beta u'\), a nonlinear restoring (spring) force \(s(u)\), and some known excitation force \(F(t)\):
The coefficient \(m\) usually represents the mass of the system. This governing equation can by discretized by centered differences:
The exact solution \({u_{\small\mbox{e}}}\) fulfills the discrete equations with a residual term:
Using (670)-(671) and (660)-(661) we get
Combining this with the previous equation, we can collect the terms
and set this sum to zero because \({u_{\small\mbox{e}}}\) solves the differential equation. We are left with the truncation error
so the scheme is of second order.
According to (711), we can add correction terms
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:
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 (708) we now consider quadratic damping \(\beta |u'|u'\):
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
The truncation error is defined through
We start with expressing the truncation error of the geometric mean. According to (676)-(677),
Using (658)-(659) for the \(D_t{u_{\small\mbox{e}}}\) factors results in
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
With
and using the differential equation on the form \(mu'' + \beta (u')^2 + s(u)=F\), we end up with
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 \({\mathcal{O}(\Delta t^2)}\).
The general model formulated as first-order ODEs¶
The second-order model (712) can be formulated as a first-order system,
The system (716)-(716) 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+\frac{1}{2}}\) in between the \(u\) points. The staggered mesh makes it easy to formulate centered differences in the system (716)-(716):
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
The truncation error in each equation fulfills
The truncation error of the centered differences is given by (658)-(659), and the geometric mean approximation analysis can be taken from (676)-(677). These results lead to
and
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:
