$$ \newcommand{\half}{\frac{1}{2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} $$

Study guide: Finite difference methods for vibration problems

Hans Petter Langtangen [1, 2]

[1] Center for Biomedical Computing, Simula Research Laboratory
[2] Department of Informatics, University of Oslo

Oct 12, 2015

A simple vibration problem

 
$$ u^{\prime\prime}(t) + \omega^2u = 0,\quad u(0)=I,\ u^{\prime}(0)=0,\ t\in (0,T] $$

 

Exact solution:

 
$$ u(t) = I\cos (\omega t) $$

 
\( u(t) \) oscillates with constant amplitude \( I \) and (angular) frequency \( \omega \). Period: \( P=2\pi/\omega \).

A centered finite difference scheme; step 1 and 2

 
$$ u^{\prime\prime}(t_n) + \omega^2u(t_n) = 0,\quad n=1,\ldots,N_t $$

 

A centered finite difference scheme; step 3

Step 3: Approximate derivative(s) by finite difference approximation(s). Very common (standard!) formula for \( u^{\prime\prime} \):

 
$$ u^{\prime\prime}(t_n) \approx \frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} $$

 

Use this discrete initial condition together with the ODE at \( t=0 \) to eliminate \( u^{-1} \):

 
$$ \frac{u^{n+1}-2u^n + u^{n-1}}{\Delta t^2} = -\omega^2 u^n $$

 

A centered finite difference scheme; step 4

Step 4: Formulate the computational algorithm. Assume \( u^{n-1} \) and \( u^n \) are known, solve for unknown \( u^{n+1} \):

 
$$ u^{n+1} = 2u^n - u^{n-1} - \Delta t^2\omega^2 u^n $$

 

Nick names for this scheme: Stormer's method or Verlet integration.

Computing the first step

Discretize \( u^{\prime}(0)=0 \) by a centered difference

 
$$ \frac{u^1-u^{-1}}{2\Delta t} = 0\quad\Rightarrow\quad u^{-1} = u^1 $$

 

Inserted in the scheme for \( n=0 \) gives

 
$$ u^1 = u^0 - \half \Delta t^2 \omega^2 u^0 $$

 

The computational algorithm

  1. \( u^0=I \)
  2. compute \( u^1 \)
  3. for \( n=1,2,\ldots,N_t-1 \):
    1. compute \( u^{n+1} \)

More precisly expressed in 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]

Note: w is consistently used for \( \omega \) in my code.

Operator notation; ODE

With \( [D_tD_t u]^n \) as the finite difference approximation to \( u^{\prime\prime}(t_n) \) we can write

 
$$ [D_tD_t u + \omega^2 u = 0]^n $$

 

\( [D_tD_t u]^n \) means applying a central difference with step \( \Delta t/2 \) twice:

 
$$ [D_t(D_t u)]^n = \frac{[D_t u]^{n+\half} - [D_t u]^{n-\half}}{\Delta t}$$

 
which is written out as

 
$$ \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} \tp $$

 

Operator notation; initial condition

 
$$ [u = I]^0,\quad [D_{2t} u = 0]^0 $$

 
where \( [D_{2t} u]^n \) is defined as

 
$$ [D_{2t} u]^n = \frac{u^{n+1} - u^{n-1}}{2\Delta t} \tp $$

 

Computing \( u^{\prime} \)

\( u \) is often displacement/position, \( u^{\prime} \) is velocity and can be computed by

 
$$ u^{\prime}(t_n) \approx \frac{u^{n+1}-u^{n-1}}{2\Delta t} = [D_{2t}u]^n $$

 

Implementation

Core algorithm

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

Plotting

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')

Main program

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)

User interface: command line

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

Running the program

vib_undamped.py:

Terminal> python vib_undamped.py --dt 0.05 --num_periods 40

Generates frames tmp_vib%04d.png in files. Can make movie:

Terminal> ffmpeg -r 12 -i tmp_vib%04d.png -c:v flv movie.flv

Can use avconv instead of ffmpeg.

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

Verification

First steps for testing and debugging

Checking convergence rates

The next function estimates convergence rates, i.e., it

Implementational details

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

Result: r contains values equal to 2.00 - as expected!

Unit test for the convergence rate

Use final r[-1] in a unit test:

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

Complete code in vib_undamped.py.

Long time simulations

Effect of the time step on long simulations

Using a moving plot window

Example:

Terminal> python vib_undamped.py --dt 0.05 --num_periods 40

Movie of the moving plot window.

!splot

Long time simulations visualized with aid of Bokeh: coupled panning of multiple graphs

