.. !split .. _trunc:vib: Vibration ODEs (2) =========================== .. _trunc:vib:undamped: Linear model without damping ---------------------------- The next example on computing the truncation error involves the following ODE for vibration problems: .. _Eq:trunc:vib:undamped:ode: .. math:: \tag{701} u''(t) + \omega^2 u(t) = 0{\thinspace .} Here, :math:`\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 .. _Eq:trunc:vib:undamped:scheme: .. math:: \tag{702} [D_tD_t u + \omega^2u=0]^n {\thinspace .} Inserting the exact solution :math:`{u_{\small\mbox{e}}}` in this equation and adding a residual :math:`R` so that :math:`{u_{\small\mbox{e}}}` can fulfill the equation results in .. _Eq:_auto261: .. math:: \tag{703} [D_tD_t {u_{\small\mbox{e}}} + \omega^2{u_{\small\mbox{e}}} =R]^n {\thinspace .} To calculate the truncation error :math:`R^n`, we use :ref:`(670) `-:ref:`(671) `, i.e., .. math:: [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 + {\mathcal{O}(\Delta t^4)}, and the fact that :math:`{u_{\small\mbox{e}}}''(t) + \omega^2{u_{\small\mbox{e}}}(t)=0`. The result is .. _Eq:_auto262: .. math:: \tag{704} 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 :math:`u'(0)` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The initial conditions for :ref:`(701) ` are :math:`u(0)=I` and :math:`u'(0)=V`. The latter involves a finite difference approximation. The standard choice .. math:: [D_{2t}u=V]^0, where :math:`u^{-1}` is eliminated with the aid of the discretized ODE for :math:`n=0`, involves a centered difference with an :math:`{\mathcal{O}(\Delta t^2)}` truncation error given by :ref:`(660) `-:ref:`(661) `. The simpler choice .. math:: [D_t^+u = V]^0, is based on a forward difference with a truncation error :math:`{\mathcal{O}(\Delta t)}`. A central question is if this initial error will impact the order of the scheme throughout the simulation. :ref:`trunc:exer:vib:ic:fw` 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 :math:`u'(0)=0` is :math:`{\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 :math:`[D_{2t}u=V]^0` to express :math:`u^{-1}=u^1 - 2\Delta t V` in order to eliminate :math:`u^{-1}` from the discretized differential equation, but the other way around: the fundamental equation is the discretized initial condition :math:`[D_{2t}u=V]^0` and we use the discretized ODE :math:`[D_tD_t + \omega^2 u=0]^0` to eliminate :math:`u^{-1}` in the discretized initial condition. From :math:`[D_tD_t + \omega^2 u=0]^0` we have .. math:: u^{-1} = 2u^0 - u^1 - \Delta t^2\omega^2 u^0, which inserted in :math:`[D_{2t}u = V]^0` gives .. _Eq:trunc:vib:undamped:ic:d2: .. math:: \tag{705} \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 .. math:: [D_t^+ u + \frac{1}{2}\omega^2\Delta t u = V]^0{\thinspace .} The truncation error is defined as .. math:: [D_t^+ {u_{\small\mbox{e}}} + \frac{1}{2}\omega^2\Delta t {u_{\small\mbox{e}}} - V = R]^0{\thinspace .} Using :ref:`(664) `-:ref:`(665) ` with one more term in the Taylor series, we get that .. math:: {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, :math:`{u_{\small\mbox{e}}}'(0)=V` and :math:`{u_{\small\mbox{e}}}''(0)=-\omega^2 {u_{\small\mbox{e}}}(0)` so we get .. math:: 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 :math:`u^{-1}` via the discretized ODE can be expressed as .. _Eq:trunc:vib:undamped:ic:d3: .. math:: \tag{706} [ D_{2t} u + \Delta t(D_tD_t u - \omega^2 u) = V]^0{\thinspace .} Writing out :ref:`(706) ` shows that the equation is equivalent to :ref:`(705) `. The truncation error is defined by .. math:: [ 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 :ref:`(660) `-:ref:`(661) ` and :ref:`(670) `-:ref:`(671) `, as well as using :math:`{u_{\small\mbox{e}}}'(0)=V` and :math:`{u_{\small\mbox{e}}}''(0) = -\omega^2{u_{\small\mbox{e}}}(0)`, gives .. math:: R^n = \frac{1}{6}{u_{\small\mbox{e}}}'''(0)\Delta t^2 + {\mathcal{O}(\Delta t^3)}{\thinspace .} Computing correction terms ~~~~~~~~~~~~~~~~~~~~~~~~~~ The idea of using correction terms to increase the order of :math:`R^n` can be applied as described in the section :ref:`trunc:decay:corr`. We look at .. math:: [D_tD_t {u_{\small\mbox{e}}} + \omega^2{u_{\small\mbox{e}}} =C + R]^n, and observe that :math:`C^n` must be chosen to cancel the :math:`\Delta t^2` term in :math:`R^n`. That is, .. math:: 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: :math:`u''=-\omega^2u`, which implies :math:`u'''' = \omega^4 u`. Adding the correction term to the ODE results in .. _Eq:trunc:vib:undamped:corr:ode: .. math:: \tag{707} u'' + \omega^2(1 - \frac{1}{12}\omega^2\Delta t^2)u = 0{\thinspace .} Solving this equation by the standard scheme .. math:: [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 :math:`{\mathcal{O}(\Delta t^4)}`. We can use another set of arguments to justify that :ref:`(707) ` leads to a higher-order method. Mathematical analysis of the scheme :ref:`(702) ` reveals that the numerical frequency :math:`\tilde\omega` is (approximately as :math:`\Delta t\rightarrow 0`) .. math:: \tilde\omega = \omega (1+\frac{1}{24}\omega^2\Delta t^2){\thinspace .} One can therefore attempt to replace :math:`\omega` in the ODE by a slightly smaller :math:`\omega` since the numerics will make it larger: .. Ref to exercise .. math:: [ 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 :math:`\Delta t^4` gives exactly the ODE :ref:`(707) `. Experiments show that :math:`u^n` is computed to 4th order in :math:`\Delta t`. You can confirm this by running a little program in the ``vib`` directory: .. code-block:: python 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. .. _trunc:vib:gen: Model with damping and nonlinearity ----------------------------------- The model :ref:`(701) ` can be extended to include damping :math:`\beta u'`, a nonlinear restoring (spring) force :math:`s(u)`, and some known excitation force :math:`F(t)`: .. _Eq:trunc:vib:gen:ode1: .. math:: \tag{708} mu'' + \beta u' + s(u) =F(t){\thinspace .} The coefficient :math:`m` usually represents the mass of the system. This governing equation can by discretized by centered differences: .. _Eq:_auto263: .. math:: \tag{709} [mD_tD_t u + \beta D_{2t} u + s(u)=F]^n {\thinspace .} The exact solution :math:`{u_{\small\mbox{e}}}` fulfills the discrete equations with a residual term: .. _Eq:_auto264: .. math:: \tag{710} [mD_tD_t {u_{\small\mbox{e}}} + \beta D_{2t} {u_{\small\mbox{e}}} + s({u_{\small\mbox{e}}})=F + R]^n {\thinspace .} Using :ref:`(670) `-:ref:`(671) ` and :ref:`(660) `-:ref:`(661) ` we get .. math:: \begin{align*} \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{align*} Combining this with the previous equation, we can collect the terms .. math:: 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 :math:`{u_{\small\mbox{e}}}` solves the differential equation. We are left with the truncation error .. _Eq:trunc:vib:gen:R: .. math:: \tag{711} 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 :ref:`(711) `, we can add correction terms .. math:: 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 :math:`u''''` and :math:`u'''` in terms of lower-order derivatives is now harder because the differential equation is more complicated: .. math:: \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''){\thinspace .} \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 :math:`\Delta t` is usually always possible in these problems to achieve the desired accuracy. Extension to quadratic damping ------------------------------ Instead of the linear damping term :math:`\beta u'` in :ref:`(708) ` we now consider quadratic damping :math:`\beta |u'|u'`: .. _Eq:trunc:vib:gen:ode2: .. math:: \tag{712} mu'' + \beta |u'|u' + s(u) =F(t){\thinspace .} A centered difference for :math:`u'` gives rise to a nonlinearity, which can be linearized using a geometric mean: :math:`[|u'|u']^n \approx |[u']^{n-\frac{1}{2}}|[u']^{n+\frac{1}{2}}`. The resulting scheme becomes .. _Eq:_auto265: .. math:: \tag{713} [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 .. _Eq:_auto266: .. math:: \tag{714} [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 :ref:`(676) `-:ref:`(677) `, .. math:: \begin{align*} |[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_{\small\mbox{e}}}'(t_n)^2\Delta t^2 +\\ &\quad \frac{1}{4}{u_{\small\mbox{e}}}(t_n){u_{\small\mbox{e}}}''(t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)}{\thinspace .} \end{align*} Using :ref:`(658) `-:ref:`(659) ` for the :math:`D_t{u_{\small\mbox{e}}}` factors results in .. math:: [|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 .. math:: [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 .. math:: 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 :math:`mu'' + \beta (u')^2 + s(u)=F`, we end up with .. math:: 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 :math:`{\mathcal{O}(\Delta t^2)}` *and* the geometric mean approximation is also :math:`{\mathcal{O}(\Delta t^2)}`. .. _trunc:vib:gen:staggered: The general model formulated as first-order ODEs ------------------------------------------------ The second-order model :ref:`(712) ` can be formulated as a first-order system, .. _Eq:trunc:vib:gen:2x2model:ode:v: .. math:: \tag{715} v' = \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right), .. _Eq:trunc:vib:gen:2x2model:ode:u: .. math:: \tag{716} u' = v{\thinspace .} The system :ref:`(716) `-:ref:`(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 :math:`u` at mesh points :math:`t_n` and :math:`v` at points :math:`t_{n+\frac{1}{2}}` in between the :math:`u` points. The staggered mesh makes it easy to formulate centered differences in the system :ref:`(716) `-:ref:`(716) `: .. _Eq:trunc:vib:gen:2x2model:ode:u:staggered: .. math:: \tag{717} \lbrack D_t u = v \rbrack^{n-\frac{1}{2}}, .. _Eq:trunc:vib:gen:2x2model:ode:v:staggered: .. math:: \tag{718} \lbrack D_t v = \frac{1}{m}( F(t) - \beta |v|v - s(u)) \rbrack^{n}{\thinspace .} The term :math:`|v^n|v^n` causes trouble since :math:`v^n` is not computed, only :math:`v^{n-\frac{1}{2}}` and :math:`v^{n+\frac{1}{2}}`. Using geometric mean, we can express :math:`|v^n|v^n` in terms of known quantities: :math:`|v^n|v^n \approx |v^{n-\frac{1}{2}}|v^{n+\frac{1}{2}}`. We then have .. _Eq:trunc:vib:gen:2x2model:ode:u:staggered2: .. math:: \tag{719} \lbrack D_t u \rbrack^{n-\frac{1}{2}} = v^{n-\frac{1}{2}}, .. _Eq:trunc:vib:gen:2x2model:ode:v:staggered2: .. math:: \tag{720} \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 .. math:: \begin{align*} \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{align*} The truncation error of the centered differences is given by :ref:`(658) `-:ref:`(659) `, and the geometric mean approximation analysis can be taken from :ref:`(676) `-:ref:`(677) `. These results lead to .. math:: {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 .. math:: {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 :math:`{u_{\small\mbox{e}}}` and :math:`{v_{\small\mbox{e}}}` are evident in these equations, and we achieve second-order accuracy for the truncation error in both equations: .. math:: R_u^{n-\frac{1}{2}}= {\mathcal{O}(\Delta t^2)}, \quad R_v^n = {\mathcal{O}(\Delta t^2)}{\thinspace .}