Integration methods¶

Generic time-stepping routine¶

All integration methods are used via the same time-stepping routine:

example_gni_project.simulate(ivp, scheme, h, duration)¶

Compute a discrete trajectory for an initial-value problem.

Parameters:

ivp (IVP) – initial value problem
scheme (function) – function defining the integration scheme
h (float) – time step size
duration (float) – duration of simulation

Returns:

simulation – simulation result

Return type:

SimulationResult

Example

Compute the evolution of a given initial value problem ivp for 10 seconds using the explicit Euler method with a time step of 0.1 seconds.

>>> import autograd.numpy as np
>>> from example_gni_project import *
>>> def vector_field(x):
...     q, p = x
...     return np.array([p, -q])
>>> ivp = IVP(("q", "p"), vector_field, np.array([1., 0.]))
>>> simulation = simulate(ivp, explicit_euler, 0.1, 0.5)
>>> simulation.trajectory[-1, :]
array([ 0.9005 , -0.49001])

Non-symplectic methods¶

Below \(X \colon \mathcal{X} \rightarrow T \mathcal{X}\) denotes the vector field and \(x_0 \in \mathcal{X}\) denotes the initial state. Together they define the initial value problem. Further, \(x_1 \in \mathcal{X}\) denotes the state of the system after one time step of size \(h\).

Explicit Euler¶

The explicit Euler method is the simplest approximation one can make:

\[x_1 = x_0 + h \, X(x_0)\]

The method is first-order and unstable in the sense that the integrator will predict a gain of energy due to numerical errors.

Implicit Euler¶

The implicit Euler method is already significantly more complicated in its implementation since it requires solving a system of possibly nonlinear equations:

\[x_1 = x_0 + h \, X(x_1)\]

The method is first-order and stable in the sense that the integrator will predict a loss of energy due to numerical errors.

Symplectic methods¶

Symplectic methods have good long-time energy behavior. They (almost) conserve the energy (Hamiltonian) of a Hamiltonian system. More precisely, symplectic methods do not conserve the physical energy but they conserve a nearby/perturbed Hamiltonian, see e.g. Reich 1999 and references therein.

Since symplectic methods mimic a (canonical) Hamiltonian structure at the discrete level, the state is partitioned into configuration variables \(q\) and conjugate momentum variables \(p\), i.e. \(x = \left[ q, \, p \right]\).

Symplectic Euler¶

The symplectic Euler method essentially combines the explicit Euler method and the implicit Euler method:

\[\begin{split}\begin{bmatrix} q_1 \\ p_1 \end{bmatrix} = \begin{bmatrix} q_0 \\ p_0 \end{bmatrix} + h \, X \biggl( \begin{bmatrix} q_0 \\ p_1 \end{bmatrix} \biggr)\end{split}\]

Another symplectic Euler method results when combining the implicit Euler method and the explicit Euler method in the other obvious way:

\[\begin{split}\begin{bmatrix} q_1 \\ p_1 \end{bmatrix} = \begin{bmatrix} q_0 \\ p_0 \end{bmatrix} + h \, X \biggl( \begin{bmatrix} q_1 \\ p_0 \end{bmatrix} \biggr)\end{split}\]

Störmer-Verlet¶

The Störmer-Verlet method results from the composition of the first and the second symplectic Euler method, which is asserted by the test_composition_störmer_verlet function in tests/test_integrators.py. To concretely describe the Störmer-Verlet method, the vector field is partitioned analogously to the state i.e. \(x = \left[ q, \, p \right]\) and \(X = \left[ X_q, \, X_p \right] = \left[ +\frac{\partial H(q, p)}{\partial p}, -\frac{\partial H(q, p)}{\partial q} \right]\). As long as \(H\) is separable, \(X_q\) only depends on \(p\) while \(X_p\) only depends on \(q\).

\begin{align} p_{1/2} &= p_0 + \frac{h}{2} \, X_p(q_0) \\ q_1 &= q_0 + h \, p_{1/2} \\ p_1 &= p_{1/2} + \frac{h}{2} \, X_p(q_1) \end{align}

The Störmer-Verlet method is the composition of the two symplectic Euler methods. As such, it is a symmetric second-order method. The symmetry is asserted by test_symmetry_störmer_verlet in test/test_integrators.py.