.. !split .. 2DO: .. _undamped -> _simple everywhere .. Long time integration by adaptive RK: will that improve the .. phase error? Do experiments where we measure the wavelength .. and plot it as function of time. Can we vectorize the .. max/min pt computation? Vibration problems lead to differential equations with solutions that oscillate in time, typically in a damped or undamped sinusoidal fashion. Such solutions put certain demands on the numerical methods compared to other phenomena whose solutions are monotone or very smooth. Both the frequency and amplitude of the oscillations need to be accurately handled by the numerical schemes. Most of the reasoning and specific building blocks introduced in the forthcoming text can be reused to construct sound methods for partial differential equations of wave nature in multiple spatial dimensions. .. _vib:model1: Finite difference discretization ================================ Many of the numerical challenges faced when computing oscillatory solutions to ODEs and PDEs can be captured by the very simple ODE :math:`u^{\prime\prime} + u =0`. This ODE is thus chosen as our starting point for method development, implementation, and analysis. A basic model for vibrations ---------------------------- .. index:: vibration ODE .. index:: oscillations .. index:: mechanical vibrations A system that vibrates without damping and external forcing can be described by the ODE problem .. _Eq:vib:ode1: .. math:: \tag{1} u^{\prime\prime} + \omega^2u = 0,\quad u(0)=I,\ u^{\prime}(0)=0,\ t\in (0,T] {\thinspace .} Here, :math:`\omega` and :math:`I` are given constants. The exact solution of :ref:`(1) ` is .. index:: period (of oscillations) .. index:: frequency (of oscillations) .. index:: Hz (unit) .. _Eq:vib:ode1:uex: .. math:: \tag{2} u(t) = I\cos (\omega t) {\thinspace .} That is, :math:`u` oscillates with constant amplitude :math:`I` and angular frequency :math:`\omega`. The corresponding period of oscillations (i.e., the time between two neighboring peaks in the cosine function) is :math:`P=2\pi/\omega`. The number of periods per second is :math:`f=\omega/(2\pi)` and measured in the unit Hz. Both :math:`f` and :math:`\omega` are referred to as frequency, but :math:`\omega` is more precisely named *angular frequency*, measured in rad/s. In vibrating mechanical systems modeled by :ref:`(1) `, :math:`u(t)` very often represents a position or a displacement of a particular point in the system. The derivative :math:`u^{\prime}(t)` then has the interpretation of velocity, and :math:`u^{\prime\prime}(t)` is the associated acceleration. The model :ref:`(1) ` is not only applicable to vibrating mechanical systems, but also to oscillations in electrical circuits. .. _vib:ode1:fdm: A centered finite difference scheme ----------------------------------- To formulate a finite difference method for the model problem :ref:`(1) ` we follow the `four steps `__ explained in [Ref1]_. .. index:: single: mesh; finite differences .. index:: mesh function Step 1: Discretizing the domain ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The domain is discretized by introducing a uniformly partitioned time mesh. The points in the mesh are :math:`t_n=n\Delta t`, :math:`n=0,1,\ldots,N_t`, where :math:`\Delta t = T/N_t` is the constant length of the time steps. We introduce a mesh function :math:`u^n` for :math:`n=0,1,\ldots,N_t`, which approximates the exact solution at the mesh points. The mesh function will be computed from algebraic equations derived from the differential equation problem. Step 2: Fulfilling the equation at discrete time points ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The ODE is to be satisfied at each mesh point: .. _Eq:vib:ode1:step2: .. math:: \tag{3} u^{\prime\prime}(t_n) + \omega^2u(t_n) = 0,\quad n=1,\ldots,N_t {\thinspace .} .. index:: centered difference .. index:: single: finite differences; centered Step 3: Replacing derivatives by finite differences ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The derivative :math:`u^{\prime\prime}(t_n)` is to be replaced by a finite difference approximation. A common second-order accurate approximation to the second-order derivative is .. _Eq:vib:ode1:step3: .. math:: \tag{4} u^{\prime\prime}(t_n) \approx \frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} {\thinspace .} Inserting :ref:`(4) ` in :ref:`(3) ` yields .. _Eq:vib:ode1:step3b: .. math:: \tag{5} \frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} = -\omega^2 u^n {\thinspace .} We also need to replace the derivative in the initial condition by a finite difference. Here we choose a centered difference, whose accuracy is similar to the centered difference we used for :math:`u^{\prime\prime}`: .. _Eq:vib:ode1:step3c: .. math:: \tag{6} \frac{u^1-u^{-1}}{2\Delta t} = 0 {\thinspace .} Step 4: Formulating a recursive algorithm ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ To formulate the computational algorithm, we assume that we have already computed :math:`u^{n-1}` and :math:`u^n` such that :math:`u^{n+1}` is the unknown value, which we can readily solve for: .. _Eq:vib:ode1:step4: .. math:: \tag{7} u^{n+1} = 2u^n - u^{n-1} - \Delta t^2\omega^2 u^n {\thinspace .} The computational algorithm is simply to apply :ref:`(7) ` successively for :math:`n=1,2,\ldots,N_t-1`. This numerical scheme sometimes goes under the name Stormer's method or `Verlet integration `__. Computing the first step ~~~~~~~~~~~~~~~~~~~~~~~~ We observe that :ref:`(7) ` cannot be used for :math:`n=0` since the computation of :math:`u^1` then involves the undefined value :math:`u^{-1}` at :math:`t=-\Delta t`. The discretization of the initial condition then comes to our rescue: :ref:`(6) ` implies :math:`u^{-1} = u^1` and this relation can be combined with :ref:`(7) ` for :math:`n=1` to yield a value for :math:`u^1`: .. math:: u^1 = 2u^0 - u^{1} - \Delta t^2 \omega^2 u^0, which reduces to .. _Eq:vib:ode1:step4b: .. math:: \tag{8} u^1 = u^0 - \frac{1}{2} \Delta t^2 \omega^2 u^0 {\thinspace .} :ref:`vib:exer:step4b:alt` asks you to perform an alternative derivation and also to generalize the initial condition to :math:`u^{\prime}(0)=V\neq 0`. The computational algorithm ~~~~~~~~~~~~~~~~~~~~~~~~~~~ The steps for solving :ref:`(1) ` becomes 1. :math:`u^0=I` 2. compute :math:`u^1` from :ref:`(8) ` 3. for :math:`n=1,2,\ldots,N_t-1`: 1. compute :math:`u^{n+1}` from :ref:`(7) ` The algorithm is more precisely expressed directly in Python: .. code-block:: python t = linspace(0, T, Nt+1) # mesh points in time dt = t[1] - t[0] # constant time step u = zeros(Nt+1) # solution u[0] = I u[1] = u[0] - 0.5*dt**2*w**2*u[0] for n in range(1, Nt): u[n+1] = 2*u[n] - u[n-1] - dt**2*w**2*u[n] .. admonition:: Remark on using ``w`` for :math:`\omega` In the code, we use ``w`` as the symbol for :math:`\omega`. The reason is that this author prefers ``w`` for readability and comparison with the mathematical :math:`\omega` instead of the full word ``omega`` as variable name. Operator notation ~~~~~~~~~~~~~~~~~ We may write the scheme using a compact difference notation (see also `examples `__ in [Ref1]_). The difference :ref:`(4) ` has the operator notation :math:`[D_tD_t u]^n` such that we can write: .. _Eq:vib:ode1:step4:op: .. math:: \tag{9} [D_tD_t u + \omega^2 u = 0]^n {\thinspace .} Note that :math:`[D_tD_t u]^n` means applying a central difference with step :math:`\Delta t/2` twice: .. math:: [D_t(D_t u)]^n = \frac{[D_t u]^{n+\frac{1}{2}} - [D_t u]^{n-\frac{1}{2}}}{\Delta t} which is written out as .. math:: \frac{1}{\Delta t}\left(\frac{u^{n+1}-u^n}{\Delta t} - \frac{u^{n}-u^{n-1}}{\Delta t}\right) = \frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} {\thinspace .} The discretization of initial conditions can in the operator notation be expressed as .. _Eq:_auto1: .. math:: \tag{10} [u = I]^0,\quad [D_{2t} u = 0]^0, where the operator :math:`[D_{2t} u]^n` is defined as .. _Eq:_auto2: .. math:: \tag{11} [D_{2t} u]^n = \frac{u^{n+1} - u^{n-1}}{2\Delta t} {\thinspace .} .. _vib:impl1: Implementation (1) =========================== .. _vib:impl1:solver: Making a solver function ------------------------ The algorithm from the previous section is readily translated to a complete Python function for computing and returning :math:`u^0,u^1,\ldots,u^{N_t}` and :math:`t_0,t_1,\ldots,t_{N_t}`, given the input :math:`I`, :math:`\omega`, :math:`\Delta t`, and :math:`T`: .. code-block:: python import numpy as np import matplotlib.pyplot as plt def solver(I, w, dt, T): """ Solve u'' + w**2*u = 0 for t in (0,T], u(0)=I and u'(0)=0, by a central finite difference method with time step dt. """ dt = float(dt) Nt = int(round(T/dt)) u = np.zeros(Nt+1) t = np.linspace(0, Nt*dt, Nt+1) u[0] = I u[1] = u[0] - 0.5*dt**2*w**2*u[0] for n in range(1, Nt): u[n+1] = 2*u[n] - u[n-1] - dt**2*w**2*u[n] return u, t We do a simple ``from module import *`` to make the code as close as possible to MATLAB, although good programming habits would prefix the ``numpy`` and ``matplotlib`` calls by (abbreviations of) the module name. A function for plotting the numerical and the exact solution is also convenient to have: .. code-block:: python def u_exact(t, I, w): return I*np.cos(w*t) def visualize(u, t, I, w): plt.plot(t, u, 'r--o') t_fine = np.linspace(0, t[-1], 1001) # very fine mesh for u_e u_e = u_exact(t_fine, I, w) plt.hold('on') plt.plot(t_fine, u_e, 'b-') plt.legend(['numerical', 'exact'], loc='upper left') plt.xlabel('t') plt.ylabel('u') dt = t[1] - t[0] plt.title('dt=%g' % dt) umin = 1.2*u.min(); umax = -umin plt.axis([t[0], t[-1], umin, umax]) plt.savefig('tmp1.png'); plt.savefig('tmp1.pdf') A corresponding main program calling these functions for a simulation of a given number of periods (``num_periods``) may take the form .. code-block:: python I = 1 w = 2*pi dt = 0.05 num_periods = 5 P = 2*pi/w # one period T = P*num_periods u, t = solver(I, w, dt, T) visualize(u, t, I, w, dt) Adjusting some of the input parameters via the command line can be handy. Here is a code segment using the ``ArgumentParser`` tool in the ``argparse`` module to define option value (``--option value``) pairs on the command line: .. code-block:: python import argparse parser = argparse.ArgumentParser() parser.add_argument('--I', type=float, default=1.0) parser.add_argument('--w', type=float, default=2*pi) parser.add_argument('--dt', type=float, default=0.05) parser.add_argument('--num_periods', type=int, default=5) a = parser.parse_args() I, w, dt, num_periods = a.I, a.w, a.dt, a.num_periods Such parsing of the command line is explained in more detailed in the "section on user interfaces": "..." in [Ref1]_. A typical execution goes like .. code-block:: text Terminal> python vib_undamped.py --num_periods 20 --dt 0.1 Computing :math:`u^{\prime}` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In mechanical vibration applications one is often interested in computing the velocity :math:`v(t)=u^{\prime}(t)` after :math:`u(t)` has been computed. This can be done by a central difference, .. _Eq:_auto3: .. math:: \tag{12} v(t_n)=u^{\prime}(t_n) \approx v^n = \frac{u^{n+1}-u^{n-1}}{2\Delta t} = [D_{2t}u]^n {\thinspace .} This formula applies for all inner mesh points, :math:`n=1,\ldots,N_t-1`. For :math:`n=0`, :math:`v(0)` is given by the initial condition on :math:`u^{\prime}(0)`, and for :math:`n=N_t` we can use a one-sided, backward difference: .. math:: v^n=[D_t^-u]^n = \frac{u^{n} - u^{n-1}}{\Delta t}{\thinspace .} Typical (scalar) code is .. code-block:: python v = np.zeros_like(u) # or v = np.zeros(len(u)) # Use central difference for internal points for i in range(1, len(u)-1): v[i] = (u[i+1] - u[i-1])/(2*dt) # Use initial condition for u'(0) when i=0 v[0] = 0 # Use backward difference at the final mesh point v[-1] = (u[-1] - u[-2])/dt We can get rid of the loop, which is slow for large :math:`N_t`, by vectorizing the central difference. The above code segment goes as follows in its vectorized version: .. code-block:: python v = np.zeros_like(u) v[1:-1] = (u[2:] - u[:-2])/(2*dt) # central difference v[0] = 0 # boundary condition u'(0) v[-1] = (u[-1] - u[-2])/dt # backward difference .. _vib:ode1:verify: Verification (1) ------------------------- Manual calculation ~~~~~~~~~~~~~~~~~~ The simplest type of verification, which is also instructive for understanding the algorithm, is to compute :math:`u^1`, :math:`u^2`, and :math:`u^3` with the aid of a calculator and make a function for comparing these results with those from the ``solver`` function. The ``test_three_steps`` function in the file `vib_undamped.py `__ shows the details how we use the hand calculations to test the code: .. code-block:: python def test_three_steps(): from math import pi I = 1; w = 2*pi; dt = 0.1; T = 1 u_by_hand = np.array([1.000000000000000, 0.802607911978213, 0.288358920740053]) u, t = solver(I, w, dt, T) diff = np.abs(u_by_hand - u[:3]).max() tol = 1E-14 assert diff < tol Testing very simple solutions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Constructing test problems where the exact solution is constant or linear helps initial debugging and verification as one expects any reasonable numerical method to reproduce such solutions to machine precision. Second-order accurate methods will often also reproduce a quadratic solution. Here :math:`[D_tD_tt^2]^n=2`, which is the exact result. A solution :math:`u=t^2` leads to :math:`u^{\prime\prime}+\omega^2 u=2 + (\omega t)^2\neq 0`. We must therefore add a source in the equation: :math:`u^{\prime\prime} + \omega^2 u = f` to allow a solution :math:`u=t^2` for :math:`f=(\omega t)^2`. By simple insertion we can show that the mesh function :math:`u^n = t_n^2` is also a solution of the discrete equations. :ref:`vib:exer:undamped:verify:linquad` asks you to carry out all details to show that linear and quadratic solutions are solutions of the discrete equations. Such results are very useful for debugging and verification. You are strongly encouraged to do this problem now! Checking convergence rates ~~~~~~~~~~~~~~~~~~~~~~~~~~ Empirical computation of convergence rates yields a good method for verification. The method and its computational are explained in detail for a simple ODE model in the section on `computing convergence rates `__ in [Ref1]_. Readers not familiar with the concept should look up this reference before proceeding. In the present problem, computing convergence rates means that we must * perform :math:`m` simulations with halved time steps: :math:`\Delta t_i=2^{-i}\Delta t_0`, :math:`i=0,\ldots,m-1`, * compute the :math:`L^2` norm of the error, :math:`E_i=\sqrt{\Delta t_i\sum_{n=0}^{N_t-1}(u^n-{u_{\small\mbox{e}}}(t_n))^2}` in each case, * estimate the convergence rates :math:`r_i` based on two consecutive experiments :math:`(\Delta t_{i-1}, E_{i-1})` and :math:`(\Delta t_{i}, E_{i})`, assuming :math:`E_i=C(\Delta t_i)^{r}` and :math:`E_{i-1}=C(\Delta t_{i-1})^{r}`. From these equations it follows that :math:`r = \ln (E_{i-1}/E_i)/\ln (\Delta t_{i-1}/\Delta t_i)`. Since this :math:`r` will vary with :math:`i`, we equip it with an index and call it :math:`r_{i-1}`, where :math:`i` runs from :math:`1` to :math:`m-1`. The computed rates :math:`r_0,r_1,\ldots,r_{m-2}` hopefully converges to a number, which hopefully is 2, the right one, in the present problem. The convergence of the rates demands that the time steps :math:`\Delta t_i` are sufficiently small for the error model :math:`E_i=(\Delta t_i)^r` to be valid. All the implementational details of computing the sequence :math:`r_0,r_1,\ldots,r_{m-2}` appear below. .. code-block:: python def convergence_rates(m, solver_function, num_periods=8): """ Return m-1 empirical estimates of the convergence rate based on m simulations, where the time step is halved for each simulation. solver_function(I, w, dt, T) solves each problem, where T is based on simulation for num_periods periods. """ from math import pi w = 0.35; I = 0.3 # just chosen values P = 2*pi/w # period dt = P/30 # 30 time step per period 2*pi/w T = P*num_periods dt_values = [] E_values = [] for i in range(m): u, t = solver_function(I, w, dt, T) u_e = u_exact(t, I, w) E = np.sqrt(dt*np.sum((u_e-u)**2)) dt_values.append(dt) E_values.append(E) dt = dt/2 r = [np.log(E_values[i-1]/E_values[i])/ np.log(dt_values[i-1]/dt_values[i]) for i in range(1, m, 1)] return r The expected convergence rate is 2, because we have used a second-order finite difference approximations :math:`[D_tD_tu]^n` to the ODE and a second-order finite difference formula for the initial condition for :math:`u^{\prime}`. Other theoretical error measures also points to :math:`r=2`. In the present problem, when :math:`\Delta t_0` corresponds to 30 time steps per period, the returned ``r`` list has all its values equal to 2.00 (if rounded to two decimals). This amazing result means that all :math:`\Delta t_i` values are well into the asymptotic regime where the error model :math:`E_i = C(\Delta t_i)^r` is valid. We can now construct a test function that computes convergence rates and checks that the final (and usually the best) estimate is sufficiently close to 2. Here, a rough tolerance of 0.1 is enough. This unit test goes like .. code-block:: python def test_convergence_rates(): r = convergence_rates(m=5, solver_function=solver, num_periods=8) # Accept rate to 1 decimal place tol = 0.1 assert abs(r[-1] - 2.0) < tol The complete code appears in the file ``vib_undamped.py``. Scaled model ------------ It is advantageous to use dimensionless variables in simulations, because fewer parameters need to be set. The present problem is made dimensionless by introducing dimensionless variables :math:`\bar t = t/t_c` and :math:`\bar u = u/u_c`, where :math:`t_c` and :math:`u_c` are characteristic scales for :math:`t` and :math:`u`, respectively. The scaled ODE problem reads .. math:: \frac{u_c}{t_c^2}\frac{d^2\bar u}{d\bar t^2} + u_c\bar u = 0,\quad u_c\bar u(0) = I,\ \frac{u_c}{t_c}\frac{d\bar u}{d\bar t}(0)=0{\thinspace .} A common choice is to take :math:`t_c` as one period of the oscillations, :math:`t_c = 2\pi/w`, and :math:`u_c=I`. This gives the dimensionless model .. _Eq:vib:ode1:model:scaled: .. math:: \tag{13} \frac{d^2\bar u}{\bar t^2} + 4\pi^2 \bar u = 0,\quad \bar u(0)=1,\ \bar u^{\prime}(0)=0{\thinspace .} Observe that there are no physical parameters in :ref:`(13) `! We can therefore perform a single numerical simulation :math:`\bar u(\bar t)` and afterwards recover any :math:`u(t; \omega, I)` by .. math:: u(t;\omega, I) = u_c\bar u(t/t_c) = I\bar u(omega t/(2\pi)){\thinspace .} We can easily check this assertion: the solution of the scaled problem is :math:`\bar u(\bar t) = \cos(2\pi\bar t)`. The formula for :math:`u` in terms of :math:`\bar u` gives :math:`u = I\cos(\omega t)`, which is nothing but the solution of the original problem with dimensions. The scaled model can by run by calling ``solver(I=1, w=2*pi, dt, T)``. Each period is now 1 and ``T`` simply counts the number of periods. Choosing ``dt`` as ``1./M`` gives ``M`` time steps per period. .. _vib:ode1:longseries: Long time simulations ===================== Figure :ref:`vib:ode1:2dt` shows a comparison of the exact and numerical solution for the scaled model :ref:`(13) ` with :math:`\Delta t=0.1, 0.05`. From the plot we make the following observations: * The numerical solution seems to have correct amplitude. * There is a angular frequency error which is reduced by reducing the time step. * The total angular frequency error grows with time. By angular frequency error we mean that the numerical angular frequency differs from the exact :math:`\omega`. This is evident by looking at the peaks of the numerical solution: these have incorrect positions compared with the peaks of the exact cosine solution. The effect can be mathematical expressed by writing the numerical solution as :math:`I\cos\tilde\omega t`, where :math:`\tilde\omega` is not exactly equal to :math:`\omega`. Later, we shall mathematically quantify this numerical angular frequency :math:`\tilde\omega`. .. _vib:ode1:2dt: .. figure:: vib_freq_err1.png :width: 800 *Effect of halving the time step* Using a moving plot window -------------------------- In vibration problems it is often of interest to investigate the system's behavior over long time intervals. Errors in the angular frequency accumulate and become more visible as time grows. We can investigate long time series by introducing a moving plot window that can move along with the :math:`p` most recently computed periods of the solution. The `SciTools `__ package contains a convenient tool for this: ``MovingPlotWindow``. Typing ``pydoc scitools.MovingPlotWindow`` shows a demo and a description of its use. The function below utilizes the moving plot window and is in fact called by the ``main`` function the ``vib_undamped`` module if the number of periods in the simulation exceeds 10. .. code-block:: python def visualize_front(u, t, I, w, savefig=False, skip_frames=1): """ Visualize u and the exact solution vs t, using a moving plot window and continuous drawing of the curves as they evolve in time. Makes it easy to plot very long time series. Plots are saved to files if savefig is True. Only each skip_frames-th plot is saved (e.g., if skip_frame=10, only each 10th plot is saved to file; this is convenient if plot files corresponding to different time steps are to be compared). """ import scitools.std as st from scitools.MovingPlotWindow import MovingPlotWindow from math import pi # Remove all old plot files tmp_*.png import glob, os for filename in glob.glob('tmp_*.png'): os.remove(filename) P = 2*pi/w # one period umin = 1.2*u.min(); umax = -umin dt = t[1] - t[0] plot_manager = MovingPlotWindow( window_width=8*P, dt=dt, yaxis=[umin, umax], mode='continuous drawing') frame_counter = 0 for n in range(1,len(u)): if plot_manager.plot(n): s = plot_manager.first_index_in_plot st.plot(t[s:n+1], u[s:n+1], 'r-1', t[s:n+1], I*cos(w*t)[s:n+1], 'b-1', title='t=%6.3f' % t[n], axis=plot_manager.axis(), show=not savefig) # drop window if savefig if savefig and n % skip_frames == 0: filename = 'tmp_%04d.png' % frame_counter st.savefig(filename) print 'making plot file', filename, 'at t=%g' % t[n] frame_counter += 1 plot_manager.update(n) We run the scaled problem (the default values for the command-line arguments ``--I`` and ``--w`` correspond to the scaled problem) for 40 periods with 20 time steps per period: .. code-block:: text Terminal> python vib_undamped.py --dt 0.05 --num_periods 40 The moving plot window is invoked, and we can follow the numerical and exact solutions as time progresses. From this demo we see that the angular frequency error is small in the beginning, but it becomes more prominent with time. A new run with :math:`\Delta t=0.1` (i.e., only 10 time steps per period) clearly shows that the phase errors become significant even earlier in the time series, deteriorating the solution further. .. _vib:ode1:anim: Making animations ----------------- .. index:: making movies .. index:: animation .. index:: WebM (video format) .. index:: Ogg (video format) .. index:: MP4 (video format) .. index:: Flash (video format) .. index:: video formats Producing standard video formats ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The ``visualize_front`` function stores all the plots in files whose names are numbered: ``tmp_0000.png``, ``tmp_0001.png``, ``tmp_0002.png``, and so on. From these files we may make a movie. The Flash format is popular, .. code-block:: text Terminal> ffmpeg -r 12 -i tmp_%04d.png -c:v flv movie.flv The ``ffmpeg`` program can be replaced by the ``avconv`` program in the above command if desired (but at the time of this writing it seems to be more momentum in the ``ffmpeg`` project). The ``-r`` option should come first and describes the number of frames per second in the movie. The ``-i`` option describes the name of the plot files. Other formats can be generated by changing the video codec and equipping the video file with the right extension: ====== ============================ Format Codec and filename ====== ============================ Flash ``-c:v flv movie.flv`` MP4 ``-c:v libx264 movie.mp4`` WebM ``-c:v libvpx movie.webm`` Ogg ``-c:v libtheora movie.ogg`` ====== ============================ .. index:: HTML5 video tag The video file can be played by some video player like ``vlc``, ``mplayer``, ``gxine``, or ``totem``, e.g., .. code-block:: text Terminal> vlc movie.webm A web page can also be used to play the movie. Today's standard is to use the HTML5 ``video`` tag: .. code-block:: html Modern browsers do not support all of the video formats. MP4 is needed to successfully play the videos on Apple devices that use the Safari browser. WebM is the preferred format for Chrome, Opera, Firefox, and Internet Explorer v9+. Flash was a popular format, but older browsers that required Flash can play MP4. All browsers that work with Ogg can also work with WebM. This means that to have a video work in all browsers, the video should be available in the MP4 and WebM formats. The proper HTML code reads .. code-block:: html The MP4 format should appear first to ensure that Apple devices will load the video correctly. .. admonition:: Caution: number the plot files correctly To ensure that the individual plot frames are shown in correct order, it is important to number the files with zero-padded numbers (0000, 0001, 0002, etc.). The printf format ``%04d`` specifies an integer in a field of width 4, padded with zeros from the left. A simple Unix wildcard file specification like ``tmp_*.png`` will then list the frames in the right order. If the numbers in the filenames were not zero-padded, the frame ``tmp_11.png`` would appear before ``tmp_2.png`` in the movie. Paying PNG files in a web browser ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. index:: scitools movie command The ``scitools movie`` command can create a movie player for a set of PNG files such that a web browser can be used to watch the movie. This interface has the advantage that the speed of the movie can easily be controlled, a feature that scientists often appreciate. The command for creating an HTML with a player for a set of PNG files ``tmp_*.png`` goes like .. code-block:: text Terminal> scitools movie output_file=vib.html fps=4 tmp_*.png The ``fps`` argument controls the speed of the movie ("frames per second"). To watch the movie, load the video file ``vib.html`` into some browser, e.g., .. code-block:: text Terminal> google-chrome vib.html # invoke web page Clicking on ``Start movie`` to see the result. Moving this movie to some other place requires moving ``vib.html`` *and all the PNG files* ``tmp_*.png``: .. code-block:: text Terminal> mkdir vib_dt0.1 Terminal> mv tmp_*.png vib_dt0.1 Terminal> mv vib.html vib_dt0.1/index.html Making animated GIF files ~~~~~~~~~~~~~~~~~~~~~~~~~ The ``convert`` program from the ImageMagick software suite can be used to produce animated GIF files from a set of PNG files: .. code-block:: text Terminal> convert -delay 25 tmp_vib*.png tmp_vib.gif The ``-delay`` option needs an argument of the delay between each frame, measured in 1/100 s, so 4 frames/s here gives 25/100 s delay. Note, however, that in this particular example with :math:`\Delta t=0.05` and 40 periods, making an animated GIF file out of the large number of PNG files is a very heavy process and not considered feasible. Animated GIFs are best suited for animations with not so many frames and where you want to see each frame and play them slowly. Using Bokeh to compare graphs ----------------------------- Instead of a moving plot frame, one can use tools that allows panning by the mouse. For example, we can show four periods of a signal in a plot and then scroll with the mouse through the rest of the simulation. The `Bokeh `__ plotting library offers such tools, but the plot must be displayed in a web browser. The documentation of Bokeh is excellent, so here we just show how the library can be used to compare a set of :math:`u` curves corresponding to long time simulations. Imagine we have performed experiments for a set of :math:`\Delta t` values. We want each curve, together with the exact solution, to appear in a plot, and then arrange all plots in a grid-like fashion: .. figure:: bokeh_gridplot1.png :width: 800 Furthermore, we want the axis to couple such that if we move into the future in one plot, all the other plots follows (note the displaced :math:`t` axes!): .. figure:: bokeh_gridplot2.png :width: 800 A function for creating a Bokeh plot, given a list of ``u`` arrays and corresponding ``t`` arrays, from different simulations, described compactly in a list of strings ``legends``, takes the following form: .. code-block:: python def bokeh_plot(u, t, legends, I, w, t_range, filename): """ Make plots for u vs t using the Bokeh library. u and t are lists (several experiments can be compared). legens contain legend strings for the various u,t pairs. """ if not isinstance(u, (list,tuple)): u = [u] # wrap in list if not isinstance(t, (list,tuple)): t = [t] # wrap in list if not isinstance(legends, (list,tuple)): legends = [legends] # wrap in list import bokeh.plotting as plt plt.output_file(filename, mode='cdn', title='Comparison') # Assume that all t arrays have the same range t_fine = np.linspace(0, t[0][-1], 1001) # fine mesh for u_e tools = 'pan,wheel_zoom,box_zoom,reset,'\ 'save,box_select,lasso_select' u_range = [-1.2*I, 1.2*I] font_size = '8pt' p = [] # list of plot objects # Make the first figure p_ = plt.figure( width=300, plot_height=250, title=legends[0], x_axis_label='t', y_axis_label='u', x_range=t_range, y_range=u_range, tools=tools, title_text_font_size=font_size) p_.xaxis.axis_label_text_font_size=font_size p_.yaxis.axis_label_text_font_size=font_size p_.line(t[0], u[0], line_color='blue') # Add exact solution u_e = u_exact(t_fine, I, w) p_.line(t_fine, u_e, line_color='red', line_dash='4 4') p.append(p_) # Make the rest of the figures and attach their axes to # the first figure's axes for i in range(1, len(t)): p_ = plt.figure( width=300, plot_height=250, title=legends[i], x_axis_label='t', y_axis_label='u', x_range=p[0].x_range, y_range=p[0].y_range, tools=tools, title_text_font_size=font_size) p_.xaxis.axis_label_text_font_size = font_size p_.yaxis.axis_label_text_font_size = font_size p_.line(t[i], u[i], line_color='blue') p_.line(t_fine, u_e, line_color='red', line_dash='4 4') p.append(p_) # Arrange all plots in a grid with 3 plots per row grid = [[]] for i, p_ in enumerate(p): grid[-1].append(p_) if (i+1) % 3 == 0: # New row grid.append([]) plot = plt.gridplot(grid, toolbar_location='left') plt.save(plot) plt.show(plot) A particular example using the ``bokeh_plot`` function appears below. .. code-block:: python def demo_bokeh(): """Solve a scaled ODE u'' + u = 0.""" from math import pi w = 1.0 # Scaled problem (frequency) P = 2*np.pi/w # Period num_steps_per_period = [5, 10, 20, 40, 80] T = 40*P # Simulation time: 40 periods u = [] # List of numerical solutions t = [] # List of corresponding meshes legends = [] for n in num_steps_per_period: dt = P/n u_, t_ = solver(I=1, w=w, dt=dt, T=T) u.append(u_) t.append(t_) legends.append('# time steps per period: %d' % n) bokeh_plot(u, t, legends, I=1, w=w, t_range=[0, 4*P], filename='tmp.html') Using a line-by-line ascii plotter ---------------------------------- Plotting functions vertically, line by line, in the terminal window using ascii characters only is a simple, fast, and convenient visualization technique for long time series. Note that the time axis then is positive downwards on the screen. The tool ``scitools.avplotter.Plotter`` makes it easy to create such plots: .. code-block:: python def visualize_front_ascii(u, t, I, w, fps=10): """ Plot u and the exact solution vs t line by line in a terminal window (only using ascii characters). Makes it easy to plot very long time series. """ from scitools.avplotter import Plotter import time from math import pi P = 2*pi/w umin = 1.2*u.min(); umax = -umin p = Plotter(ymin=umin, ymax=umax, width=60, symbols='+o') for n in range(len(u)): print p.plot(t[n], u[n], I*cos(w*t[n])), \ '%.1f' % (t[n]/P) time.sleep(1/float(fps)) The call ``p.plot`` returns a line of text, with the :math:`t` axis marked and a symbol ``+`` for the first function (``u``) and ``o`` for the second function (the exact solution). Here we append to this text a time counter reflecting how many periods the current time point corresponds to. A typical output (:math:`\omega =2\pi`, :math:`\Delta t=0.05`) looks like this: .. code-block:: text | o+ 14.0 | + o 14.0 | + o 14.1 | + o 14.1 | + o 14.2 +| o 14.2 + | 14.2 + o | 14.3 + o | 14.4 + o | 14.4 +o | 14.5 o + | 14.5 o + | 14.6 o + | 14.6 o + | 14.7 o | + 14.7 | + 14.8 | o + 14.8 | o + 14.9 | o + 14.9 | o+ 15.0 .. _vib:ode1:empirical: Empirical analysis of the solution ---------------------------------- For oscillating functions like those in Figure :ref:`vib:ode1:2dt` we may compute the amplitude and frequency (or period) empirically. That is, we run through the discrete solution points :math:`(t_n, u_n)` and find all maxima and minima points. The distance between two consecutive maxima (or minima) points can be used as estimate of the local period, while half the difference between the :math:`u` value at a maximum and a nearby minimum gives an estimate of the local amplitude. The local maxima are the points where .. _Eq:_auto4: .. math:: \tag{14} u^{n-1} < u^n > u^{n+1},\quad n=1,\ldots,N_t-1, and the local minima are recognized by .. _Eq:_auto5: .. math:: \tag{15} u^{n-1} > u^n < u^{n+1},\quad n=1,\ldots,N_t-1 {\thinspace .} In computer code this becomes .. code-block:: python def minmax(t, u): minima = []; maxima = [] for n in range(1, len(u)-1, 1): if u[n-1] > u[n] < u[n+1]: minima.append((t[n], u[n])) if u[n-1] < u[n] > u[n+1]: maxima.append((t[n], u[n])) return minima, maxima Note that the two returned objects are lists of tuples. Let :math:`(t_i, e_i)`, :math:`i=0,\ldots,M-1`, be the sequence of all the :math:`M` maxima points, where :math:`t_i` is the time value and :math:`e_i` the corresponding :math:`u` value. The local period can be defined as :math:`p_i=t_{i+1}-t_i`. With Python syntax this reads .. code-block:: python def periods(maxima): p = [extrema[n][0] - maxima[n-1][0] for n in range(1, len(maxima))] return np.array(p) The list ``p`` created by a list comprehension is converted to an array since we probably want to compute with it, e.g., find the corresponding frequencies ``2*pi/p``. Having the minima and the maxima, the local amplitude can be calculated as the difference between two neighboring minimum and maximum points: .. code-block:: python def amplitudes(minima, maxima): a = [(abs(maxima[n][1] - minima[n][1]))/2.0 for n in range(min(len(minima),len(maxima)))] return np.array(a) The code segments are found in the file `vib_empirical_analysis.py `__. Since ``a[i]`` and ``p[i]`` correspond to the :math:`i`-th amplitude estimate and the :math:`i`-th period estimate, respectively, it is most convenient to visualize the ``a`` and ``p`` values with the index ``i`` on the horizontal axis. (There is no unique time point associated with either of these estimate since values at two different time points were used in the computations.) In the analysis of very long time series, it is advantageous to compute and plot ``p`` and ``a`` instead of :math:`u` to get an impression of the development of the oscillations. Let us do this for the scaled problem and :math:`\Delta t=0.1, 0.05, 0.01`. A ready-made function .. code-block:: python plot_empirical_freq_and_amplitude(u, t, I, w) computes the empirical amplitudes and periods, and creates a plot where the amplitudes and angular frequencies are visualized together with the exact amplitude ``I`` and the exact angular frequency ``w``. We can make a little program for creating the plot: .. code-block:: python from vib_undamped import solver, plot_empirical_freq_and_amplitude from math import pi dt_values = [0.1, 0.05, 0.01] u_cases = [] t_cases = [] for dt in dt_values: # Simulate scaled problem for 40 periods u, t = solver(I=1, w=2*pi, dt=dt, T=40) u_cases.append(u) t_cases.append(t) plot_empirical_freq_and_amplitude(u_cases, t_cases, I=1, w=2*pi) Figure :ref:`vib:ode1:fig:freq_ampl` shows the result: we clearly see that lowering :math:`\Delta t` improves the angular frequency significantly, while the amplitude seems to be more accurate. The lines with :math:`\Delta t=0.01`, corresponding to 100 steps per period, can hardly be distinguished from the exact values. The next section shows how we can get mathematical insight into why amplitudes are good and frequencies are more inaccurate. .. _vib:ode1:fig:freq_ampl: .. figure:: empirical_ampl_freq.png :width: 800 *Empirical amplitude and angular frequency for three cases of time steps* .. Use it for very long time integration of CN! And of RK4! .. _vib:ode1:analysis: Analysis of the numerical scheme ================================ Deriving a solution of the numerical scheme ------------------------------------------- After having seen the phase error grow with time in the previous section, we shall now quantify this error through mathematical analysis. The key tool in the analysis will be to establish an exact solution of the discrete equations. The difference equation :ref:`(7) ` has constant coefficients and is homogeneous. Such equations are known to have solutions on the form :math:`u^n=CA^n`, where :math:`A` is some number to be determined from the difference equation and :math:`C` is found as the initial condition (:math:`C=I`). Recall that :math:`n` in :math:`u^n` is a superscript labeling the time level, while :math:`n` in :math:`A^n` is an exponent. With oscillating functions as solutions, the algebra will be considerably simplified if we seek an :math:`A` on the form .. math:: A=e^{i\tilde\omega \Delta t}, and solve for the numerical frequency :math:`\tilde\omega` rather than :math:`A`. Note that :math:`i=\sqrt{-1}` is the imaginary unit. (Using a complex exponential function gives simpler arithmetics than working with a sine or cosine function.) We have .. math:: A^n = e^{i\tilde\omega \Delta t\, n}=e^{i\tilde\omega t} = \cos (\tilde\omega t) + i\sin(\tilde \omega t) {\thinspace .} The physically relevant numerical solution can be taken as the real part of this complex expression. The calculations go as .. math:: [D_tD_t u]^n &= \frac{u^{n+1} - 2u^n + u^{n-1}}{\Delta t^2}\\ &= I\frac{A^{n+1} - 2A^n + A^{n-1}}{\Delta t^2}\\ &= \frac{I}{\Delta t^{2}}(e^{i\tilde\omega(t+\Delta t)} - 2e^{i\tilde\omega t} + e^{i\tilde\omega(t-\Delta t)})\\ &= Ie^{i\tilde\omega t}\frac{1}{\Delta t^2}\left(e^{i\tilde\omega\Delta t} + e^{i\tilde\omega(-\Delta t)} - 2\right)\\ &= Ie^{i\tilde\omega t}\frac{2}{\Delta t^2}\left(\cosh(i\tilde\omega\Delta t) -1 \right)\\ &= Ie^{i\tilde\omega t}\frac{2}{\Delta t^2}\left(\cos(\tilde\omega\Delta t) -1 \right)\\ &= -Ie^{i\tilde\omega t}\frac{4}{\Delta t^2}\sin^2(\frac{\tilde\omega\Delta t}{2}) The last line follows from the relation :math:`\cos x - 1 = -2\sin^2(x/2)` (try ``cos(x)-1`` in `wolframalpha.com `__ to see the formula). The scheme :ref:`(7) ` with :math:`u^n=Ie^{i\omega\tilde\Delta t\, n}` inserted now gives .. _Eq:_auto6: .. math:: \tag{16} -Ie^{i\tilde\omega t}\frac{4}{\Delta t^2}\sin^2(\frac{\tilde\omega\Delta t}{2}) + \omega^2 Ie^{i\tilde\omega t} = 0, which after dividing by :math:`Ie^{i\tilde\omega t}` results in .. _Eq:_auto7: .. math:: \tag{17} \frac{4}{\Delta t^2}\sin^2(\frac{\tilde\omega\Delta t}{2}) = \omega^2 {\thinspace .} The first step in solving for the unknown :math:`\tilde\omega` is .. math:: \sin^2(\frac{\tilde\omega\Delta t}{2}) = \left(\frac{\omega\Delta t}{2}\right)^2 {\thinspace .} Then, taking the square root, applying the inverse sine function, and multiplying by :math:`2/\Delta t`, results in .. _Eq:vib:ode1:tildeomega: .. math:: \tag{18} \tilde\omega = \pm \frac{2}{\Delta t}\sin^{-1}\left(\frac{\omega\Delta t}{2}\right) {\thinspace .} The first observation of :ref:`(18) ` tells that there is a phase error since the numerical frequency :math:`\tilde\omega` never equals the exact frequency :math:`\omega`. But how good is the approximation :ref:`(18) `? That is, what is the error :math:`\omega - \tilde\omega` or :math:`\tilde\omega/\omega`? Taylor series expansion for small :math:`\Delta t` may give an expression that is easier to understand than the complicated function in :ref:`(18) `: .. code-block:: ipy >>> from sympy import * >>> dt, w = symbols('dt w') >>> w_tilde_e = 2/dt*asin(w*dt/2) >>> w_tilde_series = w_tilde_e.series(dt, 0, 4) >>> print w_tilde_series w + dt**2*w**3/24 + O(dt**4) This means that .. See vib_symbolic.py for computations with sympy .. _Eq:vib:ode1:tildeomega:series: .. math:: \tag{19} \tilde\omega = \omega\left( 1 + \frac{1}{24}\omega^2\Delta t^2\right) + {\mathcal{O}(\Delta t^4)} {\thinspace .} The error in the numerical frequency is of second-order in :math:`\Delta t`, and the error vanishes as :math:`\Delta t\rightarrow 0`. We see that :math:`\tilde\omega > \omega` since the term :math:`\omega^3\Delta t^2/24 >0` and this is by far the biggest term in the series expansion for small :math:`\omega\Delta t`. A numerical frequency that is too large gives an oscillating curve that oscillates too fast and therefore "lags behind" the exact oscillations, a feature that can be seen in the left plot in Figure :ref:`vib:ode1:2dt`. Figure :ref:`vib:ode1:tildeomega:plot` plots the discrete frequency :ref:`(18) ` and its approximation :ref:`(19) ` for :math:`\omega =1` (based on the program `vib_plot_freq.py `__). Although :math:`\tilde\omega` is a function of :math:`\Delta t` in :ref:`(19) `, it is misleading to think of :math:`\Delta t` as the important discretization parameter. It is the product :math:`\omega\Delta t` that is the key discretization parameter. This quantity reflects the *number of time steps per period* of the oscillations. To see this, we set :math:`P=N_P\Delta t`, where :math:`P` is the length of a period, and :math:`N_P` is the number of time steps during a period. Since :math:`P` and :math:`\omega` are related by :math:`P=2\pi/\omega`, we get that :math:`\omega\Delta t = 2\pi/N_P`, which shows that :math:`\omega\Delta t` is directly related to :math:`N_P`. The plot shows that at least :math:`N_P\sim 25-30` points per period are necessary for reasonable accuracy, but this depends on the length of the simulation (:math:`T`) as the total phase error due to the frequency error grows linearly with time (see :ref:`vib:exer:phase:err:growth`). .. _vib:ode1:tildeomega:plot: .. figure:: discrete_freq.png :width: 400 *Exact discrete frequency and its second-order series expansion* .. _vib:ode1:analysis:sol: Exact discrete solution ----------------------- Perhaps more important than the :math:`\tilde\omega = \omega + {\cal O}(\Delta t^2)` result found above is the fact that we have an exact discrete solution of the problem: .. _Eq:vib:ode1:un:exact: .. math:: \tag{20} u^n = I\cos\left(\tilde\omega n\Delta t\right),\quad \tilde\omega = \frac{2}{\Delta t}\sin^{-1}\left(\frac{\omega\Delta t}{2}\right) {\thinspace .} We can then compute the error mesh function .. _Eq:vib:ode1:en: .. math:: \tag{21} e^n = {u_{\small\mbox{e}}}(t_n) - u^n = I\cos\left(\omega n\Delta t\right) - I\cos\left(\tilde\omega n\Delta t\right){\thinspace .} From the formula :math:`\cos 2x - \cos 2y = -2\sin(x-y)\sin(x+y)` we can rewrite :math:`e^n` so the expression is easier to interpret: .. _Eq:vib:ode1:en2: .. math:: \tag{22} e^n = -2I\sin\left(t\frac{1}{2}\left( \omega - \tilde\omega\right)\right) \sin\left(t\frac{1}{2}\left( \omega + \tilde\omega\right)\right){\thinspace .} The error mesh function is ideal for verification purposes and you are strongly encouraged to make a test based on :ref:`(20) ` by doing :ref:`vib:exer:discrete:omega`. .. _vib:ode1:analysis:conv: Convergence ----------- We can use :ref:`(19) `, :ref:`(21) `, or :ref:`(22) ` to show *convergence* of the numerical scheme, i.e., :math:`e^n\rightarrow 0` as :math:`\Delta t\rightarrow 0`. We have that .. math:: \lim_{\Delta t\rightarrow 0} \tilde\omega = \lim_{\Delta t\rightarrow 0} \frac{2}{\Delta t}\sin^{-1}\left(\frac{\omega\Delta t}{2}\right) = \omega, by L'Hopital's rule or simply asking ``sympy`` or `WolframAlpha `__ about the limit: .. code-block:: python >>> import sympy as sym >>> dt, w = sym.symbols('x w') >>> sym.limit((2/dt)*sym.asin(w*dt/2), dt, 0, dir='+') w Also :ref:`(19) ` can be used to establish this result that :math:`\tilde\omega\rightarrow\omega`. It then follows from the expression(s) for :math:`e^n` that :math:`e^n\rightarrow 0`. The global error ---------------- .. index:: single: error; global To achieve more analytical insight into the nature of the global error, we can Taylor expand the error mesh function :ref:`(21) `. Since :math:`\tilde\omega` in :ref:`(18) ` contains :math:`\Delta t` in the denominator we use the series expansion for :math:`\tilde\omega` inside the cosine function. A relevant ``sympy`` session is .. code-block:: python >>> from sympy import * >>> dt, w, t = symbols('dt w t') >>> w_tilde_e = 2/dt*asin(w*dt/2) >>> w_tilde_series = w_tilde_e.series(dt, 0, 4) >>> w_tilde_series w + dt**2*w**3/24 + O(dt**4) Series expansions in ``sympy`` have the inconvenient ``O()`` term that prevents further calculations with the series. We can use the ``removeO()`` command to get rid of the ``O()`` term: .. code-block:: python >>> w_tilde_series = w_tilde_series.removeO() >>> w_tilde_series dt**2*w**3/24 + w Using this ``w_tilde_series`` expression for :math:`\tilde w` in :ref:`(21) `, dropping :math:`I` (which is a common factor), and performing a series expansion of the error yields .. code-block:: python >>> error = cos(w*t) - cos(w_tilde_series*t) >>> error.series(dt, 0, 6) dt**2*t*w**3*sin(t*w)/24 + dt**4*t**2*w**6*cos(t*w)/1152 + O(dt**6) Since we are mainly interested in the leading-order term in such expansions (the term with lowest power in :math:`\Delta t` and goes most slowly to zero), we use the ``.as_leading_term(dt)`` construction to pick out this term: .. code-block:: python >>> error.series(dt, 0, 6).as_leading_term(dt) dt**2*t*w**3*sin(t*w)/24 The last result means that the leading order global (true) error at a point :math:`t` is proportional to :math:`\omega^3t\Delta t^2`. Now, :math:`t` is related to :math:`\Delta t` through :math:`t=n\Delta t`. The factor :math:`\sin(\omega t)` can at most be 1, so we use this value to bound the leading-order expression to its maximum value .. math:: e^n = \frac{1}{24}n\omega^3\Delta t^3{\thinspace .} This is the dominating term of the error *at a point*. We are interested in the accumulated global error, which can be taken as the :math:`\ell^2` norm of :math:`e^n`. The norm is simply computed by summing contributions from all mesh points: .. math:: ||e^n||_{\ell^2}^2 = \Delta t\sum_{n=0}^{N_t} \frac{1}{24^2}n^2\omega^6\Delta t^6 =\frac{1}{24^2}\omega^6\Delta t^7 \sum_{n=0}^{N_t} n^2{\thinspace .} The sum :math:`\sum_{n=0}^{N_t} n^2` is approximately equal to :math:`\frac{1}{3}N_t^3`. Replacing :math:`N_t` by :math:`T/\Delta t` and taking the square root gives the expression .. math:: ||e^n||_{\ell^2} = \frac{1}{24}\sqrt{\frac{T^3}{3}}\omega^3\Delta t^2{\thinspace .} This is our expression for the global (or integrated) error. The main result from this expression is that also the global error is proportional to :math:`\Delta t^2`. Stability --------- Looking at :ref:`(20) `, it appears that the numerical solution has constant and correct amplitude, but an error in the angular frequency. A constant amplitude is not necessarily the case, however! To see this, note that if only :math:`\Delta t` is large enough, the magnitude of the argument to :math:`\sin^{-1}` in :ref:`(18) ` may be larger than 1, i.e., :math:`\omega\Delta t/2 > 1`. In this case, :math:`\sin^{-1}(\omega\Delta t/2)` has a complex value and therefore :math:`\tilde\omega` becomes complex. Type, for example, ``asin(x)`` in `wolframalpha.com `__ to see basic properties of :math:`\sin^{-1} (x)`). A complex :math:`\tilde\omega` can be written :math:`\tilde\omega = \tilde\omega_r + i\tilde\omega_i`. Since :math:`\sin^{-1}(x)` has a *negative* imaginary part for :math:`x>1`, :math:`\tilde\omega_i < 0`, which means that :math:`e^{i\tilde\omega t}=e^{-\tilde\omega_i t}e^{i\tilde\omega_r t}` will lead to exponential growth in time because :math:`e^{-\tilde\omega_i t}` with :math:`\tilde\omega_i <0` has a positive exponent. .. index:: stability criterion .. admonition:: Stability criterion We do not tolerate growth in the amplitude since such growth is not present in the exact solution. Therefore, we must impose a *stability criterion* that the argument in the inverse sine function leads to real and not complex values of :math:`\tilde\omega`. The stability criterion reads .. _Eq:_auto8: .. math:: \tag{23} \frac{\omega\Delta t}{2} \leq 1\quad\Rightarrow\quad \Delta t \leq \frac{2}{\omega} {\thinspace .} With :math:`\omega =2\pi`, :math:`\Delta t > \pi^{-1} = 0.3183098861837907` will give growing solutions. Figure :ref:`vib:ode1:dt:unstable` displays what happens when :math:`\Delta t =0.3184`, which is slightly above the critical value: :math:`\Delta t =\pi^{-1} + 9.01\cdot 10^{-5}`. .. _vib:ode1:dt:unstable: .. figure:: vib_unstable.png :width: 400 *Growing, unstable solution because of a time step slightly beyond the stability limit* About the accuracy at the stability limit ----------------------------------------- An interesting question is whether the stability condition :math:`\Delta t < 2/\omega` is unfortunate, or more precisely: would it be meaningful to take larger time steps to speed up computations? The answer is a clear no. At the stability limit, we have that :math:`\sin^{-1}\omega\Delta t/2 = \sin^{-1} 1 = \pi/2`, and therefore :math:`\tilde\omega = \pi/\Delta t`. (Note that the approximate formula :ref:`(19) ` is very inaccurate for this value of :math:`\Delta t` as it predicts :math:`\tilde\omega = 2.34/pi`, which is a 25 percent reduction.) The corresponding period of the numerical solution is :math:`\tilde P=2\pi/\tilde\omega = 2\Delta t`, which means that there is just one time step :math:`\Delta t` between a peak (maximum) and a `through `__ (minimum) in the numerical solution. This is the shortest possible wave that can be represented in the mesh! In other words, it is not meaningful to use a larger time step than the stability limit. Also, the error in angular frequency when :math:`\Delta t = 2/\omega` is severe: Figure :ref:`vib:ode1:dt:stablimit` shows a comparison of the numerical and analytical solution with :math:`\omega = 2\pi` and :math:`\Delta t = 2/\omega = \pi^{-1}`. Already after one period, the numerical solution has a through while the exact solution has a peak (!). The error in frequency when :math:`\Delta t` is at the stability limit becomes :math:`\omega - \tilde\omega = \omega(1-\pi/2)\approx -0.57\omega`. The corresponding error in the period is :math:`P - \tilde P \approx 0.36P`. The error after :math:`m` periods is then :math:`0.36mP`. This error has reached half a period when :math:`m=1/(2\cdot 0.36)\approx 1.38`, which theoretically confirms the observations in Figure :ref:`vib:ode1:dt:stablimit` that the numerical solution is a through ahead of a peak already after one and a half period. Consequently, :math:`\Delta t` should be chosen much less than the stability limit to achieve meaningful numerical computations. .. _vib:ode1:dt:stablimit: .. figure:: vib_stability_limit.png :width: 400 *Numerical solution with :math:`\Delta t` exactly at the stability limit* .. admonition:: Summary From the accuracy and stability analysis we can draw three important conclusions: 1. The key parameter in the formulas is :math:`p=\omega\Delta t`. The period of oscillations is :math:`P=2\pi/\omega`, and the number of time steps per period is :math:`N_P=P/\Delta t`. Therefore, :math:`p=\omega\Delta t = 2\pi N_P`, showing that the critical parameter is the number of time steps per period. The smallest possible :math:`N_P` is 2, showing that :math:`p\in (0,\pi]`. 2. Provided :math:`p\leq 2`, the amplitude of the numerical solution is constant. 3. The ratio of the numerical angular frequency and the exact one is :math:`\tilde\omega/\omega \approx 1 + \frac{1}{24}p^2`. The error :math:`\frac{1}{24}p^2` leads to wrongly displaced peaks of the numerical solution, and the error in peak location grows linearly with time (see :ref:`vib:exer:phase:err:growth`). .. _vib:model2x2: Alternative schemes based on 1st-order equations ================================================ A standard technique for solving second-order ODEs is to rewrite them as a system of first-order ODEs and then choose a solution strategy from the vast collection of methods for first-order ODE systems. Given the second-order ODE problem .. math:: u^{\prime\prime} + \omega^2 u = 0,\quad u(0)=I,\ u^{\prime}(0)=0, we introduce the auxiliary variable :math:`v=u^{\prime}` and express the ODE problem in terms of first-order derivatives of :math:`u` and :math:`v`: .. _Eq:vib:model2x2:ueq: .. math:: \tag{24} u^{\prime} = v, .. _Eq:vib:model2x2:veq: .. math:: \tag{25} v' = -\omega^2 u {\thinspace .} The initial conditions become :math:`u(0)=I` and :math:`v(0)=0`. The Forward Euler scheme ------------------------ A Forward Euler approximation to our :math:`2\times 2` system of ODEs :ref:`(24) `-:ref:`(25) ` becomes .. _Eq:_auto9: .. math:: \tag{26} \lbrack D_t^+ u = v\rbrack^n, \lbrack D_t^+ v = -\omega^2 u\rbrack^n, or written out, .. _Eq:vib:undamped:FE1: .. math:: \tag{27} u^{n+1} = u^n + \Delta t v^n, .. _Eq:vib:undamped:FE2: .. math:: \tag{28} v^{n+1} = v^n -\Delta t \omega^2 u^n {\thinspace .} Let us briefly compare this Forward Euler method with the centered difference scheme for the second-order differential equation. We have from :ref:`(27) ` and :ref:`(28) ` applied at levels :math:`n` and :math:`n-1` that .. math:: u^{n+1} = u^n + \Delta t v^n = u^n + \Delta t (v^{n-1} -\Delta t \omega^2 u^{n-1}{\thinspace .} Since from :ref:`(27) ` .. math:: v^{n-1} = \frac{1}{\Delta t}(u^{n}-u^{n-1}), it follows that .. math:: u^{n+1} = 2u^n - u^{n-1} -\Delta t^2\omega^2 u^{n-1}, which is very close to the centered difference scheme, but the last term is evaluated at :math:`t_{n-1}` instead of :math:`t_n`. Dividing by :math:`\Delta t^2`, the left-hand side is an approximation to :math:`u^{\prime\prime}` at :math:`t_n`, while the right-hand side is sampled at :math:`t_{n-1}`. All terms should be sampled at the same mesh point, so using :math:`\omega^2 u^{n-1}` instead of :math:`\omega^2 u^n` is an inconsistency in the scheme. This inconsistency turns out to be rather crucial for the accuracy of the Forward Euler method applied to vibration problems. The Backward Euler scheme ------------------------- A Backward Euler approximation the ODE system is equally easy to write up in the operator notation: .. _Eq:_auto10: .. math:: \tag{29} \lbrack D_t^- u = v\rbrack^{n+1}, .. _Eq:_auto11: .. math:: \tag{30} \lbrack D_t^- v = -\omega u\rbrack^{n+1} {\thinspace .} This becomes a coupled system for :math:`u^{n+1}` and :math:`v^{n+1}`: .. _Eq:vib:undamped:BE1: .. math:: \tag{31} u^{n+1} - \Delta t v^{n+1} = u^{n}, .. _Eq:vib:undamped:BE2: .. math:: \tag{32} v^{n+1} + \Delta t \omega^2 u^{n+1} = v^{n} {\thinspace .} We can compare :ref:`(31) `-:ref:`(32) ` with the centered scheme :ref:`(7) ` for the second-order differential equation. To this end, we eliminate :math:`v^{n+1}` in :ref:`(31) ` using :ref:`(32) ` solved with respect to :math:`v^{n+1}`. Thereafter, we eliminate :math:`v^n` using :ref:`(31) ` solved with respect to :math:`v^{n+1}` and replacing :math:`n+1` by :math:`n`. The resulting equation involving only :math:`u^{n+1}`, :math:`u^n`, and :math:`u^{n-1}` can be ordered as .. math:: \frac{u^{n+1}-2u^n+u^{n-1}}{\Delta t^2} = -\omega^2 u^{n+1}, which has almost the same form as the centered scheme for the second-order differential equation, but the right-hand side is evaluated at :math:`u^{n+1}` and not :math:`u^n`. This inconsistent sampling of terms has a dramatic effect on the numerical solution. The Crank-Nicolson scheme ------------------------- The Crank-Nicolson scheme takes this form in the operator notation: .. _Eq:_auto12: .. math:: \tag{33} \lbrack D_t u = \overline{v}^t\rbrack^{n+\frac{1}{2}}, .. _Eq:_auto13: .. math:: \tag{34} \lbrack D_t v = -\omega \overline{u}^t\rbrack^{n+\frac{1}{2}} {\thinspace .} Writing the equations out shows that this is also a coupled system: .. _Eq:_auto14: .. math:: \tag{35} u^{n+1} - \frac{1}{2}\Delta t v^{n+1} = u^{n} + \frac{1}{2}\Delta t v^{n}, .. _Eq:_auto15: .. math:: \tag{36} v^{n+1} + \frac{1}{2}\Delta t \omega^2 u^{n+1} = v^{n} - \frac{1}{2}\Delta t \omega^2 u^{n} {\thinspace .} To see the nature of this approximation, and that it is actually very promising, we write the equations as follows .. _Eq:vib:undamped:CN3a: .. math:: \tag{37} u^{n+1} - u^n = \frac{1}{2}\Delta t(v^{n+1} + v^n), .. _Eq:vib:undamped:CN4a: .. math:: \tag{38} v^{n+1} = v^n -\frac{1}{2}\Delta t(u^{n+1} + u^n), and add the latter at the previous time level as well: .. _Eq:vib:undamped:CN4b1: .. math:: \tag{39} v^{n} = v^{n-1} -\frac{1}{2}\Delta t(u^{n} + u^{n-1}) We can also rewrite :ref:`(37) ` at the previous time level as .. _Eq:vib:undamped:CN4b: .. math:: \tag{40} v^{n+1} + v^n = \frac{2}{\Delta t}(u^{n+1} - u^n){\thinspace .} Inserting :ref:`(38) ` for :math:`v^{n+1}` in :ref:`(37) ` and :ref:`(39) ` for :math:`v^{n}` in :ref:`(37) ` yields after some reordering: .. math:: u^{n+1} - n^n = \frac{1}{2}(-\frac{1}{2}\Delta t\omega^2 (u^{n+1} + 2u^n + u^{n-1}) + v^ + v^{n-1}){\thinspace .} Now, :math:`v^n + v^{n-1}` can be eliminated by means of :ref:`(40) `. The result becomes .. _Eq:vib:undamped:CN5: .. math:: \tag{41} u^{n+1} - 2u^n + u^{n-1} = \Delta t^2\omega^2 \frac{1}{4}(u^{n+1} + 2u^n + u^{n-1}){\thinspace .} We have that .. math:: \frac{1}{4}(u^{n+1} + 2u^n + u^{n-1}) \approx u^n + {\mathcal{O}(\Delta t^2)}, meaning that :ref:`(41) ` is an approximation to the centered scheme :ref:`(7) ` for the second-order ODE where the sampling error in the term :math:`\Delta t^2\omega^2 u^n` is of the same order as the approximation errors in the finite differences, i.e., :math:`{\mathcal{O}(\Delta t^2)}`. The Crank-Nicolson scheme written as :ref:`(41) ` therefore has consistent sampling of all terms at the same time point :math:`t_n`. The implication is a much better method than the Forward and Backward Euler schemes. .. _vib:model2x2:compare: Comparison of schemes --------------------- We can easily compare methods like the ones above (and many more!) with the aid of the `Odespy `__ package. Below is a sketch of the code. .. code-block:: python import odespy import numpy as np def f(u, t, w=1): u, v = u # u is array of length 2 holding our [u, v] return [v, -w**2*u] def run_solvers_and_plot(solvers, timesteps_per_period=20, num_periods=1, I=1, w=2*np.pi): P = 2*np.pi/w # duration of one period dt = P/timesteps_per_period Nt = num_periods*timesteps_per_period T = Nt*dt t_mesh = np.linspace(0, T, Nt+1) legends = [] for solver in solvers: solver.set(f_kwargs={'w': w}) solver.set_initial_condition([I, 0]) u, t = solver.solve(t_mesh) There is quite some more code dealing with plots also, and we refer to the source file `vib_undamped_odespy.py `__ for details. Observe that keyword arguments in ``f(u,t,w=1)`` can be supplied through a solver parameter ``f_kwargs`` (dictionary of additional keyword arguments to ``f``). Specification of the Forward Euler, Backward Euler, and Crank-Nicolson schemes is done like this: .. code-block:: python solvers = [ odespy.ForwardEuler(f), # Implicit methods must use Newton solver to converge odespy.BackwardEuler(f, nonlinear_solver='Newton'), odespy.CrankNicolson(f, nonlinear_solver='Newton'), ] .. index:: phase plane plot The ``vib_undamped_odespy.py`` program makes two plots of the computed solutions with the various methods in the ``solvers`` list: one plot with :math:`u(t)` versus :math:`t`, and one *phase plane plot* where :math:`v` is plotted against :math:`u`. That is, the phase plane plot is the curve :math:`(u(t),v(t))` parameterized by :math:`t`. Analytically, :math:`u=I\cos(\omega t)` and :math:`v=u^{\prime}=-\omega I\sin(\omega t)`. The exact curve :math:`(u(t),v(t))` is therefore an ellipse, which often looks like a circle in a plot if the axes are automatically scaled. The important feature, however, is that exact curve :math:`(u(t),v(t))` is closed and repeats itself for every period. Not all numerical schemes are capable of doing that, meaning that the amplitude instead shrinks or grows with time. Figure :ref:`vib:ode1:1st:odespy:theta:phaseplane` show the results. Note that Odespy applies the label MidpointImplicit for what we have specified as ``CrankNicolson`` in the code (``CrankNicolson`` is just a synonym for class ``MidpointImplicit`` in the Odespy code). The Forward Euler scheme in Figure :ref:`vib:ode1:1st:odespy:theta:phaseplane` has a pronounced spiral curve, pointing to the fact that the amplitude steadily grows, which is also evident in Figure :ref:`vib:ode1:1st:odespy:theta`. The Backward Euler scheme has a similar feature, except that the spriral goes inward and the amplitude is significantly damped. The changing amplitude and the sprial form decreases with decreasing time step. The Crank-Nicolson scheme looks much more accurate. In fact, these plots tell that the Forward and Backward Euler schemes are not suitable for solving our ODEs with oscillating solutions. .. _vib:ode1:1st:odespy:theta:phaseplane: .. figure:: vib_theta_1_pp.png :width: 800 *Comparison of classical schemes in the phase plane for two time step values* .. _vib:ode1:1st:odespy:theta: .. figure:: vib_theta_1_u.png :width: 800 *Comparison of solution curves for classical schemes* Runge-Kutta methods ------------------- We may run two popular standard methods for first-order ODEs, the 2nd- and 4th-order Runge-Kutta methods, to see how they perform. Figures :ref:`vib:ode1:1st:odespy:RK:phaseplane` and :ref:`vib:ode1:1st:odespy:RK` show the solutions with larger :math:`\Delta t` values than what was used in the previous two plots. .. _vib:ode1:1st:odespy:RK:phaseplane: .. figure:: vib_RK_1_pp.png :width: 800 *Comparison of Runge-Kutta schemes in the phase plane* .. _vib:ode1:1st:odespy:RK: .. figure:: vib_RK_1_u.png :width: 800 *Comparison of Runge-Kutta schemes* The visual impression is that the 4th-order Runge-Kutta method is very accurate, under all circumstances in these tests, while the 2nd-order scheme suffers from amplitude errors unless the time step is very small. The corresponding results for the Crank-Nicolson scheme are shown in Figure :ref:`vib:ode1:1st:odespy:CN:long:phaseplane`. It is clear that the Crank-Nicolson scheme outperforms the 2nd-order Runge-Kutta method. Both schemes have the same order of accuracy :math:`{\mathcal{O}(\Delta t^2)}`, but their differences in the accuracy that matters in a real physical application is very clearly pronounced in this example. :ref:`vib:exer:undamped:odespy` invites you to investigate how the amplitude is computed by a series of famous methods for first-order ODEs. .. _vib:ode1:1st:odespy:CN:long:phaseplane: .. figure:: vib_CN_10_pp.png :width: 800 *Long-time behavior of the Crank-Nicolson scheme in the phase plane* Analysis of the Forward Euler scheme ------------------------------------ We may try to find exact solutions of the discrete equations :ref:`(27) `-:ref:`(28) ` in the Forward Euler method. An "ansatz" is .. math:: u^n &= IA^n,\\ v^n &= qIA^n, where :math:`q` and :math:`A` are unknown numbers. We could have used a complex exponential form :math:`e^{i\tilde\omega n\Delta t}` since we get oscillatory form, but the oscillations grow in the Forward Euler method, so the numerical frequency :math:`\tilde\omega` will be complex anyway (producing an exponentially growing amplitude). Therefore, it is easier to just work with potentially complex :math:`A` and :math:`q` as introduced above. The Forward Euler scheme leads to .. math:: A &= 1 + \Delta t q,\\ A &= 1 - \Delta t\omega^2 q^{-1}{\thinspace .} We can easily eliminate :math:`A`, get :math:`q^2 + \omega^2=0`, and solve for .. math:: q = \pm i\omega, which gives .. math:: A = 1 \pm \Delta t i\omega{\thinspace .} We shall take the real part of :math:`A^n` as the solution. The two values of :math:`A` are complex conjugates, and the real part of :math:`A^n` will be the same for both roots. This is easy to realize if we rewrite the complex numbers in polar form, which is also convenient for further analysis and understanding. The polar form :math:`re^{i\theta}` of a complex number :math:`x+iy` has :math:`r=\sqrt{x^2+y^2}` and :math:`\theta = \tan^{-1}(y/x)`. Hence, the polar form of the two values for :math:`A` become .. math:: 1 \pm \Delta t i\omega = \sqrt{1+\omega^2\Delta t^2}e^{\pm i\tan^{-1}(\omega\Delta t)}{\thinspace .} Now it is very easy to compute :math:`A^n`: .. math:: (1 \pm \Delta t i\omega)^n = (1+\omega^2\Delta t^2)^{n/2}e^{\pm ni\tan^{-1}(\omega\Delta t)}{\thinspace .} Since :math:`\cos (\theta n) = \cos (-\theta n)`, the real part of the two numbers become the same. We therefore continue with the solution that has the plus sign. The general solution is :math:`u^n = CA^n`, where :math:`C` is a constant determined from the initial condition: :math:`u^0=C=I`. We have :math:`u^n=IA^n` and :math:`v^n=qIA^n`. The final solutions are just the real part of the expressions in polar form: .. _Eq:_auto16: .. math:: \tag{42} u^n = I(1+\omega^2\Delta t^2)^{n/2}\cos (n\tan^{-1}(\omega\Delta t)), .. _Eq:_auto17: .. math:: \tag{43} v^n =- \omega I(1+\omega^2\Delta t^2)^{n/2}\sin (n\tan^{-1}(\omega\Delta t)){\thinspace .} The expression :math:`(1+\omega^2\Delta t^2)^{n/2}` causes growth of the amplitude, since a number greater than one is raised to a positive exponent :math:`n/2`. We can develop a series expression to better understand the formula for the amplitude. Introducing :math:`p=\omega\Delta t` as the key variable and using ``sympy`` gives .. code-block:: python >>> from sympy import * >>> p = symbols('p', real=True) >>> n = symbols('n', integer=True, positive=True) >>> amplitude = (1 + p**2)**(n/2) >>> amplitude.series(p, 0, 4) 1 + n*p**2/2 + O(p**4) The amplitude goes like :math:`1 + \frac{1}{2} n\omega^2\Delta t^2`, clearly growing linearly in time (with :math:`n`). We can also investigate the error in the angular frequency by a series expansion: .. code-block:: python >>> n*atan(p).series(p, 0, 4) n*(p - p**3/3 + O(p**4)) This means that the solution for :math:`u^n` can be written as .. math:: u^n = (1 + \frac{1}{2} n\omega^2\Delta t^2 + \Oof(\Delta t^4)) \cos\left(\omega t - \frac{1}{3}\omega t\Delta t^2 + {\mathcal{O}(\Delta t^4)}\right) {\thinspace .} The error in the angular frequency is of the same order as in the scheme :ref:`(7) ` for the second-order ODE, but error in the amplitude is severe. .. _vib:model1:energy: Energy considerations ===================== .. index:: mechanical energy .. index:: energy principle The observations of various methods in the previous section can be better interpreted if we compute a quantity reflecting the total *energy of the system*. It turns out that this quantity, .. math:: E(t) = \frac{1}{2}(u^{\prime})^2 + \frac{1}{2}\omega^2u^2, is *constant* for all :math:`t`. Checking that :math:`E(t)` really remains constant brings evidence that the numerical computations are sound. It turns out that :math:`E` is proportional to the mechanical energy in the system. Conservation of energy is much used to check numerical simulations. Derivation of the energy expression ----------------------------------- We start out with multiplying .. math:: u^{\prime\prime} + \omega^2 u = 0, by :math:`u^{\prime}` and integrating from :math:`0` to :math:`T`: .. math:: \int_0^T u^{\prime\prime}u^{\prime} dt + \int_0^T\omega^2 u u^{\prime} dt = 0{\thinspace .} Observing that .. math:: u^{\prime\prime}u^{\prime} = \frac{d}{dt}\frac{1}{2}(u^{\prime})^2,\quad uu^{\prime} = \frac{d}{dt} {\frac{1}{2}}u^2, we get .. math:: \int_0^T (\frac{d}{dt}\frac{1}{2}(u^{\prime})^2 + \frac{d}{dt} \frac{1}{2}\omega^2u^2)dt = E(T) - E(0)=0, where we have introduced .. _Eq:vib:model1:energy:balance1: .. math:: \tag{44} E(t) = \frac{1}{2}(u^{\prime})^2 + \frac{1}{2}\omega^2u^2{\thinspace .} The important result from this derivation is that the total energy is constant: .. math:: E(t) = E(0){\thinspace .} .. admonition:: :math:`E(t)` is closely related to the system's energy The quantity :math:`E(t)` derived above is physically not the mechanical energy of a vibrating mechanical system, but the energy per unit mass. To see this, we start with Newton's second law :math:`F=ma` (:math:`F` is the sum of forces, :math:`m` is the mass of the system, and :math:`a` is the acceleration). The displacement :math:`u` is related to :math:`a` through :math:`a=u^{\prime\prime}`. With a spring force as the only force we have :math:`F=-ku`, where :math:`k` is a spring constant measuring the stiffness of the spring. Newton's second law then implies the differential equation .. math:: -ku = mu^{\prime\prime}\quad\Rightarrow mu^{\prime\prime} + ku = 0{\thinspace .} This equation of motion can be turned into an energy balance equation by finding the work done by each term during a time interval :math:`[0,T]`. To this end, we multiply the equation by :math:`du=u^{\prime}dt` and integrate: .. math:: \int_0^T muu^{\prime}dt + \int_0^T kuu^{\prime}dt = 0{\thinspace .} The result is .. math:: \tilde E(t) = E_k(t) + E_p(t) = 0, where .. _Eq:vib:model1:energy:kinetic: .. math:: \tag{45} E_k(t) = \frac{1}{2}mv^2,\quad v=u^{\prime}, is the *kinetic energy* of the system, and .. _Eq:vib:model1:energy:potential: .. math:: \tag{46} E_p(t) = {\frac{1}{2}}ku^2 is the *potential energy*. The sum :math:`\tilde E(t)` is the total mechanical energy. The derivation demonstrates the famous energy principle that, under the right physical circumstances, any change in the kinetic energy is due to a change in potential energy and vice versa. (This principle breaks down when we introduce damping in system, as we do in the section :ref:`vib:model2`.) The equation :math:`mu^{\prime\prime}+ku=0` can be divided by :math:`m` and written as :math:`u^{\prime\prime} + \omega^2u=0` for :math:`\omega=\sqrt{k/m}`. The energy expression :math:`E(t)=\frac{1}{2}(u^{\prime})^2 + \frac{1}{2}\omega^2u^2` derived earlier is then :math:`\tilde E(t)/m`, i.e., mechanical energy per unit mass. Energy of the exact solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Analytically, we have :math:`u(t)=I\cos\omega t`, if :math:`u(0)=I` and :math:`u^{\prime}(0)=0`, so we can easily check that the energy evolution and confirm that :math:`E(t)` is constant: .. math:: E(t) = {\frac{1}{2}}I^2 (-\omega\sin\omega t)^2 + \frac{1}{2}\omega^2 I^2 \cos^2\omega t = \frac{1}{2}\omega^2 (\sin^2\omega t + \cos^2\omega t) = \frac{1}{2}\omega^2 {\thinspace .} An error measure based on energy -------------------------------- The constant energy is well expressed by its initial value :math:`E(0)`, so that the error in mechanical energy can be computed as a mesh function by .. _Eq:_auto18: .. math:: \tag{47} e_E^n = \frac{1}{2}\left(\frac{u^{n+1}-u^{n-1}}{2\Delta t}\right)^2 + \frac{1}{2}\omega^2 (u^n)^2 - E(0), \quad n=1,\ldots,N_t-1, where .. math:: E(0) = {\frac{1}{2}}V^2 + \frac{1}{2}\omega^2I^2, if :math:`u(0)=I` and :math:`u^{\prime}(0)=V`. Note that we have used a centered approximation to :math:`u^{\prime}`: :math:`\boldsymbol{u}^{\prime}(t_n)\approx [D_{2t}u]^n`. A useful norm of the mesh function :math:`e_E^n` for the discrete mechanical energy can be the maximum absolute value of :math:`e_E^n`: .. math:: ||e_E^n||_{\ell^\infty} = \max_{1\leq n `, it may be used directly in :ref:`(51) ` to compute :math:`v^{n+1}`. In physics, it is more common to update the :math:`v` equation first, with a forward difference, and thereafter the :math:`u` equation, with a backward difference that applies the most recently computed :math:`v` value: .. _Eq:vib:model2x2:EulerCromer:veq1b: .. math:: \tag{52} v^{n+1} = v^n -\Delta t \omega^2u^{n}, .. _Eq:vib:model2x2:EulerCromer:ueq1b: .. math:: \tag{53} u^{n+1} = u^n + \Delta t v^{n+1}{\thinspace .} The advantage of ordering the ODEs as in :ref:`(52) `-:ref:`(53) ` becomes evident when consider complicated models. Such models are included if we write our vibration ODE more generally as .. math:: \ddot u + g(u, u^{\prime}, t)=0{\thinspace .} We can rewrite this second-order ODE as two first-order ODEs, .. math:: v' &= -g(u,v,t),\\ u^{\prime} &= v{\thinspace .} This rewrite allows the following scheme to be used: .. math:: v^{n+1} &= v^n -\Delta t\, g(u^n,v^n,t),\\ u^{n+1} &= u^n + \Delta t\, v^{n+1}{\thinspace .} We realize that the first update works well with any :math:`g` since old values :math:`u^n` and :math:`v^n` are used. Switching the equations would demand :math:`u^n{+1}` and :math:`v^{n+1}` values in :math:`g`. .. Despite using a backward difference, there is no need to solve a coupled .. system for :math:`u^{n+1}` and :math:`v^{n+1}` because the structure of the ODEs .. allows :ref:`(50) ` The scheme :ref:`(52) `-:ref:`(53) ` goes under several names: forward-backward scheme, `semi-implicit Euler method `__, semi-explicit Euler, symplectic Euler, Newton-Stormer-Verlet, and Euler-Cromer. We shall stick to the latter name. Since both time discretizations are based on first-order difference approximation, one may think that the scheme is only of first-order, but this is not true: the use of a forward and then a backward difference make errors cancel so that the overall error in the scheme is :math:`{\mathcal{O}(\Delta t^2)}`. This is explained below. .. _vib:model2x2:EulerCromer:equiv: Equivalence with the scheme for the second-order ODE ---------------------------------------------------- We may eliminate the :math:`v^n` variable from :ref:`(50) `-:ref:`(51) ` or :ref:`(52) `-:ref:`(53) `. The :math:`v^{n+1}` term in :ref:`(52) ` can be eliminated from :ref:`(53) `: .. _Eq:vib:model2x2:EulerCromer:elim1: .. math:: \tag{54} u^{n+1} = u^n + \Delta t (v^n - \omega^2\Delta t^2 u^n){\thinspace .} The :math:`v^{n}` quantity can be expressed by :math:`u^n` and :math:`u^{n-1}` using :ref:`(53) `: .. math:: v^{n} = \frac{u^n - u^{n-1}}{\Delta t}, and when this is inserted in :ref:`(54) ` we get .. _Eq:_auto21: .. math:: \tag{55} u^{n+1} = 2u^n - u^{n-1} - \Delta t^2 \omega^2u^{n}, which is nothing but the centered scheme :ref:`(7) `! The two seemingly different numerical methods are mathematically equivalent. Consequently, the previous analysis of :ref:`(7) ` also applies to the Euler-Cromer method. In particular, the amplitude is constant, given that the stability criterion is fulfilled, but there is always an angular frequency error :ref:`(19) `. :ref:`vib:exer:EulerCromer:analysis` gives guidance on how to derive the exact discrete solution of the two equations in the Euler-Cromer method. Although the Euler-Cromer scheme and the method :ref:`(7) ` are equivalent, there could be differences in the way they handle the initial conditions. Let is look into this topic. The initial condition :math:`u^{\prime}=0` means :math:`u^{\prime}=v=0`. From :ref:`(53) ` we get :math:`v^1=-\omega^2 u^0` and :math:`u^1=u^0 - \omega^2\Delta t^2 u^0`. When using a centered approximation of :math:`u^{\prime}(0)=0` combined with the discretization :ref:`(7) ` of the second-order ODE, we get :math:`u^1=u^0 - \frac{1}{2}\omega^2\Delta t^2 u^0`. The difference is :math:`\frac{1}{2}\omega^2\Delta t^2 u^0`, which is of second order in :math:`\Delta t`, seemingly consistent with the overall error in the scheme for the differential equation model. A different view can also be taken. If we approximate :math:`u^{\prime}(0)=0` by a backward difference, :math:`(u^0-u^{-1})/\Delta t =0`, we get :math:`u^{-1}=u^0`, and when combined with :ref:`(7) `, it results in :math:`u^1=u^0 - \omega^2\Delta t^2 u^0`. This means that the Euler-Cromer method based on :ref:`(53) `-:ref:`(52) ` corresponds to using only a first-order approximation to the initial condition in the method from the section :ref:`vib:ode1:fdm`. Correspondingly, using the formulation :ref:`(50) `-:ref:`(51) ` with :math:`v^n=0` leads to :math:`u^1=u^0`, which can be interpreted as using a forward difference approximation for the initial condition :math:`u^{\prime}(0)=0`. Both Euler-Cromer formulations lead to slightly different values for :math:`u^1` compared to the method in the section :ref:`vib:ode1:fdm`. The error is :math:`\frac{1}{2}\omega^2\Delta t^2 u^0` and of the same order as the overall scheme. .. _vib:model2x2:EulerCromer:impl: Implementation (2) --------------------------- The function below, found in `vib_EulerCromer.py `__ implements the Euler-Cromer scheme :ref:`(52) `-:ref:`(53) `: .. code-block:: python import numpy as np def solver(I, w, dt, T): """ Solve v' = - w**2*u, u'=v for t in (0,T], u(0)=I and v(0)=0, by an Euler-Cromer method. """ dt = float(dt) Nt = int(round(T/dt)) u = np.zeros(Nt+1) v = np.zeros(Nt+1) t = np.linspace(0, Nt*dt, Nt+1) v[0] = 0 u[0] = I for n in range(0, Nt): v[n+1] = v[n] - dt*w**2*u[n] u[n+1] = u[n] + dt*v[n+1] return u, v, t Since the Euler-Cromer scheme is equivalent to the finite difference method for the second-order ODE :math:`u^{\prime\prime}+\omega^2u=0` (see the section :ref:`vib:model2x2:EulerCromer:equiv`), the performance of the above ``solver`` function is the same as for the ``solver`` function in the section :ref:`vib:impl1`. The only difference is the formula for the first time step, as discussed above. This deviation in the Euler-Cromer scheme means that the discrete solution listed in the section :ref:`vib:ode1:analysis:sol` is not a solution of the Euler-Cromer scheme! To verify the implementation of the Euler-Cromer method we can adjust ``v[1]`` so that the computer-generated values can be compared with the formula :ref:`(20) ` from in the section :ref:`vib:ode1:analysis:sol`. This adjustment is done in an alternative solver function, ``solver_ic_fix`` in ``vib_EulerCromer.py``. Since we now have an exact solution of the discrete equations available, we can write a test function ``test_solver`` for checking the equality of computed values with the formula :ref:`(20) `: .. code-block:: python def test_solver(): """ Test solver with fixed initial condition against equivalent scheme for the 2nd-order ODE u'' + u = 0. """ I = 1.2; w = 2.0; T = 5 dt = 2/w # longest possible time step u, v, t = solver_ic_fix(I, w, dt, T) from vib_undamped import solver as solver2 # 2nd-order ODE u2, t2 = solver2(I, w, dt, T) error = np.abs(u - u2).max() tol = 1E-14 assert error < tol Another function, ``demo``, visualizes the difference between Euler-Cromer scheme and the scheme :ref:`(7) ` for the second-oder ODE, arising from the mismatch in the first time level. .. is anything gained? is v of higher order than D_2t u from the .. other approach, i.e., if we need v, is this alg better? Probably not .. since v is related u through a difference .. make exercises: .. investigate how important the u^1 wrong formula really is on .. convergence rate .. new file: genealizations, systems, .. new file: apps .. exercise: damping analysis, see geophysics book first... The velocity Verlet algorithm ----------------------------- Another very popular algorithm for vibration problems :math:`u^{\prime\prime}+\omega^2u=0` can be derived as follows. First, we step :math:`u` forward from :math:`t_n` to :math:`t_{n+1}` using a three-term Taylor series, .. math:: u(t_{n+1}) = u(t_n) + u^{\prime}(t_n)\Delta t + \frac{1}{2}u^{\prime\prime}(t_n)\Delta t^2{\thinspace .} Using :math:`u^{\prime}=v` and :math:`u^{\prime\prime}=-\omega^2u`, we get the updating formula .. math:: u^{n+1} = u^n + v^n\Delta t - \frac{1}{2}\Delta^2\omega^2u^n{\thinspace .} Second, the first-order equation for :math:`v`, .. math:: v'=-\omega^2u, is discretized by a centered difference in a Crank-Nicolson fashion at :math:`t_{n+\frac{1}{2}}`: .. math:: \frac{v^{n+1}-v^n}{\Delta t} = -\omega^2\frac{1}{2}(u^n + u^{n+1}){\thinspace .} To summarize, we have the scheme .. _Eq:vib:model2x2:Verlet:dueq: .. math:: \tag{56} u^{n+1} = u^n + v^n\Delta t - \frac{1}{2}\Delta^2\omega^2u^n .. _Eq:vib:model2x2:Verlet:dveq: .. math:: \tag{57} v^{n+1} = v^n -\frac{1}{2}\Delta t\omega^2 (u^n + u^{n+1}), known as the *velocity Verlet* algorithm. Observe that this scheme is explicit since :math:`u^{n+1}` in :ref:`(57) ` is already computed from :ref:`(56) `. The algorithm can be straightforwardly implemented as shown below (the code appears in the file `vib_undamped_velocity_Verlet.py `__). .. code-block:: python from vib_undamped import convergence_rates, main def solver(I, w, dt, T, return_v=False): """ Solve u'=v, v'=-w**2*u for t in (0,T], u(0)=I and v(0)=0, by the velocity Verlet method with time step dt. """ dt = float(dt) Nt = int(round(T/dt)) u = np.zeros(Nt+1) v = np.zeros(Nt+1) t = np.linspace(0, Nt*dt, Nt+1) u[0] = I v[0] = 0 for n in range(Nt): u[n+1] = u[n] + v[n]*dt - 0.5*dt**2*w**2*u[n] v[n+1] = v[n] - 0.5*dt*w**2*(u[n] + u[n+1]) if return_v: return u, v, t else: # Return just u and t as in the vib_undamped.py's solver return u, t We provide the option that this ``solver`` function returns the same data as the ``solver`` function from the section :ref:`vib:impl1:solver` (if ``return_v`` is ``False``), but we may return ``v`` along with ``u`` and ``t``. The error in the Taylor series expansion behind :ref:`(56) ` is :math:`{\mathcal{O}(\Delta t^3)}`, while the error in the central difference for :math:`v` is :math:`{\mathcal{O}(\Delta t^2)}`. The overall error is then no better than :math:`{\mathcal{O}(\Delta t^2)}`, which can be verified empirically using the ``convergence_rates`` function from :ref:`vib:ode1:verify`: .. code-block:: python >>> import vib_undamped_velocity_Verlet as m >>> m.convergence_rates(4, solver_function=m.solver) [2.0036366687367346, 2.0009497328124835, 2.000240105995295] .. The output confirms that the overall convergence rate is 2.