$$ \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]
Svein Linge [3, 1]

[1] Center for Biomedical Computing, Simula Research Laboratory
[2] Department of Informatics, University of Oslo
[3] Department of Process, Energy and Environmental Technology, University College of Southeast Norway

Jul 14, 2016










Table of contents

A simple vibration problem
      A centered finite difference scheme; step 1 and 2
      A centered finite difference scheme; step 3
      A centered finite difference scheme; step 4
      Computing the first step
      The computational algorithm
      Operator notation; ODE
      Operator notation; initial condition
      Computing \( u^{\prime} \)
Implementation
      Core algorithm
      Plotting
      Main program
      User interface: command line
      Running the program
Verification
      First steps for testing and debugging
      Checking convergence rates
      Implementational details
      Unit test for the convergence rate
Long time simulations
      Effect of the time step on long simulations
      Using a moving plot window
      Long time simulations visualized with aid of Bokeh: coupled panning of multiple graphs
      How does Bokeh plotting code look like?
Analysis of the numerical scheme
      Movie of the angular frequency error
      We can derive an exact solution of the discrete equations
      Calculations of an exact solution of the discrete equations
      Solving for the numerical frequency
      Polynomial approximation of the frequency error
      Plot of the frequency error
      Exact discrete solution
      Convergence of the numerical scheme
      Stability
      The stability criterion
      Summary of the analysis
Alternative schemes based on 1st-order equations
      Rewriting 2nd-order ODE as system of two 1st-order ODEs
      The Forward Euler scheme
      The Backward Euler scheme
      The Crank-Nicolson scheme
      Comparison of schemes via Odespy
      Forward and Backward Euler and Crank-Nicolson
      Phase plane plot of the numerical solutions
      Plain solution curves
      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
      Observations of RK and CN methods
      Energy conservation property
      Derivation of the energy conservation property
      Remark about \( E(t) \)
      The Euler-Cromer method; idea
      The Euler-Cromer method; complete formulas
      Euler-Cromer is equivalent to the scheme for \( u^{\prime\prime}+\omega^2u=0 \)
      The schemes are not equivalent wrt the initial conditions
Generalization: damping, nonlinear spring, and external excitation
      A centered scheme for linear damping
      Initial conditions
      Linearization via a geometric mean approximation
      A centered scheme for quadratic damping
      Initial condition for quadratic damping
      Algorithm
      Implementation
      Verification
      Demo program
      Euler-Cromer formulation
      Staggered grid
      Linear damping
      Quadratic damping
      Initial conditions









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: 

1
2
3
4
5
6
7
8
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 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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

def solver_adjust_w(I, w, dt, T, adjust_w=True):
    """
    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)
    w_adj = w*(1 - w**2*dt**2/24.) if adjust_w else w

    u[0] = I
    u[1] = u[0] - 0.5*dt**2*w_adj**2*u[0]
    for n in range(1, Nt):
        u[n+1] = 2*u[n] - u[n-1] - dt**2*w_adj**2*u[n]
    return u, t









Plotting 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
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 

1
2
3
4
5
6
7
8
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 

1
2
3
4
5
6
7
8
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: 

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

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

1
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 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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, E_values, dt_values

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









Unit test for the convergence rate

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def test_convergence_rates():
    r, E, dt = 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
    # Test that adjusted w obtains 4th order convergence
    r, E, dt = convergence_rates(
        m=5, solver_function=solver_adjust_w, num_periods=8)
    print 'adjust w rates:', r
    assert abs(r[-1] - 4.0) < tol

def plot_convergence_rates():
    r2, E2, dt2 = convergence_rates(
        m=5, solver_function=solver, num_periods=8)
    plt.loglog(dt2, E2)
    r4, E4, dt4 = convergence_rates(
        m=5, solver_function=solver_adjust_w, num_periods=8)
    plt.loglog(dt4, E4)
    plt.legend(['original scheme', r'adjusted $\omega$'],
               loc='upper left')
    plt.title('Convergence of finite difference methods')
    from plotslopes import slope_marker
    slope_marker((dt2[1], E2[1]), (2,1))
    slope_marker((dt4[1], E4[1]), (4,1))
    plt.savefig('tmp_convrate.png'); plt.savefig('tmp_convrate.pdf')
    plt.show()

Complete code in vib_undamped.py.









Long time simulations









Effect of the time step on long simulations









Using a moving plot window

Example: 

1
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? 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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: 

1
2
3
4
5
>>> 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: 

1
2
3
4
>>> 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, \label{vib:model2x2:ueq}\\ v^{\prime} &= -\omega^2 u \label{vib:model2x2:veq} \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, \label{_auto1}\\ v^{n+1} &= v^n -\Delta t \omega^2 u^n \tp \label{_auto2} \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}, \label{_auto3}\\ v^{n+1} + \Delta t \omega^2 u^{n+1} = v^{n} \tp \label{_auto4} \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}, \label{_auto5}\\ v^{n+1} + \half\Delta t \omega^2 u^{n+1} &= v^{n} - \half\Delta t \omega^2 u^{n} \tp \label{_auto6} \end{align} $$









Comparison of schemes via Odespy

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
import odespy
import numpy as np

def f(u, t, w=1):
    # v, u numbering for EulerCromer to work well
    v, u = u  # u is array of length 2 holding our [v, u]
    return [-w**2*u, v]

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([0, I])
        u, t = solver.solve(t_mesh)









Forward and Backward Euler and Crank-Nicolson 

1
2
3
4
5
6
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 \label{_auto7}\\ [D_t^-u &= v]^{n+1} \label{_auto8} \end{align} $$









The Euler-Cromer method; complete formulas

Written out: $$ \begin{align} u^0 &= I, \label{_auto9}\\ v^0 &= 0, \label{_auto10}\\ v^{n+1} &= v^n -\Delta t \omega^2u^{n} \label{vib:model2x2:EulerCromer:veq1}\\ u^{n+1} &= u^n + \Delta t v^{n+1} \label{vib:model2x2:EulerCromer:ueq1} \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 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
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: 

1
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*} $$
© 2016, Hans Petter Langtangen, Svein Linge. Released under CC Attribution 4.0 license