!splot

How does Bokeh plotting code look like?

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)

Analysis of the numerical scheme

Can we understand the frequency error?

Movie of the angular frequency error

\( u^{\prime\prime} + \omega^2 u = 0 \), \( u(0)=1 \), \( u^{\prime}(0)=0 \), \( \omega=2\pi \), \( \uex(t)=\cos (2\pi t) \), \( \Delta t = 0.05 \) (20 intervals per period)

We can derive an exact solution of the discrete equations

Calculations of an exact solution of the discrete equations

 
$$ u^n = IA^n = I\exp{(\tilde\omega \Delta t\, n)}=I\exp{(\tilde\omega t)} = I\cos (\tilde\omega t) + iI\sin(\tilde \omega t) \tp $$

 

 
$$ \begin{align*} [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}\\ &= I\frac{\exp{(i\tilde\omega(t+\Delta t))} - 2\exp{(i\tilde\omega t)} + \exp{(i\tilde\omega(t-\Delta t))}}{\Delta t^2}\\ &= I\exp{(i\tilde\omega t)}\frac{1}{\Delta t^2}\left(\exp{(i\tilde\omega(\Delta t))} + \exp{(i\tilde\omega(-\Delta t))} - 2\right)\\ &= I\exp{(i\tilde\omega t)}\frac{2}{\Delta t^2}\left(\cosh(i\tilde\omega\Delta t) -1 \right)\\ &= I\exp{(i\tilde\omega t)}\frac{2}{\Delta t^2}\left(\cos(\tilde\omega\Delta t) -1 \right)\\ &= -I\exp{(i\tilde\omega t)}\frac{4}{\Delta t^2}\sin^2(\frac{\tilde\omega\Delta t}{2}) \end{align*} $$

 

Solving for the numerical frequency

The scheme with \( u^n=I\exp{(i\omega\tilde\Delta t\, n)} \) inserted gives

 
$$ -I\exp{(i\tilde\omega t)}\frac{4}{\Delta t^2}\sin^2(\frac{\tilde\omega\Delta t}{2}) + \omega^2 I\exp{(i\tilde\omega t)} = 0 $$

 
which after dividing by \( I\exp{(i\tilde\omega t)} \) results in

 
$$ \frac{4}{\Delta t^2}\sin^2(\frac{\tilde\omega\Delta t}{2}) = \omega^2 $$

 
Solve for \( \tilde\omega \):

 
$$ \tilde\omega = \pm \frac{2}{\Delta t}\sin^{-1}\left(\frac{\omega\Delta t}{2}\right) $$

 

Polynomial approximation of the frequency error

Taylor series expansion for small \( \Delta t \) gives a formula that is easier to understand:

>>> from sympy import *
>>> dt, w = symbols('dt w')
>>> w_tilde = asin(w*dt/2).series(dt, 0, 4)*2/dt
>>> print w_tilde
(dt*w + dt**3*w**3/24 + O(dt**4))/dt  # note the final "/dt"

 
$$ \tilde\omega = \omega\left( 1 + \frac{1}{24}\omega^2\Delta t^2\right) + {\cal O}(\Delta t^3) $$

 
The numerical frequency is too large (to fast oscillations).

Plot of the frequency error

Recommendation: 25-30 points per period.

Exact discrete solution

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

 

The error mesh function,

 
$$ e^n = \uex(t_n) - u^n = I\cos\left(\omega n\Delta t\right) - I\cos\left(\tilde\omega n\Delta t\right) $$

 
is ideal for verification and further analysis!

 
$$ e^n = I\cos\left(\omega n\Delta t\right) - I\cos\left(\tilde\omega n\Delta t\right) = -2I\sin\left(t\half\left( \omega - \tilde\omega\right)\right) \sin\left(t\half\left( \omega + \tilde\omega\right)\right) $$

 

Convergence of the numerical scheme

Can easily show convergence:

 
$$ e^n\rightarrow 0 \hbox{ as }\Delta t\rightarrow 0,$$

 
because

 
$$ \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:

>>> import sympy as sym
>>> dt, w = sym.symbols('x w')
>>> sym.limit((2/dt)*sym.asin(w*dt/2), dt, 0, dir='+')
w

Stability

Observations:

What is the consequence of complex \( \tilde\omega \)?

The stability criterion

Cannot tolerate growth and must therefore demand a stability criterion

 
$$ \frac{\omega\Delta t}{2} \leq 1\quad\Rightarrow\quad \Delta t \leq \frac{2}{\omega} $$

 

