.. !split .. _vib:model2: Generalization: damping, nonlinearities, and excitation ======================================================= .. index:: nonlinear restoring force .. index:: nonlinear spring .. index:: forced vibrations We shall now generalize the simple model problem from the section :ref:`vib:model1` to include a possibly nonlinear damping term :math:`f(u^{\prime})`, a possibly nonlinear spring (or restoring) force :math:`s(u)`, and some external excitation :math:`F(t)`: .. _Eq:vib:ode2: .. math:: \tag{71} mu^{\prime\prime} + f(u^{\prime}) + s(u) = F(t),\quad u(0)=I,\ u^{\prime}(0)=V,\ t\in (0,T] {\thinspace .} We have also included a possibly nonzero initial value of :math:`u^{\prime}(0)`. The parameters :math:`m`, :math:`f(u^{\prime})`, :math:`s(u)`, :math:`F(t)`, :math:`I`, :math:`V`, and :math:`T` are input data. There are two main types of damping (friction) forces: linear :math:`f(u^{\prime})=bu`, or quadratic :math:`f(u^{\prime})=bu^{\prime}|u^{\prime}|`. Spring systems often feature linear damping, while air resistance usually gives rise to quadratic damping. Spring forces are often linear: :math:`s(u)=cu`, but nonlinear versions are also common, the most famous is the gravity force on a pendulum that acts as a spring with :math:`s(u)\sim \sin(u)`. .. _vib:ode2:fdm:flin: A centered scheme for linear damping ------------------------------------ Sampling :ref:`(71) ` at a mesh point :math:`t_n`, replacing :math:`u^{\prime\prime}(t_n)` by :math:`[D_tD_tu]^n`, and :math:`u^{\prime}(t_n)` by :math:`[D_{2t}u]^n` results in the discretization .. _Eq:_auto32: .. math:: \tag{72} [mD_tD_t u + f(D_{2t}u) + s(u) = F]^n, which written out means .. _Eq:vib:ode2:step3b: .. math:: \tag{73} m\frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} + f(\frac{u^{n+1}-u^{n-1}}{2\Delta t}) + s(u^n) = F^n, where :math:`F^n` as usual means :math:`F(t)` evaluated at :math:`t=t_n`. Solving :ref:`(73) ` with respect to the unknown :math:`u^{n+1}` gives a problem: the :math:`u^{n+1}` inside the :math:`f` function makes the equation *nonlinear* unless :math:`f(u^{\prime})` is a linear function, :math:`f(u^{\prime})=bu^{\prime}`. For now we shall assume that :math:`f` is linear in :math:`u^{\prime}`. Then .. _Eq:vib:ode2:step3b2: .. math:: \tag{74} m\frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} + b\frac{u^{n+1}-u^{n-1}}{2\Delta t} + s(u^n) = F^n, which gives an explicit formula for :math:`u` at each new time level: .. _Eq:vib:ode2:step4: .. math:: \tag{75} u^{n+1} = (2mu^n + (\frac{b}{2}\Delta t - m)u^{n-1} + \Delta t^2(F^n - s(u^n)))(m + \frac{b}{2}\Delta t)^{-1} {\thinspace .} For the first time step we need to discretize :math:`u^{\prime}(0)=V` as :math:`[D_{2t}u = V]^0` and combine with :ref:`(75) ` for :math:`n=0`. The discretized initial condition leads to .. _Eq:vib:ode2:ic:du: .. math:: \tag{76} u^{-1} = u^{1} - 2\Delta t V, which inserted in :ref:`(75) ` for :math:`n=0` gives an equation that can be solved for :math:`u^1`: .. _Eq:vib:ode2:step4b: .. math:: \tag{77} u^1 = u^0 + \Delta t\, V + \frac{\Delta t^2}{2m}(-bV - s(u^0) + F^0) {\thinspace .} .. _vib:ode2:fdm:fquad: A centered scheme for quadratic damping --------------------------------------- When :math:`f(u^{\prime})=bu^{\prime}|u^{\prime}|`, we get a quadratic equation for :math:`u^{n+1}` in :ref:`(73) `. This equation can be straightforwardly solved by the well-known formula for the roots of a quadratic equation. However, we can also avoid the nonlinearity by introducing an approximation with an error of order no higher than what we already have from replacing derivatives with finite differences. .. index:: geometric mean .. index:: single: averaging; geometric We start with :ref:`(71) ` and only replace :math:`u^{\prime\prime}` by :math:`D_tD_tu`, resulting in .. _Eq:vib:ode2:quad:idea1: .. math:: \tag{78} [mD_tD_t u + bu^{\prime}|u^{\prime}| + s(u) = F]^n{\thinspace .} Here, :math:`u^{\prime}|u^{\prime}|` is to be computed at time :math:`t_n`. The idea is now to introduce a *geometric mean*, defined by .. math:: (w^2)^n \approx w^{n-\frac{1}{2}}w^{n+\frac{1}{2}}, for some quantity :math:`w` depending on time. The error in the geometric mean approximation is :math:`{\mathcal{O}(\Delta t^2)}`, the same as in the approximation :math:`u^{\prime\prime}\approx D_tD_tu`. With :math:`w=u^{\prime}` it follows that .. math:: [u^{\prime}|u^{\prime}|]^n \approx u^{\prime}(t_{n+\frac{1}{2}})|u^{\prime}(t_{n-\frac{1}{2}})|{\thinspace .} The next step is to approximate :math:`u^{\prime}` at :math:`t_{n\pm 1/2}`, and fortunately a centered difference fits perfectly into the formulas since it involves :math:`u` values at the mesh points only. With the approximations .. _Eq:vib:ode2:quad:idea2: .. math:: \tag{79} u^{\prime}(t_{n+1/2})\approx [D_t u]^{n+\frac{1}{2}},\quad u^{\prime}(t_{n-1/2})\approx [D_t u]^{n-\frac{1}{2}}, we get .. _Eq:_auto33: .. math:: \tag{80} [u^{\prime}|u^{\prime}|]^n \approx [D_tu]^{n+\frac{1}{2}}|[D_tu]^{n-\frac{1}{2}}| = \frac{u^{n+1}-u^n}{\Delta t} \frac{|u^n-u^{n-1}|}{\Delta t} {\thinspace .} The counterpart to :ref:`(73) ` is then .. _Eq:vib:ode2:step3b:quad: .. math:: \tag{81} m\frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} + b\frac{u^{n+1}-u^n}{\Delta t}\frac{|u^n-u^{n-1}|}{\Delta t} + s(u^n) = F^n, which is linear in the unknown :math:`u^{n+1}`. Therefore, we can easily solve :ref:`(81) ` with respect to :math:`u^{n+1}` and achieve the explicit updating formula .. math:: u^{n+1} = \left( m + b|u^n-u^{n-1}|\right)^{-1}\times \nonumber .. _Eq:vib:ode2:step4:quad: .. math:: \tag{82} \qquad \left(2m u^n - mu^{n-1} + bu^n|u^n-u^{n-1}| + \Delta t^2 (F^n - s(u^n)) \right) {\thinspace .} .. Make exercise to solve complicated u^1 equation with Bisection/Newton In the derivation of a special equation for the first time step we run into some trouble: inserting :ref:`(76) ` in :ref:`(82) ` for :math:`n=0` results in a complicated nonlinear equation for :math:`u^1`. By thinking differently about the problem we can easily get away with the nonlinearity again. We have for :math:`n=0` that :math:`b[u^{\prime}|u^{\prime}|]^0 = bV|V|`. Using this value in :ref:`(78) ` gives .. _Eq:_auto34: .. math:: \tag{83} [mD_tD_t u + bV|V| + s(u) = F]^0 {\thinspace .} Writing this equation out and using :ref:`(76) ` results in the special equation for the first time step: .. _Eq:vib:ode2:step4b:quad: .. math:: \tag{84} u^1 = u^0 + \Delta t V + \frac{\Delta t^2}{2m}\left(-bV|V| - s(u^0) + F^0\right) {\thinspace .} A forward-backward discretization of the quadratic damping term --------------------------------------------------------------- The previous section first proposed to discretize the quadratic damping term :math:`|u^{\prime}|u^{\prime}` using centered differences: :math:`[|D_{2t}|D_{2t}u]^n`. As this gives rise to a nonlinearity in :math:`u^{n+1}`, it was instead proposed to use a geometric mean combined with centered differences. But there are other alternatives. To get rid of the nonlinearity in :math:`[|D_{2t}|D_{2t}u]^n`, one can think differently: apply a backward difference to :math:`|u^{\prime}|`, such that the term involves known values, and apply a forward difference to :math:`u^{\prime}` to make the term linear in the unknown :math:`u^{n+1}`. With mathematics, .. _Eq:vib:ode2:nonlin:fbdiff: .. math:: \tag{85} [\beta |u^{\prime}|u^{\prime}]^n \approx \beta |[D_t^-u]^n|[D_t^+ u]^n = \beta\left\vert\frac{u^n-u^{n-1}}{\Delta t}\right\vert \frac{u^{n+1}-u^n}{\Delta t}{\thinspace .} The forward and backward differences have both an error proportional to :math:`\Delta t` so one may think the discretization above leads to a first-order scheme. However, by looking at the formulas, we realize that the forward-backward differences in :ref:`(85) ` result in exactly the same scheme as in :ref:`(81) ` where we used a geometric mean and centered differences and committed errors of size :math:`{\mathcal{O}(\Delta t^2)}`. Therefore, the forward-backward differences in :ref:`(85) ` act in a symmetric way and actually produce a second-order accurate discretization of the quadratic damping term. .. _vib:ode2:solver: Implementation (3) --------------------------- The algorithm arising from the methods in the sections :ref:`vib:ode2:fdm:flin` and :ref:`vib:ode2:fdm:fquad` is very similar to the undamped case in the section :ref:`vib:ode1:fdm`. The difference is basically a question of different formulas for :math:`u^1` and :math:`u^{n+1}`. This is actually quite remarkable. The equation :ref:`(71) ` is normally impossible to solve by pen and paper, but possible for some special choices of :math:`F`, :math:`s`, and :math:`f`. On the contrary, the complexity of the nonlinear generalized model :ref:`(71) ` versus the simple undamped model is not a big deal when we solve the problem numerically! The computational algorithm takes the form 1. :math:`u^0=I` 2. compute :math:`u^1` from :ref:`(77) ` if linear damping or :ref:`(84) ` if quadratic damping 3. for :math:`n=1,2,\ldots,N_t-1`: a. compute :math:`u^{n+1}` from :ref:`(75) ` if linear damping or :ref:`(82) ` if quadratic damping Modifying the ``solver`` function for the undamped case is fairly easy, the big difference being many more terms and if tests on the type of damping: .. code-block:: python def solver(I, V, m, b, s, F, dt, T, damping='linear'): """ Solve m*u'' + f(u') + s(u) = F(t) for t in (0,T], u(0)=I and u'(0)=V, by a central finite difference method with time step dt. If damping is 'linear', f(u')=b*u, while if damping is 'quadratic', f(u')=b*u'*abs(u'). F(t) and s(u) are Python functions. """ dt = float(dt); b = float(b); m = float(m) # avoid integer div. Nt = int(round(T/dt)) u = np.zeros(Nt+1) t = np.linspace(0, Nt*dt, Nt+1) u[0] = I if damping == 'linear': u[1] = u[0] + dt*V + dt**2/(2*m)*(-b*V - s(u[0]) + F(t[0])) elif damping == 'quadratic': u[1] = u[0] + dt*V + \ dt**2/(2*m)*(-b*V*abs(V) - s(u[0]) + F(t[0])) for n in range(1, Nt): if damping == 'linear': u[n+1] = (2*m*u[n] + (b*dt/2 - m)*u[n-1] + dt**2*(F(t[n]) - s(u[n])))/(m + b*dt/2) elif damping == 'quadratic': u[n+1] = (2*m*u[n] - m*u[n-1] + b*u[n]*abs(u[n] - u[n-1]) + dt**2*(F(t[n]) - s(u[n])))/\ (m + b*abs(u[n] - u[n-1])) return u, t The complete code resides in the file `vib.py `__. .. _vib:ode2:verify: Verification (3) ------------------------- Constant solution ~~~~~~~~~~~~~~~~~ For debugging and initial verification, a constant solution is often very useful. We choose :math:`{u_{\small\mbox{e}}}(t)=I`, which implies :math:`V=0`. Inserted in the ODE, we get :math:`F(t)=s(I)` for any choice of :math:`f`. Since the discrete derivative of a constant vanishes (in particular, :math:`[D_{2t}I]^n=0`, :math:`[D_tI]^n=0`, and :math:`[D_tD_t I]^n=0`), the constant solution also fulfills the discrete equations. The constant should therefore be reproduced to machine precision. The function ``test_constant`` in ``vib.py`` implements this test. Linear solution ~~~~~~~~~~~~~~~ Now we choose a linear solution: :math:`{u_{\small\mbox{e}}} = ct + d`. The initial condition :math:`u(0)=I` implies :math:`d=I`, and :math:`u^{\prime}(0)=V` forces :math:`c` to be :math:`V`. Inserting :math:`{u_{\small\mbox{e}}}=Vt+I` in the ODE with linear damping results in .. math:: 0 + bV + s(Vt+I) = F(t), while quadratic damping requires the source term .. math:: 0 + b|V|V + s(Vt+I) = F(t){\thinspace .} Since the finite difference approximations used to compute :math:`u^{\prime}` all are exact for a linear function, it turns out that the linear :math:`{u_{\small\mbox{e}}}` is also a solution of the discrete equations. :ref:`vib:exer:verify:gen:linear` asks you to carry out all the details. Quadratic solution ~~~~~~~~~~~~~~~~~~ Choosing :math:`{u_{\small\mbox{e}}} = bt^2 + Vt + I`, with :math:`b` arbitrary, fulfills the initial conditions and fits the ODE if :math:`F` is adjusted properly. The solution also solves the discrete equations with linear damping. However, this quadratic polynomial in :math:`t` does not fulfill the discrete equations in case of quadratic damping, because the geometric mean used in the approximation of this term introduces an error. Doing :ref:`vib:exer:verify:gen:linear` will reveal the details. One can fit :math:`F^n` in the discrete equations such that the quadratic polynomial is reproduced by the numerical method (to machine precision). Catching bugs ~~~~~~~~~~~~~ How good are the constant and quadratic solutions at catching bugs in the implementation? * Use ``m`` instead of ``2*m`` in the denominator of ``u[1]``: constant works, while quadratic fails. * Use ``b*dt`` instead of ``b*dt/2`` in the updating formula for ``u[n+1]`` in case of linear damping: constant and quadratic fail. * Use ``F[n+1]`` instead of ``F[n]`` in case of linear or quadratic damping: constant solution works, quadratic fails. We realize that the constant solution is very useful to catch bugs because of its simplicity (easy to predict what the different terms in the formula should evaluate to), while it seems the quadratic solution is capable of detecting all other types of typos in the scheme (?). This results demonstrates why we focus so much on exact, simple polynomial solutions of the numerical schemes in these writings. .. More: classes, cases with pendulum approx u vs sin(u), .. making UI via parampool .. _vib:ode2:viz: Visualization (1) -------------------------- The functions for visualizations differ significantly from those in the undamped case in the ``vib_undamped.py`` program because, in the present general case, we do not have an exact solution to include in the plots. Moreover, we have no good estimate of the periods of the oscillations as there will be one period determined by the system parameters, essentially the approximate frequency :math:`\sqrt{s'(0)/m}` for linear :math:`s` and small damping, and one period dictated by :math:`F(t)` in case the excitation is periodic. This is, however, nothing that the program can depend on or make use of. Therefore, the user has to specify :math:`T` and the window width to get a plot that moves with the graph and shows the most recent parts of it in long time simulations. The ``vib.py`` code contains several functions for analyzing the time series signal and for visualizing the solutions. .. _vib:ode2:ui: User interface -------------- .. index:: ArgumentParser (Python class) .. index:: argparse (Python module) The ``main`` function is changed substantially from the ``vib_undamped.py`` code, since we need to specify the new data :math:`c`, :math:`s(u)`, and :math:`F(t)`. In addition, we must set :math:`T` and the plot window width (instead of the number of periods we want to simulate as in ``vib_undamped.py``). To figure out whether we can use one plot for the whole time series or if we should follow the most recent part of :math:`u`, we can use the ``plot_empricial_freq_and_amplitude`` function's estimate of the number of local maxima. This number is now returned from the function and used in ``main`` to decide on the visualization technique. .. code-block:: python def main(): import argparse parser = argparse.ArgumentParser() parser.add_argument('--I', type=float, default=1.0) parser.add_argument('--V', type=float, default=0.0) parser.add_argument('--m', type=float, default=1.0) parser.add_argument('--c', type=float, default=0.0) parser.add_argument('--s', type=str, default='u') parser.add_argument('--F', type=str, default='0') parser.add_argument('--dt', type=float, default=0.05) parser.add_argument('--T', type=float, default=140) parser.add_argument('--damping', type=str, default='linear') parser.add_argument('--window_width', type=float, default=30) parser.add_argument('--savefig', action='store_true') a = parser.parse_args() from scitools.std import StringFunction s = StringFunction(a.s, independent_variable='u') F = StringFunction(a.F, independent_variable='t') I, V, m, c, dt, T, window_width, savefig, damping = \ a.I, a.V, a.m, a.c, a.dt, a.T, a.window_width, a.savefig, \ a.damping u, t = solver(I, V, m, c, s, F, dt, T) num_periods = empirical_freq_and_amplitude(u, t) if num_periods <= 15: figure() visualize(u, t) else: visualize_front(u, t, window_width, savefig) show() The program ``vib.py`` contains the above code snippets and can solve the model problem :ref:`(71) `. As a demo of ``vib.py``, we consider the case :math:`I=1`, :math:`V=0`, :math:`m=1`, :math:`c=0.03`, :math:`s(u)=\sin(u)`, :math:`F(t)=3\cos(4t)`, :math:`\Delta t = 0.05`, and :math:`T=140`. The relevant command to run is .. code-block:: text Terminal> python vib.py --s 'sin(u)' --F '3*cos(4*t)' --c 0.03 This results in a `moving window following the function `__ on the screen. Figure :ref:`vib:ode2:fig:demo` shows a part of the time series. .. _vib:ode2:fig:demo: .. figure:: vib_gen_demo.png :width: 600 *Damped oscillator excited by a sinusoidal function* .. _vib:ode2:EulerCromer: The Euler-Cromer scheme for the generalized model ------------------------------------------------- The ideas of the Euler-Cromer method from the section :ref:`vib:model2x2:EulerCromer` carry over to the generalized model. We write :ref:`(71) ` as two equations for :math:`u` and :math:`v=u^{\prime}`. The first equation is taken as the one with :math:`v^{\prime}` on the left-hand side: .. _Eq:vib:ode2:EulerCromer:veq: .. math:: \tag{86} v^{\prime} = \frac{1}{m}(F(t)-s(u)-f(v)), .. _Eq:vib:ode2:EulerCromer:ueq: .. math:: \tag{87} u^{\prime} = v{\thinspace .} The idea is to step :ref:`(86) ` forward using a standard Forward Euler method, while we update :math:`u` from :ref:`(87) ` with a Backward Euler method, utilizing the recent, computed :math:`v^{n+1}` value. In detail, .. _Eq:vib:ode2:EulerCromer:dveq0a: .. math:: \tag{88} \frac{v^{n+1}-v^n}{\Delta t} = \frac{1}{m}(F(t_n)-s(u^n)-f(v^n)), .. _Eq:vib:ode2:EulerCromer:dueq0a: .. math:: \tag{89} \frac{u^{n+1}-u^n}{\Delta t} = v^{n+1}, resulting in the explicit scheme .. _Eq:vib:ode2:EulerCromer:dveq: .. math:: \tag{90} v^{n+1} = v^n + \Delta t\frac{1}{m}(F(t_n)-s(u^n)-f(v^n)), .. _Eq:vib:ode2:EulerCromer:dueq0: .. math:: \tag{91} u^{n+1} = u^n + \Delta t\,v^{n+1}{\thinspace .} We immediately note one very favorable feature of this scheme: all the nonlinearities in :math:`s(u)` and :math:`f(v)` are evaluated at a previous time level. This makes the Euler-Cromer method easier to apply and hence much more convenient than the centered scheme for the second-order ODE :ref:`(71) `. The initial conditions are trivially set as .. _Eq:_auto35: .. math:: \tag{92} v^0 = V, .. _Eq:_auto36: .. math:: \tag{93} u^0 = I{\thinspace .} .. Discussion of which method is best: .. ``_ .. _vib:model2x2:gen:StormerVerlet: The Stoermer-Verlet algorithm for the generalized model ------------------------------------------------------- We can easily apply the ideas from the section :ref:`vib:model2x2:StormerVerlet` to extend that method to the generalized model .. math:: \begin{align*} v^{\prime} &= \frac{1}{m}(F(t)-s(u)-f(v)),\\ u^{\prime} &= v{\thinspace .} \end{align*} However, since the scheme is essentially centered differences for the ODE system on a staggered mesh, we do not go into detail here, but refer to the section :ref:`vib:ode2:staggered`. .. _vib:ode2:staggered: A staggered Euler-Cromer scheme for a generalized model ------------------------------------------------------- The more general model for vibration problems, .. _Eq:_auto37: .. math:: \tag{94} mu'' + f(u') + s(u) = F(t),\quad u(0)=I,\ u'(0)=V,\ t\in (0,T], can be rewritten as a first-order ODE system .. _Eq:vib:ode2:staggered:veq: .. math:: \tag{95} v' = m^{-1}\left(F(t) - f(v) - s(u)\right), .. _Eq:vib:ode2:staggered:ueq: .. math:: \tag{96} u' = v{\thinspace .} It is natural to introduce a staggered mesh (see the section :ref:`vib:model2x2:staggered`) and seek :math:`u` at mesh points :math:`t_n` (the numerical value is denoted by :math:`u^n`) and :math:`v` between mesh points at :math:`t_{n+1/2}` (the numerical value is denoted by :math:`v^{n+\frac{1}{2}}`). A centered difference approximation to :ref:`(96) `-:ref:`(95) ` can then be written in operator notation as .. _Eq:vib:ode2:staggered:dveq: .. math:: \tag{97} \lbrack D_tv = m^{-1}\left(F(t) - f(v) - s(u)\right)\rbrack^n, .. _Eq:vib:ode2:staggered:dueq: .. math:: \tag{98} \lbrack D_t u = v\rbrack^{n+\frac{1}{2}}{\thinspace .} Written out, .. _Eq:vib:ode2:staggered:dveq2: .. math:: \tag{99} \frac{v^{n+\frac{1}{2}} - v^{n-\frac{1}{2}}}{\Delta t} = m^{-1}\left(F^n - f(v^n) - s(u^n)\right), .. _Eq:vib:ode2:staggered:dueq2: .. math:: \tag{100} \frac{u^n - u^{n-1}}{\Delta t} = v^{n+\frac{1}{2}}{\thinspace .} With linear damping, :math:`f(v)=bv`, we can use an arithmetic mean for :math:`f(v^n)`: :math:`f(v^n)\approx = \frac{1}{2}(f(v^{n-\frac{1}{2}}) + f(v^{n+\frac{1}{2}}))`. The system :ref:`(99) `-:ref:`(100) ` can then be solved with respect to the unknowns :math:`u^n` and :math:`v^{n+\frac{1}{2}}`: .. _Eq:vib:ode2:staggered:v:scheme:lin: .. math:: \tag{101} v^{n+\frac{1}{2}} = \left(1 + \frac{b}{2m}\Delta t\right)^{-1}\left( v^{n-\frac{1}{2}} + {\Delta t} m^{-1}\left(F^n - {\frac{1}{2}}f(v^{n-\frac{1}{2}}) - s(u^n)\right)\right), .. _Eq:vib:ode2:staggered:u:scheme:lin: .. math:: \tag{102} u^n = u^{n-1} + {\Delta t}v^{n-\frac{1}{2}}{\thinspace .} In case of quadratic damping, :math:`f(v)=b|v|v`, we can use a geometric mean: :math:`f(v^n)\approx b|v^{n-\frac{1}{2}}|v^{n+\frac{1}{2}}`. Inserting this approximation in :ref:`(99) `-:ref:`(100) ` and solving for the unknowns :math:`u^n` and :math:`v^{n+\frac{1}{2}}` results in .. _Eq:vib:ode2:staggered:v:scheme:quad: .. math:: \tag{103} v^{n+\frac{1}{2}} = (1 + \frac{b}{m}|v^{n-\frac{1}{2}}|\Delta t)^{-1}\left( v^{n-\frac{1}{2}} + {\Delta t} m^{-1}\left(F^n - s(u^n)\right)\right), .. _Eq:vib:ode2:staggered:u:scheme:quad: .. math:: \tag{104} u^n = u^{n-1} + {\Delta t}v^{n-\frac{1}{2}}{\thinspace .} The initial conditions are derived at the end of the section :ref:`vib:model2x2:staggered`: .. _Eq:vib:ode2:staggered:u02: .. math:: \tag{105} u^0 = I, .. _Eq:vib:ode2:staggered:v02: .. math:: \tag{106} v^\frac{1}{2} = V - \frac{1}{2}\Delta t\omega^2I {\thinspace .} .. _vib:ode2:PEFRL: The PEFRL 4th-order accurate algorithm -------------------------------------- A variant of the Euler-Cromer type of algorithm, which provides an error :math:`{\mathcal{O}(\Delta t^4)}` if :math:`f(v)=0`, is called PEFRL [Ref05]_. This algorithm is very well suited for integrating dynamic systems (especially those without damping) over very long time periods. Define .. math:: g(u,v) = \frac{1}{m}(F(t)-s(u)-f(v)){\thinspace .} The algorithm is explicit and features these simple steps: .. _Eq:_auto38: .. math:: \tag{107} u^{n+1,1} = u^n + \xi\Delta t v^n, .. _Eq:_auto39: .. math:: \tag{108} v^{n+1,1} = v^n + \frac{1}{2}(1-2\lambda)\Delta t g(u^{n+1,1},v^n), .. _Eq:_auto40: .. math:: \tag{109} u^{n+1,2} = u^{n+1,1} + \chi\Delta t v^{n+1,1}, .. _Eq:_auto41: .. math:: \tag{110} v^{n+1,2} = v^{n+1,1} + \lambda\Delta t g(u^{n+1,2}, v^{n+1,1}), .. _Eq:_auto42: .. math:: \tag{111} u^{n+1,3} = u^{n+1,2} + (1-2(\chi + \xi))\Delta t v^{n+1,2}, .. _Eq:_auto43: .. math:: \tag{112} v^{n+1,3} = v^{n+1,2} + \lambda\Delta t g(u^{n+1,3}, v^{n+1,2}), .. _Eq:_auto44: .. math:: \tag{113} u^{n+1,4} = u^{n+1,3} + \chi\Delta t v^{n+1,3}, .. _Eq:_auto45: .. math:: \tag{114} v^{n+1} = v^{n+1,3} + \frac{1}{2}(1-2\lambda)\Delta t g(u^{n+1,4},v^{n+1,3}), .. _Eq:_auto46: .. math:: \tag{115} u^{n+1} = u^{n+1,4} + \xi\Delta t v^{n+1} .. Compare with eq 22 in article The parameters :math:`\xi`, :math:`\lambda`, and :math:`\xi` have the values .. _Eq:_auto47: .. math:: \tag{116} \xi = 0.1786178958448091, .. _Eq:_auto48: .. math:: \tag{117} \lambda = -0.2123418310626054, .. _Eq:_auto49: .. math:: \tag{118} \chi = -0.06626458266981849 Exercises and Problems (2) =================================== .. --- begin exercise --- .. _vib:exer:gen:class: Exercise 1.19: Implement the solver via classes ----------------------------------------------- Reimplement the ``vib.py`` program using a class ``Problem`` to hold all the physical parameters of the problem, a class ``Solver`` to hold the numerical parameters and compute the solution, and a class ``Visualizer`` to display the solution. .. --- begin hint in exercise --- **Hint.** Use the ideas and examples for an `ODE model `__ in [Ref02]_. More specifically, make a superclass ``Problem`` for holding the scalar physical parameters of a problem and let subclasses implement the :math:`s(u)` and :math:`F(t)` functions as methods. Try to call up as much existing functionality in ``vib.py`` as possible. .. --- end hint in exercise --- .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``vib_class``. .. --- end exercise --- .. --- begin exercise --- .. _vib:exer:quad:damping:bw: Problem 1.20: Use a backward difference for the damping term ------------------------------------------------------------ As an alternative to discretizing the damping terms :math:`\beta u^{\prime}` and :math:`\beta |u^{\prime}|u^{\prime}` by centered differences, we may apply backward differences: .. math:: \begin{align*} [u^{\prime}]^n &\approx [D_t^-u]^n,\\ & [|u^{\prime}|u^{\prime}]^n\\ &\approx [|D_t^-u|D_t^-u]^n\\ &= |[D_t^-u]^n|[D_t^-u]^n{\thinspace .} \end{align*} The advantage of the backward difference is that the damping term is evaluated using known values :math:`u^n` and :math:`u^{n-1}` only. Extend the `vib.py `__ code with a scheme based on using backward differences in the damping terms. Add statements to compare the original approach with centered difference and the new idea launched in this exercise. Perform numerical experiments to investigate how much accuracy that is lost by using the backward differences. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``vib_gen_bwdamping``. .. --- end exercise --- .. --- begin exercise --- .. _vib:exer:quad:damping:fwbw: Exercise 1.21: Use the forward-backward scheme with quadratic damping --------------------------------------------------------------------- We consider the generalized model with quadratic damping, expressed as a system of two first-order equations as in the section :ref:`vib:ode2:staggered`: .. math:: \begin{align*} u^{\prime} &= v,\\ v' &= \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right){\thinspace .} \end{align*} However, contrary to what is done in the section :ref:`vib:ode2:staggered`, we want to apply the idea of a forward-backward discretization: :math:`u` is marched forward by a one-sided Forward Euler scheme applied to the first equation, and thereafter :math:`v` can be marched forward by a Backward Euler scheme in the second equation, see in the section :ref:`vib:model2x2:EulerCromer`. Express the idea in operator notation and write out the scheme. Unfortunately, the backward difference for the :math:`v` equation creates a nonlinearity :math:`|v^{n+1}|v^{n+1}`. To linearize this nonlinearity, use the known value :math:`v^n` inside the absolute value factor, i.e., :math:`|v^{n+1}|v^{n+1}\approx |v^n|v^{n+1}`. Show that the resulting scheme is equivalent to the one in the section :ref:`vib:ode2:staggered` for some time level :math:`n\geq 1`. What we learn from this exercise is that the first-order differences and the linearization trick play together in "the right way" such that the scheme is as good as when we (in the section :ref:`vib:ode2:staggered`) carefully apply centered differences and a geometric mean on a staggered mesh to achieve second-order accuracy. There is a difference in the handling of the initial conditions, though, as explained at the end of the section :ref:`vib:model2x2:EulerCromer`. Filename: ``vib_gen_bwdamping``. .. --- end exercise ---