Runge–Kutta Schemes

Runge-Kutta Schemes

rk_symbolic_weight(t: Tree | Forest | ForestSum, s: int, explicit: bool = False, a_mask: list = None, b_mask: list = None, mathematica_code: bool = False, rationalise: bool = True) → sympy.core.add.Add | str

Returns the elementary weight of a Tree, Forest or ForestSum \(t\) as a SymPy symbolic expression.

Parameters:

• t – A Tree, Forest or ForestSum

• s (int) – The number of Runge–Kutta stages

• explicit (bool) – If true, assumes the Runge–Kutta scheme is explicit, i.e. \(a_{ij} = 0\) for \(i \leq j\).

• a_mask – A two-dimensional array specifying which coefficients of the Runge–Kutta parameter matrix \(A\) are non-zero. If not None, sets \(a_{ij} = 0\) for all \(i,j\) such that a_mask[i][j] = 0.

• b_mask – A one-dimensional array or list specifying which coefficients of the Runge–Kutta parameter vector \(b\) are non-zero. If not None, sets \(b_i = 0\) for all \(i\) such that b_mask[i] = 0.

• mathematica_code (bool) – If true, outputs the expression as mathematica code.

• rationalise (bool) – If true, will attempt to rationalise the coefficients in the expression

Returns:

The elementary weight of \(t\), as a SymPy symbolic expression if mathematica_code is False or as a string if mathematica_code is True.

Return type:

sympy.core.add.Add | string

Example usage:

rk_order_cond(t: Tree | Forest | ForestSum, s: int, explicit: bool = False, a_mask: list = None, b_mask: list = None, mathematica_code: bool = False, rationalise: bool = True) → sympy.core.add.Add | str

Returns the Runge–Kutta order condition associated with tree \(t\) as a SymPy symbolic expression.

Parameters:

• t – A Tree

• s (int) – The number of Runge–Kutta stages

• explicit (bool) – If true, assumes the Runge–Kutta scheme is explicit, i.e. \(a_{ij} = 0\) for \(i \leq j\).

• a_mask – A two-dimensional array specifying which coefficients of the Runge–Kutta parameter matrix \(A\) are non-zero. If not None, sets \(a_{ij} = 0\) for all \(i,j\) such that a_mask[i][j] = 0.

• b_mask – A one-dimensional array or list specifying which coefficients of the Runge–Kutta parameter vector \(b\) are non-zero. If not None, sets \(b_i = 0\) for all \(i\) such that b_mask[i] = 0.

• mathematica_code (bool) – If true, outputs the expression as mathematica code.

• rationalise (bool) – If true, will attempt to rationalise the coefficients in the expression

Returns:

The order condition associated with the tree \(t\), as a SymPy symbolic expression if mathematica_code is False or as a string if mathematica_code is True.

Return type:

sympy.core.add.Add | string

Example usage:

class RK(a, b, name=None)

A Runge–Kutta method with the Butcher tableau:

\[\begin{split}\begin{array}{c|c} c & A \\ \hline & b^T \end{array}\end{split}\]

where \(c_i = \sum_{j=1}^s a_{ij}\).

Parameters:

• a – The Runge–Kutta parameter matrix \(A\).

• b – The Runge–Kutta parameter vector \(b\).

__mul__(other: RK) → RK

Returns the composition of two RK schemes, \((A_1, b_1)\) and \((A_2, b_2)\), with Butcher tableau:

\[\begin{split}\begin{array}{c|cc} c_1 & A_1 & 0 \\ c_2 & e b_1^T & A_2\\ \hline & b_1 & b_2 \end{array}\end{split}\]

where \(e\) is the vector of 1s.

Return type:

RK

__pow__(exponent: int) → RK

Returns the compositional power of the Runge–Kutta scheme. In particular, self ** (-1) returns the scheme with Butcher tableau:

\[\begin{split}\begin{array}{c|c} & A - e b^T \\ \hline & -b^T \end{array}\end{split}\]

where \(e\) is the vector of 1s.

Parameters:

exponent (int) – Exponent

