Ruth4Integrator#

class gala.integrate.Ruth4Integrator(func, func_args=(), func_units=None, progress=False, store_all=True)[source]#

Bases: Integrator

A 4th order symplectic integrator.

Given a function for computing time derivatives of the phase-space coordinates, this object computes the orbit at specified times.

See also

https://en.wikipedia.org/wiki/Symplectic_integrator#A_fourth-order_example

Naming convention for variables:

im1 = i-1
im1_2 = i-1/2
ip1 = i+1
ip1_2 = i+1/2

Parameters:

funcfunc: A callable object that computes the phase-space time derivatives at a time and point in phase space.
func_argstuple (optional): Any extra arguments for the derivative function.
func_unitsUnitSystem (optional): If using units, this is the unit system assumed by the integrand function.

Examples

Using q as our coordinate variable and p as the conjugate momentum, we want to numerically solve for an orbit in the potential (Hamiltonian)

\[\begin{split}\Phi &= \frac{1}{2}q^2\\ H(q, p) &= \frac{1}{2}(p^2 + q^2)\end{split}\]

In this system,

\[\begin{split}\dot{q} &= \frac{\partial \Phi}{\partial p} = p \\ \dot{p} &= -\frac{\partial \Phi}{\partial q} = -q\end{split}\]

We will use the variable w to represent the full phase-space vector, \(w = (q, p)\). We define a function that computes the time derivates at any given time, t, and phase-space position, w:

def F(t, w):
    dw = [w[1], -w[0]]
    return dw

Note

The force here is not time dependent, but this function always has to accept the independent variable (e.g., time) as the first argument.

To create an integrator object, just pass this acceleration function in to the constructor, and then we can integrate orbits from a given vector of initial conditions:

integrator = Ruth4Integrator(acceleration)
times, ws = integrator(w0=[1., 0.], dt=0.1, n_steps=1000)

Note

When integrating a single vector of initial conditions, the return array will have 2 axes. In the above example, the returned array will have shape (2, 1001). If an array of initial conditions are passed in, the return array will have 3 axes, where the last axis is for the individual orbits.

Methods Summary

`__call__`(w0[, mmap])	Run the integrator starting from the specified phase-space position.
`run`(w0[, mmap])
`step`(t, w, dt)	Step forward the positions and velocities by the given timestep.

Methods Documentation

__call__(w0, mmap=None, **time_spec)[source]#

Run the integrator starting from the specified phase-space position. The initial conditions w0 should be a PhaseSpacePosition instance.

There are a few combinations of keyword arguments accepted for specifying the timestepping. For example, you can specify a fixed timestep (dt) and a number of steps (n_steps), or an array of times:

dt, n_steps[, t1] : (numeric, int[, numeric])
    A fixed timestep dt and a number of steps to run for.
dt, t1, t2 : (numeric, numeric, numeric)
    A fixed timestep dt, an initial time, and a final time.
t : array-like
    An array of times to solve on.

Parameters:

w0PhaseSpacePosition: Initial conditions.
**time_spec: Timestep information passed to parse_time_specification.

Returns:

orbitOrbit

run(w0, mmap=None, **time_spec)#: Deprecated since version 1.9: The run function is deprecated and may be removed in a future version. Use Integrator call method instead.

Run the integrator starting from the specified phase-space position.

Deprecated since version 1.9: Use the __call__ method instead.

step(t, w, dt)[source]#

Step forward the positions and velocities by the given timestep.

Parameters:

dtnumeric: The timestep to move forward.