Try \( \Delta t = \frac{2}{\omega} + 9.01\cdot 10^{-5} \) (slightly too big!):

Summary of the analysis

We can draw three important conclusions:

  1. The key parameter in the formulas is \( p=\omega\Delta t \) (dimensionless)
    1. Period of oscillations: \( P=2\pi/\omega \)
    2. Number of time steps per period: \( N_P=P/\Delta t \)
    3. \( \Rightarrow\ p=\omega\Delta t = 2\pi/ N_P \sim 1/N_P \)
    4. The smallest possible \( N_P \) is 2 \( \Rightarrow \) $p\in (0,\pi]$

  2. For \( p\leq 2 \) the amplitude of \( u^n \) is constant (stable solution)
  3. \( u^n \) has a relative frequency error \( \tilde\omega/\omega \approx 1 + \frac{1}{24}p^2 \), making numerical peaks occur too early

Alternative schemes based on 1st-order equations

Rewriting 2nd-order ODE as system of two 1st-order ODEs

The vast collection of ODE solvers (e.g., in Odespy) cannot be applied to

 
$$ u^{\prime\prime} + \omega^2 u = 0$$

 
unless we write this higher-order ODE as a system of 1st-order ODEs.

Introduce an auxiliary variable \( v=u^{\prime} \):

 
$$ \begin{align} u^{\prime} &= v, \tag{1}\\ v^{\prime} &= -\omega^2 u \tag{2} \tp \end{align} $$

 

Initial conditions: \( u(0)=I \) and \( v(0)=0 \).

The Forward Euler scheme

We apply the Forward Euler scheme to each component equation:

 
$$ [D_t^+ u = v]^n,$$

 

 
$$ [D_t^+ v = -\omega^2 u]^n,$$

 
or written out,

 
$$ \begin{align} u^{n+1} &= u^n + \Delta t v^n, \tag{3}\\ v^{n+1} &= v^n -\Delta t \omega^2 u^n \tp \tag{4} \end{align} $$

 

The Backward Euler scheme

We apply the Backward Euler scheme to each component equation:

 
$$ [D_t^- u = v]^{n+1},$$

 

 
$$ [D_t^- v = -\omega u]^{n+1} \tp $$

 
Written out:

 
$$ \begin{align} u^{n+1} - \Delta t v^{n+1} = u^{n}, \tag{5}\\ v^{n+1} + \Delta t \omega^2 u^{n+1} = v^{n} \tp \tag{6} \end{align} $$

 
This is a coupled \( 2\times 2 \) system for the new values at \( t=t_{n+1} \)!

The Crank-Nicolson scheme

 
$$ [D_t u = \overline{v}^t]^{n+\half},$$

 

 
$$ [D_t v = -\omega \overline{u}^t]^{n+\half}$$

 
The result is also a coupled system:

 
$$ \begin{align} u^{n+1} - \half\Delta t v^{n+1} &= u^{n} + \half\Delta t v^{n}, \tag{7}\\ v^{n+1} + \half\Delta t \omega^2 u^{n+1} &= v^{n} - \half\Delta t \omega^2 u^{n} \tp \tag{8} \end{align} $$

 

Comparison of schemes via Odespy

Can use Odespy to compare many methods for first-order schemes:

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)

Forward and Backward Euler and Crank-Nicolson

solvers = [
    odespy.ForwardEuler(f),
    # Implicit methods must use Newton solver to converge
    odespy.BackwardEuler(f, nonlinear_solver='Newton'),
    odespy.CrankNicolson(f, nonlinear_solver='Newton'),
    ]

Two plot types:

Phase plane plot of the numerical solutions

Note: CrankNicolson in Odespy leads to the name MidpointImplicit in plots.

Plain solution curves


Figure 1: Comparison of classical schemes.

Observations from the figures

Runge-Kutta methods of order 2 and 4; short time series

Runge-Kutta methods of order 2 and 4; longer time series

Crank-Nicolson; longer time series

(MidpointImplicit means CrankNicolson in Odespy)

Observations of RK and CN methods

Energy conservation property

The model

 
$$ u^{\prime\prime} + \omega^2 u = 0,\quad u(0)=I,\ u^{\prime}(0)=V,$$

 
has the nice energy conservation property that

 
$$ E(t) = \half(u^{\prime})^2 + \half\omega^2u^2 = \hbox{const}\tp$$

 
This can be used to check solutions.

Derivation of the energy conservation property