Return type:

RK

adjoint() → RK

Returns the adjoint Runge–Kutta method, given by the Butcher tableau:

\[\begin{split}\begin{array}{c|c} \widetilde{c} & e \widetilde{b}^T - \widetilde{A} \\ \hline & \widetilde{b}^T \end{array}\end{split}\]

where \(\widetilde{b}_i := b_{s+1-i}\) and \(\widetilde{A}_{ij} := A_{s+1 - i, s+ 1 - j}\) for all \(1 \leq i,j \leq s\).

Return type:

RK

antisymmetric_order(tol: float = 1e-10, limit: int = 10) → int

Returns the antisymmetric order of the RK scheme. See [SEFTS25] for details.

Parameters:

• tol (float) – Tolerance for evaluating order conditions.

• limit (int) – Highest admissible order. If the order equals or exceeds this limit, a runtime error will be raised.

Return type:

int

elementary_weights_map() → Map

Returns the elementary weight function of the Runge-Kutta method as an instance of the Map class.

Return type:

Map

modified_equation_map() → Map

Returns the map corresponding to the elementary weights function of the modified (B-series) vector field, \(\widetilde{\phi}\), defined by

\[(\widetilde{\phi} \star e)(t) = \phi(t)\]

where \(\phi\) is the elementary weights function of the Runge-Kutta scheme and \(e(t) = 1 / t!\) is the elementary weights function of the exact solution. Equivalently,

\[\widetilde{\phi}(t) = (\phi \star e^{\star (-1)})(t).\]

Returns:

Elementary weights map of the modified vector field

Return type:

Map

order(tol: float = 1e-10, limit: int = 10) → int

Returns the order of the RK scheme.

Parameters:

• tol (float) – Tolerance for evaluating order conditions. An order condition of the form self.elementary_weights(t) = 1./t.factorial() is considered to be satisfied if abs( self.elementary_weights(t) - 1./t.factorial() ) < tol

• limit (int) – Highest admissible order. If the order equals or exceeds this limit, a runtime error will be raised.

Return type:

int

planar_antisymmetric_order(tol: float = 1e-10, limit: int = 10) → int

Returns the antisymmetric order of the RK scheme on ordered (planar) trees, using the NCK Hopf algebra.

Checks D(tau) = ((sign . Phi) *_nck Phi)(tau) - epsilon(tau) = 0 for all ordered trees tau.

Parameters:

• tol (float) – Tolerance for evaluating order conditions.

• limit (int) – Highest admissible order.

Return type:

int

planar_order(tol: float = 1e-10, limit: int = 10) → int

Returns the order of the RK scheme on ordered (planar) trees.

Checks Phi(tau) = 1/gamma(tau) for all ordered trees tau.

Parameters:

• tol (float) – Tolerance for evaluating order conditions.

• limit (int) – Highest admissible order.

Return type:

int

reverse() → RK

Returns the RK scheme given by reversing the step size h to -h, with Butcher tableau:

\[\begin{split}\begin{array}{c|c} -c & -A \\ \hline & -b^T \end{array}\end{split}\]

Return type:

RK

run(y0: list | np.ndarray, t0: float, t_end: float, f: Callable[[float, float], list | np.ndarray], n: int, tol: float = 1e-10, max_iter: int = 100, plot: bool = False, plot_dims: list | np.ndarray = None, plot_kwargs: dict = None) → Tuple[list, list]

Runs the Runge–Kutta method.

Parameters:

• y0 (list | array) – Initial condition for y

• t0 (float) – Initial condition for t

• t_end (float) – End point for t

• f (callable) – Function defining the ODE \(dy / dt = f(t,y)\).

• n (int) – Number of steps

• tol (float) – Tolerance for the root solving algorithm. Only applicable if the scheme is implicit.

• max_iter (int) – Maximum number of iterations for the root solving algorithm. Only applicable if the scheme is implicit.

• plot (bool) – If true, will plot the solution

• plot_dims (list | array) – List of dimensions of the solution to plot

