Lecture 05 - Time integration

SSP and symplectic methods

Georgii Oblapenko

Objectives

In this lecture, we will learn about advanced methods for solving time-dependent PDEs and ODEs. In particular, we will be discussing

  1. Stability
  2. Monotonicity
  3. “Structure” preservation

ODEs: why?

Want to solve ODE

\[ d\frac{u}{dt} = f(u,t),\quad u \in \mathbb{R}^n, \: t > 0, \quad u(0) = u^0. \]

ODEs come from - direct description of system (i.e. biology, chemical reactors) - semi-discretization of system (i.e. \(u\) — vector of cell values from a FVM discretization, RHS are the fluxes + sources)

Note

Method of lines: discretize in space, solve in time.

Take PDE \[ \partial_t{u(x,t)} + \sum_{\alpha} \partial_{\alpha} a_{\alpha} u(x,t) = f(u(x,t),t), \]

discretize in space (FD/FVM/FEM/DG): \(u(x_i,t)=u_i(t)\), discretize operators in terms of \(\{u_j\}\), obtain ODEs \[ \partial_t{\mathbf{u}} + \mathbf{A} \mathbf{u}(t) = \mathbf{f}(t). \]

PDE discretization - remark

Alternatives to Method of Lines

Are all discretizations of time-dependent PDEs in time done via Method of Lines?

Lax-Wendroff: couples time and space discretization. For linear advection \[ \partial_t u + a_0 \partial_x u = 0,\:C = \frac{a_0 \Delta t}{\Delta x} \]

\[ u_j^{n+1} = u_j^n + \frac{C}{2}(1+C)u_{j-1}^n + (1-C^2)u_{j}^n - \frac{C}{2}(1+C)u_{j+1}^n \]

ODE solvers: computational cost

We want best accuracy and minimum computational cost.

  • Cost related to number of evaluations of \(f(u,t)\), but also to timestep size
  • Harder for adaptive methods, since they change \(\Delta t\) based on error estimate
  • Even harder for implicit methods, since they involve solution of non-linear systems

ODE solvers: accuracy

Example 1: you want to model a system of chemical reactions. Which method would you choose?

  1. Method with \(\mathcal{O}(\Delta t^4)\) error, but still might lead to negative densities
  2. Method with \(\mathcal{O}(\Delta t^2)\) error, but preserves non-negativity of densities

Example 2: you want to model a system of particles. Which method would you choose?

  1. Method with \(\mathcal{O}(\Delta t^4)\) error, but does not conserve energy
  2. Method with \(\mathcal{O}(\Delta t^2)\) error, but conserves energy

We therefore will focus on different (not all!) notions of error/accuracy and talk about some problem-specific methods.

Recap of some standard methods

Let’s recall some textbook methods first.

\[ t_n = n \Delta t, \quad u^n = u(t_n). \]

Tip

Which methods would you name?

Forward Euler

Simplest approach: \[ u^{n+1} = u^n + \Delta t f(u^n, t^n) \]

Properties?

\(\mathcal{O}(\Delta t)\) error; what else?

RK4

\[ k_1 = \Delta t f(u^n, t^n), \]

\[ k_2 = \Delta t f(u^n + \frac{1}{2}k_1, t^n + \frac{1}{2}\Delta t), \]

\[ k_3 = \Delta t f(u^n + \frac{1}{2}k_2, t^n + \frac{1}{2}\Delta t), \]

\[ k_4 = \Delta t f(u^n + k_3, t^n \Delta t), \]

\[ u^{n+1} = u^n + \frac{1}{6}k_1 + \frac{1}{3}k_2 + \frac{1}{3}k_3 + \frac{1}{6}k_4 \]

Properties?

\(\mathcal{O}((\Delta t)^4)\) error; what else?

Stability, monotonicity

Let’s discuss stability and monotonicity. These are not identical properties!

Question

Which one is stronger?

Stability

Linear analysis:

\[ u' = -\alpha u,\: \alpha \in \mathbb{C} > 0 \Rightarrow u(t) = u_0 e^{-\alpha t}. \]

Expect numerically for \(\alpha: Re(\alpha) < 0\):

\[ \lim_{n \to \infty} u^{n} = 0 \]

I.e. solution does not “blow up”.

Stability

For forward Euler: \(\Delta t < \frac{2}{\alpha}\)

Monotonicity

\[ u(t) = u_0 e^{-\alpha t}: t_2 > t_1 \Rightarrow |u(t_1)| > |u(t_2)| \]

Expect numerically: \[ |u^{n}| \leq |u^{n-1}| \]

Monotonicity

For forward Euler: \(\Delta t < \frac{1}{\alpha}\)

Example where monotonicity may be desirable:

\[ \mathbf{u} = (\rho, \rho v_1, \rho v_2, \rho e), \] if \(\rho\) or \(\rho e\) become negative - solver might crash or produce meaningless results!

Generalized stability/monotonicity

Consider \(G(u^{n})\) instead of \(|u^{n}|\), where \(G\) is a 1) norm 2) semi-norm or 3) convex functional