Multiply \( u^{\prime\prime}+\omega^2u=0 \) by \( u^{\prime} \) and integrate:

 
$$ \int_0^T u^{\prime\prime}u^{\prime} dt + \int_0^T\omega^2 u u^{\prime} dt = 0\tp$$

 
Observing that

 
$$ u^{\prime\prime}u^{\prime} = \frac{d}{dt}\half(u^{\prime})^2,\quad uu^{\prime} = \frac{d}{dt} {\half}u^2,$$

 
we get

 
$$ \int_0^T (\frac{d}{dt}\half(u^{\prime})^2 + \frac{d}{dt} \half\omega^2u^2)dt = E(T) - E(0), $$

 
where

 
$$ E(t) = \half(u^{\prime})^2 + \half\omega^2u^2 $$

 

Remark about \( E(t) \)

\( E(t) \) does not measure energy, energy per mass unit.

Starting with an ODE coming directly from Newton's 2nd law \( F=ma \) with a spring force \( F=-ku \) and \( ma=mu^{\prime\prime} \) (\( a \): acceleration, \( u \): displacement), we have

 
$$ mu^{\prime\prime} + ku = 0$$

 
Integrating this equation gives a physical energy balance:

 
$$ E(t) = \underbrace{{\half}mv^2}_{\hbox{kinetic energy} } + \underbrace{{\half}ku^2}_{\hbox{potential energy}} = E(0),\quad v=u^{\prime} $$

 
Note: the balance is not valid if we add other terms to the ODE.

The Euler-Cromer method; idea

2x2 system for \( u^{\prime\prime}+\omega^2u=0 \):

 
$$ \begin{align*} v^{\prime} &= -\omega^2u\\ u^{\prime} &= v \end{align*} $$

 

Forward-backward discretization:

 
$$ \begin{align} [D_t^+v &= -\omega^2u]^n \tag{9}\\ [D_t^-u &= v]^{n+1} \tag{10} \end{align} $$

 

The Euler-Cromer method; complete formulas

Written out:

 
$$ \begin{align} u^0 &= I, \tag{11}\\ v^0 &= 0, \tag{12}\\ v^{n+1} &= v^n -\Delta t \omega^2u^{n} \tag{13}\\ u^{n+1} &= u^n + \Delta t v^{n+1} \tag{14} \end{align} $$

 

Names: Forward-backward scheme, Semi-implicit Euler method, symplectic Euler, semi-explicit Euler, Newton-Stormer-Verlet, and Euler-Cromer.

Euler-Cromer is equivalent to the scheme for \( u^{\prime\prime}+\omega^2u=0 \)

We can eliminate \( v^n \) and \( v^{n+1} \), resulting in

 
$$ u^{n+1} = 2u^n - u^{n-1} - \Delta t^2 \omega^2u^{n} $$

 

which is the centered finite differrence scheme for \( u^{\prime\prime}+\omega^2u=0 \)!

The schemes are not equivalent wrt the initial conditions

 
$$ u^{\prime}=v=0\quad\Rightarrow\quad v^0=0,$$

 
so

 
$$ \begin{align*} v^1 &= v^0 - \Delta t\omega^2 u^0 = - \Delta t\omega^2 u^0\\ u^1 &= u^0 + \Delta t v^1 = u^0 - \Delta t\omega^2 u^0 != \underbrace{u^0 - \frac{1}{2}\Delta t\omega^2 u^0}_{\mbox{from }[D_tD_t u +\omega^2 u=0]^n\mbox{ and }[D_{2t}u=0]^0} \end{align*} $$

 

The exact discrete solution derived earlier does not fit the Euler-Cromer scheme because of mismatch for \( u^1 \).