• plot_kwargs (dict) – kwargs to pass to pyplot.plot() if plotting the solution.

Returns:

t_vals, y_vals - the lists of values of t and y respectively

Return type:

tuple[list, list]

step(y0: list | np.ndarray, t0: float, f: Callable[[float, float], list | np.ndarray], h: float, tol: float = 1e-10, max_iter: int = 100) → list | np.ndarray

Runs one step of the Runge–Kutta method.

Parameters:

• y0 (list | array) – Initial condition for y

• t0 (float) – Initial condition for t

• f (callable) – Function defining the ODE \(dy / dt = f(t,y)\).

• h (float) – Step size

• tol (float) – Tolerance for the root solving algorithm. Only applicable if the scheme is implicit.

• max_iter (int) – Maximum number of iterations for the root solving algorithm. Only applicable if the scheme is implicit.

Returns:

Next point, y1

Return type:

list | array

Instances

We provide a handful of Runge-Kutta schemes for convenience, as instances of the RK class. We list these below, along with their Butcher tableux.

Explicit Methods

euler

The Euler method

\[\begin{split}\begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array}\end{split}\]

heun_rk2

Heun’s RK2 method

\[\begin{split}\begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array}\end{split}\]

midpoint

The midpoint method

\[\begin{split}\begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \end{array}\end{split}\]

kutta_rk3

Kutta’s RK3 method

\[\begin{split}\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array}\end{split}\]

heun_rk3

Heun’s RK3 method

\[\begin{split}\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 2/3 & 0 & 2/3 & 0 \\ \hline & 1/4 & 0 & 3/4 \end{array}\end{split}\]

ralston_rk3

Ralston’s RK3 method

\[\begin{split}\begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 3/4 & 0 & 3/4 & 0 \\ \hline & 2/9 & 1/3 & 4/9 \end{array}\end{split}\]

rk4

The RK4 method

\[\begin{split}\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & 1/6 & 1/3 & 1/3 & 1/6 \end{array}\end{split}\]

ralston_rk4

Ralston’s RK4 method

\[\begin{split}\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ \frac{2}{5} & \frac{2}{5} & 0 & 0 & 0 \\ \frac{14 - 3\sqrt{5}}{16} & \frac{-2889 + 1428\sqrt{5}}{1024} & \frac{3785 - 1620\sqrt{5}}{1024} & 0 & 0 \\ 1 & \frac{-3365 + 2094\sqrt{5}}{6040} & \frac{-975-3046\sqrt{5}}{2552} & \frac{467040+203968\sqrt{5}}{240845} & 0 \\ \hline & \frac{263+24\sqrt{5}}{1812} & \frac{125 - 1000\sqrt{5}}{3828} & \frac{3426304+1661952\sqrt{5}}{5924787} & \frac{30-4\sqrt{5}}{123} \end{array}\end{split}\]

nystrom_rk5

Nyström’s RK5 method

\[\begin{split}\begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 & 0 & 0 & 0 \\ 2/5 & 4/25 & 6/25 & 0 & 0 & 0 & 0 \\ 1 & 1/4 & -3 & 15/4 & 0 & 0 & 0 \\ 2/3 & 2/27 & 10/9 & -50/81 & 8/81 & 0 & 0 \\ 4/5 & 2/25 & 12/25 & 2/15 & 8/75 & 0 & 0 \\ \hline & 23/192 & 0 & 125/192 & 0 & -27/64 & 125/192 \end{array}\end{split}\]

EES25(x)

The Explicit and Effectively Symmetric scheme of order 2 and antisymmetric order 5, \(EES(2,5;x)\) [SEFTS25].

\[\begin{split}\begin{array}{c|cccc} 0 & 0 & 0 & 0 \\ \displaystyle\frac{1+2x}{4(1-x)} & \displaystyle\frac{1+2x}{4(1-x)} & 0 & 0 \\ \displaystyle\frac{3}{4(1-x)} & \displaystyle\frac{(4x-1)^2}{4(x-1)(1-4x^2)} & \displaystyle\frac{1-x}{(1-4x^2)} & 0 \\ \hline & x & \displaystyle\frac{1}{2} & \displaystyle\frac{1}{2} - x \end{array}\end{split}\]