Numerical example

\[ u' = -\alpha u,\: \alpha > 0,\: u(0) = u_0 \]

Note

For forward Euler: Monotonicity: \(\Delta t < \frac{1}{\alpha}\), stability: \(\Delta t < \frac{2}{\alpha}\).

Code 
import numpy as np
from matplotlib import pyplot as plt


def forward_euler(alpha, u0, times):
    u = np.zeros(len(times))
    u[0] = u0
    for (i, t) in enumerate(times[1:]):
        u[i+1] = u[i] + (times[i+1] - times[i]) * (-alpha * u[i])
    return u


fig, ax = plt.subplots(figsize=(6.5, 4.2))

u0 = 2.0

alpha = 10.0
times = np.linspace(0, 1, 100)

times_fe1 = np.linspace(0, 1, 5)
dt1 = times_fe1[1] - times_fe1[0]

times_fe2 = np.linspace(0, 1, 10)
dt2 = times_fe2[1] - times_fe2[0]

times_fe3 = np.linspace(0, 1, 12)
dt3 = times_fe3[1] - times_fe3[0]

dts = r"$\Delta t$"
oas = r"$1/\alpha$"
tas = r"$2/\alpha$"


ax.plot(times, u0 * np.exp(-alpha * times),
        'k-', label=f"Analytical")
ax.plot(times_fe1, forward_euler(alpha, u0, times_fe1),
        color='tab:red', label=f"F.E. ({dts} = {dt1:.3f}; {tas} = {2.0/alpha:.2f})")
ax.plot(times_fe2, forward_euler(alpha, u0, times_fe2),
        color='tab:purple', label=f"F.E. ({dts} = {dt2:.3f}; {oas} = {1.0/alpha:.2f})")
ax.plot(times_fe3, forward_euler(alpha, u0, times_fe3),
        color='tab:blue',
        label=f"F.E. ({dts} = {dt3:.3f})")

ax.legend(fontsize=12)
ax.set_xlabel("t", fontsize=10)
ax.set_ylabel("u", fontsize=10)
ax.tick_params(axis='both', labelsize=8)

ax.set_yscale("log")
plt.show()

Euler’s method for \(du/dt = -a u\)

Region of absolute stability

Re-write discretization of ODE (\(-\alpha \to \lambda\))

\[ u'(t) = \lambda u,\: u^n = R(\lambda \Delta t) u^{n-1}, \]

\(R(z), z \in \mathbb{C}\) is called the stability function.

Region of absolute stability: \[ \mathcal{S} = \{ z \in \mathbb{C} | |R(z)| < 1 \} \]

Region of absolute stability: examples

Forward Euler: \[ R(z) = 1 + z \]

Backward (implicit Euler): \[ u^{n+1} = u^n + \Delta \lambda t u^{n+1} \implies R(z) = \frac{1}{1 + z} \]

RK4: \[ R(z) = 1 + z + \frac{1}{2!} z^2 + \frac{1}{3!} z^3 + \frac{1}{4!} z^4 \]

SSPRK methods

Strong stability-preserving (Total Variation Diminishing) Runge-Kutta methods:

  • ensure strong stability/monotonicity of solutions and have higher order than forward Euler: \(||u^{n+1}|| \leq ||u^n||\) (or a generalized \(G(u)\)) for all \(n\)
  • usually written as convex combination of forward Euler steps

SSPRK methods: formulation

An m-stage method for \(u' = f(u)\) is written as (Shu-Osher form) \[ u^{(0)} = u^n, \]

\[ u^{(i)} = \sum_{j=0}^{i-1} \left(\alpha_{ij} u^{(j)} + \Delta t \beta_{ij} f(u^{(j)}) \right), \quad i = 1, \dots, m; \: \alpha_{ij} \geq 0, \]

\[ u^{n+1} = u^{(m)}. \]

Theorem

Theorem

If \(||u^{n+1}|| \leq ||u^n||\) in the forward Euler method with \(\Delta \leq \Delta_{FE}\), then the SSPRK method is stable for \(\Delta t \leq c_{SSP} \Delta_{FE}\).

  • \(c_{SSP}\) is called the SSP coefficient, \(c_{SSP} = \min \frac{\alpha_{ij}}{\beta_{ij}}\)
  • Most research focused on finding SSPRK methods with largest \(c_{SSP}\)
  • \(\Delta_{FE}\) - depends on spatial discretization, \(c_{SSP}\) - on time-stepping scheme

SSPRK examples

SSPRK2: \[ u^{(1)} = u^n + \Delta t f(u^n), \] \[ u^{(n+1)} = \frac{1}{2}u^n + \frac{1}{2}u^{(1)} + \frac{1}{2}\Delta t f(u^{(1)}). \]

SSPRK3: \[ u^{(1)} = u^n + \Delta t f(u^n), \] \[ u^{(2)} = \frac{3}{4}u^n + \frac{1}{4}u^{(1)} + \frac{1}{4}\Delta t f(u^{(1)}), \] \[ u^{(n+1)} = \frac{1}{3}u^n + \frac{2}{3}u^{(2)} + \frac{2}{3}\Delta t f(u^{(2)}). \]

SSP coefficient, effective CFL

Question

What does \(c_{SSP} = 1\) (for example, SSPRK3) mean?

Makes more sense to talk about effective CFL: \(c_{eff} = c_{SSP} / k\), \(k\) - number of evaluations of RHS. For example, SSPRK3 has \(c_{eff} = 1/3\) (but third-order convergence!).

Stability vs. monotonicity vs. accuracy

SSPRK methods may have same restrictions on \(\Delta t\) to achieve monotonicity, but

  • are more accurate
  • have a larger stability region (monotone \(\neq\) stable!)

Stability for oscillatory problems

Consider harmonic oscillator

\[ \frac{d}{dt} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} = \begin{pmatrix} u_2 \\ - \omega^2 u_1 \end{pmatrix} \]

Analytical solution: \[ \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} = c_1 \begin{pmatrix} \cos(\omega t) \\ -\omega\sin(\omega t) \end{pmatrix} + c_2 \begin{pmatrix} \frac{1}{\omega}\sin(\omega t) \\ \cos(\omega t) \end{pmatrix} \]

Forward Euler

Energy is given by \[ E = (\omega^2 u_1^2 + u_1^2)/2 \] and is conserved for the analytical solution.

Energy blow-up

For forward Euler we can show that \(E^{n+1} = \frac{1}{2} (1 + \omega^2 \Delta t^2) \left( \omega^2 (u_1^n)^1 + (u_2^n)^2 \right)\). Eigenvalues of harmonic oscillator system are purely imaginary (\(\pm i\omega\)) and do not lie in stability region of forward Euler.

Forward Euler: analysis

Let’s write the forward Euler scheme as a matrix equation \[ \begin{pmatrix} u_1^{n+1} \\ u_2^{n+1} \end{pmatrix} = \begin{pmatrix} 1 & \Delta t \\ -\omega^2 \Delta t & 1 \end{pmatrix} \begin{pmatrix} u_1^{n} \\ u_2^{n} \end{pmatrix} =: J \mathbf{u} \]

Determinant of \(J\)

Note that \(det(J) = (1 + \omega^2 \Delta t^2)\) is exactly the energy growth factor!

Re-write integration step via growth factor \(\gamma\): \[ u_1^{n+1} = \gamma u_1^n,\:u_2^{n+1} = \gamma u_2^n \implies \begin{pmatrix} \gamma - 1 & -\Delta t \\ \omega^2 \Delta t & \gamma - 1 \end{pmatrix} \begin{pmatrix} u_1^{n} \\ u_2^{n} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \]

Exercise

Non-trivial solution \(\implies\) \(det(LHS) = 0\); can show that \(|\gamma| > 1\), i.e. method is unstable.

Example - harmonic oscillator

Code 
def forward_euler_generic(rhs, u0, times, omega):
    u = np.zeros((len(times), len(u0)))
    u[0,:] = u0
    for (i, t) in enumerate(times[1:]):
        u[i+1,:] = u[i,:] + (times[i+1] - times[i]) * rhs(u[i,:], omega)
    return u

def ssprk22_generic(rhs, u0, times, omega):
    u = np.zeros((len(times), len(u0)))
    u[0,:] = u0
    for (i, t) in enumerate(times[1:]):
        u_1 = u[i,:] + (times[i+1] - times[i]) * rhs(u[i,:], omega)
        u[i+1,:] = 0.5 * u[i,:]+ 0.5 * u_1 + 0.5 * (times[i+1] - times[i]) * rhs(u_1, omega)
    return u

def rk2_generic(rhs, u0, times, omega):
    # Ralston's method
    u = np.zeros((len(times), len(u0)))
    u[0,:] = u0
    for (i, t) in enumerate(times[1:]):
        dt = (times[i+1] - times[i])
        k1 = rhs(u[i,:], omega)
        k2 = rhs(u[i,:] + 0.75 * k1 * dt, omega)
        u[i+1,:] = u[i,:] +  dt * (0.25 * k1 + 0.75 * k2)
    return u

def harmonic_rhs(u, omega):
    return np.array([u[1], -omega**2 * u[0]])

def energy(u, omega):
    return 0.5 * (u[:,1]**2 + omega**2 * u[:,0]**2)

fig, axes = plt.subplots(figsize=(12.0, 4.2), ncols=2)

ax = axes[0]
ax_e = axes[1]

t_max = 5.0
omega = 3.0
times = np.linspace(0, t_max, 100)

u0 = [0.5, 0.5]

analytical = np.array([u0[0] * np.cos(omega * times) + (u0[0]/omega) * np.sin(omega * times),
                       -u0[0] * omega * np.sin(omega * times) + u0[0] * np.cos(omega * times)]).T



times_fe1 = np.linspace(0, t_max, 100)
dt1 = times_fe1[1] - times_fe1[0]
sol_fe1 = forward_euler_generic(harmonic_rhs, u0, times_fe1, omega)

times_fe2 = np.linspace(0, t_max, 400)
dt2 = times_fe2[1] - times_fe2[0]
sol_fe2 = forward_euler_generic(harmonic_rhs, u0, times_fe2, omega)

times_ssprk2 = np.linspace(0, t_max, 80)
dt3 = times_ssprk2[1] - times_ssprk2[0]
sol_ssprk2 = ssprk22_generic(harmonic_rhs, u0, times_ssprk2, omega)

times_rk2 = np.linspace(0, t_max, 80)
dt4 = times_rk2[1] - times_rk2[0]
sol_rk2 = rk2_generic(harmonic_rhs, u0, times_rk2, omega)


ax.plot(sol_fe1[:,0], sol_fe1[:,1],
        color='tab:red', label=f"F.E. dt = {dt1:.3f}")

ax.plot(sol_fe2[:,0], sol_fe2[:,1],
        color='tab:purple', label=f"F.E. dt = {dt2:.3f}")

ax.plot(sol_ssprk2[:,0], sol_ssprk2[:,1],
        color='tab:blue', label=f"SSPRK2 dt = {dt3:.3f}")

ax.plot(sol_rk2[:,0], sol_rk2[:,1],
        color='tab:green', label=f"RK2 dt = {dt4:.3f}")


ax.plot(analytical[:,0], analytical[:,1],
        'k-', label=f"Analytical")

ax_e.plot(times, energy(analytical, omega), 'k')
ax_e.plot(times_fe1, energy(sol_fe1, omega), 'tab:red')
ax_e.plot(times_fe2, energy(sol_fe2, omega), 'tab:purple')
ax_e.plot(times_ssprk2, energy(sol_ssprk2, omega), 'tab:blue')
ax_e.plot(times_rk2, energy(sol_rk2, omega), 'tab:green')

ax.legend(fontsize=12)
ax.set_xlabel("x(t)", fontsize=10)
ax.set_ylabel("v(t)", fontsize=10)
ax.tick_params(axis='both', labelsize=8)


print(f"dE F.E. (dt={dt1:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_fe1, omega)[-1]):.3e}")
print(f"dE F.E. (dt={dt2:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_fe2, omega)[-1]):.3e}")
print(f"dE SSPRK2 (dt={dt3:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_ssprk2, omega)[-1]):.3e}")
print(f"dE RK2 (dt={dt4:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_rk2, omega)[-1]):.3e}")

ax_e.set_yscale("log")
plt.show()
dE F.E. (dt=0.051) at t=5.0: -1.057e+01
dE F.E. (dt=0.013) at t=5.0: -9.460e-01
dE SSPRK2 (dt=0.063) at t=5.0: -3.250e-02
dE RK2 (dt=0.063) at t=5.0: 3.461e-01

Harmonic oscillator phase space (left) and energy (right), omega=3

Question

Why does SSPRK2 still cause a growth of energy?

Stiffness

We don’t go into detail about stiff ODEs in this lecture, but still mention some definitions.

  • Definition 1 Step size required for stability is much smaller than the step size required for accuracy
  • Definition 2 Computational cost of an explicit method becomes larger than that of an implicit one

Arise in biology, chemistry, low-Mach number compressible flows.

Symplectic integrators

Consider ODE system given by

\[ \frac{d q_i}{d t} = -\frac{\partial H}{\partial p_i}, \quad \frac{d p_i}{d t} = \frac{\partial H}{\partial q_i}, \] where \(H\) is the Hamiltonian of the system.

  • Why do are we interested in specific methods for such systems?
  • What can these systems describe?

Naming

\(q_i\) are called “generalized coordinates”, \(p_i\) - “generalized momenta”

Example: interacting particles

\[ H(p,q) = \frac{1}{2} \frac{1}{m}\sum_i p_i^2 + \sum_{i<j} V(|q_i - q_j|), \]

\(V\) - interaction potential. Examples: gravity, electrostatic potential (both \(V(r)\) are proportional to \(1/r\)). Potential gives force:

\[ F = -\nabla V. \]

Hamiltonian of interacting particles?

What does Hamiltonian resemble?

Hamiltonian is kinetic + potential energy! Time integration of particle movement should conserve total energy.

First integrals

Cosider the system \[ \frac{d y}{dt} = f(y). \]

\(I(y)\) is called a first integral of a system, if

\[ \frac{d I}{dt} = \frac{d I}{d y} f(y) = 0. \]

H as first integral

The Hamiltonian is a first integral of a Hamiltonian system!

Symplectic integrators

Symplectic integrators preserve Hamiltonian structure. Energy conservation by default!

  • Exact solution to inexact (but Hamiltonian) approximation of Hamiltonian system

Define mapping \(\psi\):

\[ \psi: (\mathbf{q}(t), \mathbf{p}(t)) \to (\mathbf{q}(t+ \Delta t \mathbf{p}), \mathbf{p}(t + \Delta t)) \]

Theorem Symplectic iff.

\[ (D\psi)^T \begin{pmatrix} 0 & I \\ -I & 0 \end{pmatrix} (D\psi) = \begin{pmatrix} 0 & I \\ -I & 0 \end{pmatrix}, \] \(D \psi\) is the Jacobi matrix of the mapping (\(\psi\) - “canonical transformation”.).

Phase space volume-preserving integrators

A weaker condition: \(\mu(\psi(\mathcal{S})) \equiv \mu(\mathcal{S})\), where \(\mu\) is a measure (volume). This implies that the Jacobi matrix \(D\psi\) has determinant 1.

What does this give?

  • Ensures no spurious sources/sinks/attractors
  • Ensures that time-averaging remains equivalent to ensemble averaging for a certain class of systems (ergodic systems): important for particle methods!

Difference between symplectic and volume-preserving mappings:

The non-squeezing theorem explains the difference nicely and shows that symplecticity is a stronger condition.

Leapfrog integrator

If \(H(\mathbf{p},\mathbf{q}) = T(\mathbf{p}) + V(\mathbf{q})\) (separable Hamiltonian), then we can integrate with leapfrog:

\[ \mathbf{p}^{n+1/2} = \mathbf{p}^{n} - \Delta t\nabla{ V}(\mathbf{q}^n), \] \[ \mathbf{q}^{n+1} = \mathbf{q}^{n} + \Delta t\nabla{ T}(\mathbf{p}^{n+1/2}), \]

Leapfrog is symplectic and has error \(\mathcal{O}((\Delta t)^2)\).

Starting leapfrog?

Leapfrog requires time offset between coordinates and momenta (it is essentially integrating using the midpoint rule). How do we start a simulation if we know only \(\mathbf{q}(t=0)\) and \(\mathbf{p}(t=0)\)?

Harmonic oscillator

One can show that stability for harmonic oscillator requires \(\Delta t \leq 2/\omega\); this can be done using the growth factor analysis (keeping in mind that a step of \(\Delta t/2\) corresponds to a growth factor of \(\sqrt{\gamma}\) and re-writing \(u_2^{n+1/2}\) and \(u_2^{n-1/2}\) in terms of \(u_2^n\)).

Verlet integrator

What if we need \(\mathbf{q}\), \(\mathbf{p}\) at the same point in time?

\[ \mathbf{p}^{n+1/2} = \mathbf{p}^{n} - \frac{\Delta t}{2}\nabla{V}(\mathbf{q}^n), \] \[ \mathbf{q}^{n+1} = \mathbf{q}^{n} + \Delta t\nabla{T}(\mathbf{p}^{n+1/2}), \] \[ \mathbf{p}^{n+1} = \mathbf{p}^{n+1/2} - \frac{\Delta t}{2}\nabla{V}(\mathbf{q}^{n+1}), \]

Naming

Often this is called “Leapfrog in kick-drift-kick” form or “Störmer-Verlet”.

Example - harmonic oscillator revisited

Code 
def verlet_generic(rhs, u0, times, omega):
    u = np.zeros((len(times), len(u0)))
    u[0,:] = u0
    for (i, t) in enumerate(times[1:]):
        dt = (times[i+1] - times[i])
        p_hs = u[i,1] + 0.5 * dt * rhs(u[i,:], omega)[1]
        u[i+1,0] = u[i,0] + dt * p_hs
        u[i+1,1] = p_hs + 0.5 * dt * rhs(u[i+1,:], omega)[1]
    return u

fig, axes = plt.subplots(figsize=(12.0, 4.2), ncols=2)

ax = axes[0]
ax_e = axes[1]

times_verlet = np.linspace(0, t_max, 20)
dt5 = times_verlet[1] - times_verlet[0]
sol_verlet = verlet_generic(harmonic_rhs, u0, times_verlet, omega)

ax.plot(sol_fe2[:,0], sol_fe2[:,1],
        color='tab:purple', label=f"F.E. dt = {dt2:.3f}")

ax.plot(sol_ssprk2[:,0], sol_ssprk2[:,1],
        color='tab:blue', label=f"SSPRK2 dt = {dt3:.3f}")

ax.plot(sol_verlet[:,0], sol_verlet[:,1],
        color='tab:red', label=f"Verlet dt = {dt5:.3f}; 2/omega={2/omega:.3f}")

ax.plot(analytical[:,0], analytical[:,1],
        'k-', label=f"Analytical")

ax_e.plot(times, energy(analytical, omega), 'k')
ax_e.plot(times_fe2, energy(sol_fe2, omega), 'tab:purple')
ax_e.plot(times_ssprk2, energy(sol_ssprk2, omega), 'tab:blue')
ax_e.plot(times_verlet, energy(sol_verlet, omega), 'tab:red')

ax.legend(fontsize=12, loc="upper right")
ax.set_xlabel("x(t)", fontsize=10)
ax.set_ylabel("v(t)", fontsize=10)
ax.tick_params(axis='both', labelsize=8)


print(f"dE F.E. (dt={dt2:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_fe2, omega)[-1]):.3e}")
print(f"dE SSPRK2 (dt={dt3:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_ssprk2, omega)[-1]):.3e}")
print(f"dE Verlet (dt={dt5:.3f}) at t={t_max}: {(energy(analytical, omega)[-1] - energy(sol_verlet, omega)[-1]):.3e}")

# ax_e.set_yscale("log")
plt.show()
dE F.E. (dt=0.013) at t=5.0: -9.460e-01
dE SSPRK2 (dt=0.063) at t=5.0: -3.250e-02
dE Verlet (dt=0.263) at t=5.0: 4.697e-02

Harmonic oscillator phase space (left) and energy (right), omega=3

Code

Summary

Methods not covered in this lecture but also useful/actively developed:

Exercises

  1. Derive the stability polynomial \(R(z)\) for the RK4 method
  2. Prove the SSP property of a method in Shu-Osher form (hint: use the TVD property of forward Euler, the triangle inequality and look at the definition of \(c_{SSP}\))
  3. Derive energy growth rate for SSPRK2 for the harmonic oscillator
  4. Show that forward Euler is always unstable for the harmonic oscillator via growth factor analysis
  5. Show that the forward Euler method is not a symplectic method

Questions

  1. What is the difference between monotonicity and stability of a method?
  2. Why are implicit methods oftentimes more computationally expensive?
  3. What is the difference between a symplectic and a phase space volume-preserving integrator? Which condition is stronger? Why?
  4. What two formulations of explicit symplectic integrators do you know and what is the difference between them?

References

  1. S. Gottlieb, On High Order Strong Stability Preserving Runge–Kutta and Multi Step Time Discretizations
  2. Kubatko, Yeager et al., Optimal Strong-Stability-Preserving Runge–Kutta Time Discretizations for Discontinuous Galerkin Methods
  3. A. Giorgilli, Notes on Hamiltonian Dynamical Systems