Generalization: damping, nonlinear spring, and external excitation

 
$$ mu'' + f(u') + s(u) = F(t),\quad u(0)=I,\ u'(0)=V,\ t\in (0,T] $$

 
Input data: \( m \), \( f(u') \), \( s(u) \), \( F(t) \), \( I \), \( V \), and \( T \).

Typical choices of \( f \) and \( s \):

A centered scheme for linear damping

 
$$ [mD_tD_t u + f(D_{2t}u) + s(u) = F]^n $$

 
Written out

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

 
Assume \( f(u') \) is linear in \( u'=v \):

 
$$ u^{n+1} = \left(2mu^n + (\frac{b}{2}\Delta t - m)u^{n-1} + \Delta t^2(F^n - s(u^n)) \right)(m + \frac{b}{2}\Delta t)^{-1} $$

 

Initial conditions

\( u(0)=I \), \( u'(0)=V \):

 
$$ \begin{align*} \lbrack u &=I\rbrack^0\quad\Rightarrow\quad u^0=I\\ \lbrack D_{2t}u &=V\rbrack^0\quad\Rightarrow\quad u^{-1} = u^{1} - 2\Delta t V \end{align*} $$

 
End result:

 
$$ u^1 = u^0 + \Delta t\, V + \frac{\Delta t^2}{2m}(-bV - s(u^0) + F^0) $$

 
Same formula for \( u^1 \) as when using a centered scheme for \( u''+\omega u=0 \).

Linearization via a geometric mean approximation

In general, the geometric mean approximation reads

 
$$ (w^2)^n \approx w^{n-\half}w^{n+\half}\tp$$

 
For \( |u'|u' \) at \( t_n \):

 
$$ [u'|u'|]^n \approx u'(t_n+{\half})|u'(t_n-{\half})|\tp$$

 
For \( u' \) at \( t_{n\pm 1/2} \) we use centered difference:

 
$$ u'(t_{n+1/2})\approx [D_t u]^{n+\half},\quad u'(t_{n-1/2})\approx [D_t u]^{n-\half} $$

 

A centered scheme for quadratic damping

After some algebra:

 
$$ \begin{align*} u^{n+1} &= \left( m + b|u^n-u^{n-1}|\right)^{-1}\times \\ & \qquad \left(2m u^n - mu^{n-1} + bu^n|u^n-u^{n-1}| + \Delta t^2 (F^n - s(u^n)) \right) \end{align*} $$

 

Initial condition for quadratic damping

Simply use that \( u'=V \) in the scheme when \( t=0 \) (\( n=0 \)):

 
$$ [mD_tD_t u + bV|V| + s(u) = F]^0 $$

 

which gives

 
$$ u^1 = u^0 + \Delta t V + \frac{\Delta t^2}{2m}\left(-bV|V| - s(u^0) + F^0\right) $$

 

Algorithm

  1. \( u^0=I \)
  2. compute \( u^1 \) (formula depends on linear/quadratic damping)
  3. for \( n=1,2,\ldots,N_t-1 \):
    1. compute \( u^{n+1} \) from formula (depends on linear/quadratic damping)

Implementation

def solver(I, V, m, b, s, F, dt, T, damping='linear'):
    dt = float(dt); b = float(b); m = float(m) # avoid integer div.
    Nt = int(round(T/dt))
    u = zeros(Nt+1)
    t = 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

Verification

Demo program

vib.py supports input via the command line:

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.

Euler-Cromer formulation

We rewrite

 
$$ mu'' + f(u') + s(u) = F(t),\quad u(0)=I,\ u'(0)=V,\ t\in (0,T] $$

 
as a first-order ODE system

 
$$ \begin{align*} u' &= v \\ v' &= m^{-1}\left(F(t) - f(v) - s(u)\right) \end{align*} $$

 

Staggered grid

 
$$ \begin{align*} \lbrack D_t u &= v\rbrack^{n-\half} \\ \lbrack D_tv &= m^{-1}\left(F(t) - f(v) - s(u)\right)\rbrack^n \end{align*} $$

 

Written out,

 
$$ \begin{align*} \frac{u^n - u^{n-1}}{\Delta t} &= v^{n-\half}\\ \frac{v^{n+\half} - v^{n-\half}}{\Delta t} &= m^{-1}\left(F^n - f(v^n) - s(u^n)\right) \end{align*} $$

 

Problem: \( f(v^n) \)

Linear damping

With \( f(v)=bv \), we can use an arithmetic mean for \( bv^n \) a la Crank-Nicolson schemes.

 
$$ \begin{align*} u^n & = u^{n-1} + {\Delta t}v^{n-\half},\\ v^{n+\half} &= \left(1 + \frac{b}{2m}\Delta t\right)^{-1}\left( v^{n-\half} + {\Delta t} m^{-1}\left(F^n - {\half}f(v^{n-\half}) - s(u^n)\right)\right)\tp \end{align*} $$

 

Quadratic damping

With \( f(v)=b|v|v \), we can use a geometric mean

 
$$ b|v^n|v^n\approx b|v^{n-\half}|v^{n+\half}, $$

 
resulting in

 
$$ \begin{align*} u^n & = u^{n-1} + {\Delta t}v^{n-\half},\\ v^{n+\half} &= (1 + \frac{b}{m}|v^{n-\half}|\Delta t)^{-1}\left( v^{n-\half} + {\Delta t} m^{-1}\left(F^n - s(u^n)\right)\right)\tp \end{align*} $$

 

Initial conditions

 
$$ \begin{align*} u^0 &= I\\ v^{\half} &= V - \half\Delta t\omega^2I \end{align*} $$