Parameters:

x – Parameter of the \(EES(2,5)\) scheme.

Return type:

kauri.RK

EES27(x, plus=True)

The Explicit and Effectively Symmetric scheme of order 2 and antisymmetric order 7, \(EES(2,7;x)\). For the Butcher tableau, see [SEFTS25].

Parameters:

• x – Parameter of the \(EES(2,7)\) scheme.

• plus – If True, takes \(+\sqrt{2}\) in the Butcher tableau, otherwise \(-\sqrt{2}\).

Return type:

kauri.RK

Implicit Methods

backward_euler

The backward Euler method

\[\begin{split}\begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array}\end{split}\]

implicit_midpoint

The implicit midpoint method

\[\begin{split}\begin{array}{c|c} 1/2 & 1/2 \\ \hline & 1 \end{array}\end{split}\]

crank_nicolson

The Crank Nicolson method

\[\begin{split}\begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1/2 & 1/2 \\ \hline & 1/2 & 1/2 \end{array}\end{split}\]

gauss6

Kuntzmann & Butcher method of order 6, based on Gaussian quadrature

\[\begin{split}\begin{array}{c|ccc} \frac{1}{2} - \frac{\sqrt{15}}{10} & \frac{5}{36} & \frac{2}{9} - \frac{\sqrt{15}}{15} & \frac{5}{36} - \frac{\sqrt{15}}{30} \\ \frac{1}{2} & \frac{5}{36} + \frac{\sqrt{15}}{24} & \frac{2}{9} & \frac{5}{36} - \frac{\sqrt{15}}{24} \\ \frac{1}{2} + \frac{\sqrt{15}}{10} & \frac{5}{36} + \frac{\sqrt{15}}{30} & \frac{2}{9} + \frac{\sqrt{15}}{15} & \frac{5}{36} \\ \hline & \frac{5}{18} & \frac{4}{9} & \frac{5}{18} \end{array}\end{split}\]

radau_iia

The Radau IIA method of order 5

\[\begin{split}\begin{array}{c|ccc} \frac{4 - \sqrt{6}}{10} & \frac{88 - 7\sqrt{6}}{360} & \frac{296 - 169\sqrt{6}}{1800} & \frac{-2 + 3\sqrt{6}}{225} \\ \frac{4 + \sqrt{6}}{10} & \frac{296 + 169\sqrt{6}}{1800} & \frac{88 + 7\sqrt{6}}{360} & \frac{-2 - 3\sqrt{6}}{225} \\ 1 & \frac{16 - \sqrt{6}}{36} & \frac{16 + \sqrt{6}}{36} & \frac{1}{9} \\ \hline & \frac{16 - \sqrt{6}}{36} & \frac{16 + \sqrt{6}}{36} & \frac{1}{9} \\ \end{array}\end{split}\]

lobatto6

Butcher’s Lobatto formula of order 6.

Note

This uses an equivalent but non-standard Butcher tableau where the last column of A is zero. The elementary weights and numerical solution are identical to the standard Lobatto IIIA(4) tableau, but the A matrix does not satisfy the standard Lobatto IIIA property \(a_{s,j} = b_j\).

\[\begin{split}\begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ \frac{5 - \sqrt{5}}{10} & \frac{5 + \sqrt{5}}{60} & \frac{1}{6} & \frac{15 - 7\sqrt{5}}{60} & 0 \\ \frac{5 + \sqrt{5}}{10} & \frac{5 - \sqrt{5}}{60} & \frac{15 + 7\sqrt{5}}{60} & \frac{1}{6} & 0 \\ 1 & \frac{1}{6} & \frac{5 - \sqrt{5}}{12} & \frac{5 + \sqrt{5}}{12} & 0 \\ \hline & \frac{1}{12} & \frac{5}{12} & \frac{5}{12} & \frac{1}{12} \end{array}\end{split}\]