solprop.stepper.tableau#

Explicit Runge-Kutta tableaux with embedded error estimates.

Two order-5(4) pairs are provided: DOPRI5 [1] (the default scheme) and TSIT5 [2]. Both are first-same-as-last (FSAL), so a step costs six evaluations of the nonlinear term and the embedded companion supplies the local error estimate for free (docs/dev/stepper.md sections 2 and 3).

A tableau is data, not behaviour: it carries the nodes, the Butcher matrix, the propagated weights and the error weights, and it validates its own structure. check_order_conditions() measures the order the coefficients actually attain, from the rooted-tree conditions [3]. Applying a tableau to a state is a backend kernel and does not live here.

Every rational coefficient is written as the exact rational of its reference and evaluated once, at import, so each constant is a single correctly-rounded division and no division happens per step.

References#

[1]

J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae”, J. Comput. Appl. Math. 6, 19 (1980). doi:10.1016/0771-050X(80)90013-3

[2]

Ch. Tsitouras, “Runge-Kutta pairs of order 5(4) satisfying only the first column simplifying assumption”, Comput. Math. Appl. 62, 770 (2011). doi:10.1016/j.camwa.2011.06.002

[3]

J. C. Butcher, Numerical Methods for Ordinary Differential Equations, 3rd ed., Wiley (2016), chapter 3.

Module Attributes

`MAX_CHECKED_ORDER`	Highest order for which `check_order_conditions()` knows the rooted trees.
`DOPRI5`	Dormand-Prince RK5(4)7M, the default scheme.
`DOPRI5_DENSE`	DOPRI5's dense-output weights `d`, Hairer, Noersett and Wanner's (II.6.12).
`TSIT5`	Tsitouras RK5(4), the alternative scheme.

Functions

`check_order_conditions`(tableau[, tolerance])	Measure the order the tableau's coefficients attain.
`continuous_extension`(tableau)	Return a tableau's dense-output weights, or `None` if it has none here.
`extension_coefficients`(tableau, start, ...)	Contract an accepted DOPRI5 step into its continuous extension's coefficients.
`extension_value`(coefficients, theta)	Evaluate a continuous extension at a fraction of its step.

Classes

`OrderConditionReport`(tolerance, ...)	The order a tableau's coefficients attain, measured tree by tree.
`Tableau`(name, nodes, matrix, ...)	An explicit Runge-Kutta tableau with an embedded companion.

solprop.stepper.tableau.DOPRI5 = Tableau(name='DOPRI5', nodes=array([0. , 0.2 , 0.3 , 0.8 , 0.88888889, 1. , 1. ]), matrix=array([[ 0. , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0.2 , 0. , 0. , 0. , 0. , 0. , 0. ], [ 0.075 , 0.225 , 0. , 0. , 0. , 0. , 0. ], [ 0.97777778, -3.73333333, 3.55555556, 0. , 0. , 0. , 0. ], [ 2.95259869, -11.59579332, 9.82289285, -0.29080933, 0. , 0. , 0. ], [ 2.84627525, -10.75757576, 8.90642272, 0.27840909, -0.2735313 , 0. , 0. ], [ 0.09114583, 0. , 0.4492363 , 0.65104167, -0.32237618, 0.13095238, 0. ]]), solution_weights=array([ 0.09114583, 0. , 0.4492363 , 0.65104167, -0.32237618, 0.13095238, 0. ]), error_weights=array([ 0.00123264, 0. , -0.00425277, 0.03697917, -0.0508638 , 0.04190476, -0.025 ]), order=5, embedded_order=4, first_same_as_last=True)#: Dormand-Prince RK5(4)7M, the default scheme.

solprop.stepper.tableau.DOPRI5_DENSE = array([-1.12701757, 0. , 2.67542448, -5.68552696, 3.52193237, -1.76728126, 2.38246893])#

DOPRI5’s dense-output weights d, Hairer, Noersett and Wanner’s (II.6.12).

With r1 = y0, r2 = y1 - y0, r3 = h k1 - r2, r4 = r2 - h k7 - r3 and r5 = h sum_i d_i k_i, the fourth-order continuous extension at theta in [0, 1] is r1 + theta (r2 + (1 - theta) (r3 + theta (r4 + (1 - theta) r5))) (continuous_extension(); SEP 14 phase 3).

solprop.stepper.tableau.MAX_CHECKED_ORDER = 5#: Highest order for which check_order_conditions() knows the rooted trees.

class solprop.stepper.tableau.OrderConditionReport(tolerance, solution_residuals, embedded_residuals, solution_order, embedded_order, max_residual, matches_declared_orders)[source]#

Bases: object

The order a tableau’s coefficients attain, measured tree by tree.

Attributes:

tolerancefloat: Residual below which a condition counts as satisfied, dimensionless.
solution_residualscollections.abc.Mapping: Absolute residual |sum_i b_i Phi_i(t) - 1/gamma(t)| of every rooted tree through order MAX_CHECKED_ORDER, keyed by tree label, for the propagated weights. Dimensionless.
embedded_residualscollections.abc.Mapping: The same for the embedded companion b_hat.
solution_orderint: Highest order, at most MAX_CHECKED_ORDER, whose conditions the propagated weights all satisfy within tolerance.
embedded_orderint: The same for the embedded companion.
max_residualfloat: Largest residual among the conditions the tableau declares it satisfies (up to Tableau.order for the propagated weights and Tableau.embedded_order for the companion). Independent of tolerance, so it is a measurement rather than a restatement of the threshold.
matches_declared_ordersbool: Whether the attained orders equal the declared ones.

embedded_order: int#

embedded_residuals: Mapping[str, float]#

matches_declared_orders: bool#

max_residual: float#

solution_order: int#

solution_residuals: Mapping[str, float]#

tolerance: float#

solprop.stepper.tableau.TSIT5 = Tableau(name='Tsit5', nodes=array([0. , 0.161 , 0.327 , 0.9 , 0.98002554, 1. , 1. ]), matrix=array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 1.61000000e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [-8.48065549e-03, 3.35480655e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 2.89715306e+00, -6.35944849e+00, 4.36229543e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 5.32586483e+00, -1.17488836e+01, 7.49553934e+00, -9.24950664e-02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [ 5.86145544e+00, -1.29209693e+01, 8.15936790e+00, -7.15849733e-02, -2.82690504e-02, 0.00000000e+00, 0.00000000e+00], [ 9.64607668e-02, 1.00000000e-02, 4.79889650e-01, 1.37900857e+00, -3.29006952e+00, 2.32471052e+00, 0.00000000e+00]]), solution_weights=array([ 0.09646077, 0.01 , 0.47988965, 1.37900857, -3.29006952, 2.32471052, 0. ]), error_weights=array([-0.00178001, -0.00081643, 0.00788088, -0.14471101, 0.58235717, -0.45808211, 0.01515152]), order=5, embedded_order=4, first_same_as_last=True)#: Tsitouras RK5(4), the alternative scheme.

class solprop.stepper.tableau.Tableau(name, nodes, matrix, solution_weights, error_weights, order, embedded_order, first_same_as_last)[source]#

Bases: object

An explicit Runge-Kutta tableau with an embedded companion.

The stage derivatives are k_i = f(z + c_i h, y + h sum_j a_ij k_j); the step is y_new = y + h sum_i b_i k_i and the local error estimate is e = h sum_i (b_i - b_hat_i) k_i.

Equality is by identity (eq=False): the fields are arrays, for which == is elementwise and would make a generated __eq__ raise.

Parameters:

namestr: Name of the scheme, used in error messages.
nodesnumpy.ndarray: Stage abscissae c of shape (n_stages,), as fractions of the step; c[0] is 0 for an explicit method.
matrixnumpy.ndarray: Butcher matrix A of shape (n_stages, n_stages), strictly lower triangular and dimensionless.
solution_weightsnumpy.ndarray: Weights b of shape (n_stages,) of the propagated solution, which is the one of order order. Dimensionless, summing to 1.
error_weightsnumpy.ndarray: Weights b - b_hat of shape (n_stages,) of the error estimate, dimensionless, summing to 0. The difference is stored rather than b_hat because it is the only combination ever used, and forming it once avoids a cancellation between two nearly equal vectors at every step.
orderint: Order of the propagated solution.
embedded_orderint: Order of the embedded companion b_hat, which sets the controller exponent.
first_same_as_lastbool: Whether the last stage derivative equals the first derivative of the next step (FSAL), so it can be reused instead of recomputed.

Raises:

StepperError: If the shapes disagree, an entry is not finite, A is not strictly lower triangular, a row of A does not sum to its node, b does not sum to 1, b - b_hat does not sum to 0 or is zero throughout, the orders are not 0 < embedded < order, or the FSAL flag is set but the last row of A is not b.

Notes

Two transcription errors are easy to make: swapping b and b_hat, which propagates the lower-order solution, and storing an error vector in the slot labelled as the higher-order weights. The two sum rules checked here catch both: an error vector does not sum to 1 and a weight vector does not sum to 0. See docs/dev/stepper.md section 10.

Examples

>>> from solprop.stepper import DOPRI5
>>> DOPRI5.n_stages, DOPRI5.order, DOPRI5.embedded_order
(7, 5, 4)
>>> DOPRI5.first_same_as_last
True
>>> bool(abs(DOPRI5.embedded_weights.sum() - 1.0) < 1e-15)
True

embedded_order: int#

property embedded_weights: Float[ndarray, 'n_stages'][source]#

Return the weights of the embedded companion.

Returns:

numpy.ndarray: Weights b_hat = b - (b - b_hat) of shape (n_stages,), dimensionless, summing to 1.

error_weights: Float[ndarray, 'n_error']#

first_same_as_last: bool#

matrix: Float[ndarray, 'n_rows n_cols']#

property n_stages: int[source]#

Return the number of stages, counting the reused FSAL stage.

Returns:

int: Number of stages.

name: str#

nodes: Float[ndarray, 'n_stages']#

order: int#

solution_weights: Float[ndarray, 'n_solution']#

solprop.stepper.tableau.check_order_conditions(tableau, tolerance=1e-15)[source]#

Measure the order the tableau’s coefficients attain.

A tableau of order p satisfies sum_i b_i Phi_i(t) = 1/gamma(t) for every rooted tree t of order up to p. There are 17 such trees through order 5. The check is applied to the propagated weights and to the embedded companion, so the labelling of the pair is verified rather than assumed, which catches a swapped or mislabelled pair.

This function reports; it does not raise. A caller that wants a hard failure asserts on matches_declared_orders.

Parameters:

tableauTableau: The tableau to measure.
tolerancefloat, optional: Residual below which a condition counts as satisfied, dimensionless. The default 1e-15 is above the worst residual of either shipped tableau (8.0e-16, in Tsit5’s order-5 conditions, whose coefficients are published decimals rather than exact rationals) and eleven orders below the residual of an order that is genuinely not attained (8.7e-4 for both embedded companions at order 5).

Returns:

OrderConditionReport: The residuals, the attained orders and whether they match the declared ones.

Examples

>>> from solprop.stepper import DOPRI5, check_order_conditions
>>> report = check_order_conditions(DOPRI5)
>>> report.solution_order, report.embedded_order
(5, 4)
>>> report.matches_declared_orders
True
>>> bool(report.max_residual < 1e-15)
True

solprop.stepper.tableau.continuous_extension(tableau)[source]#

Return a tableau’s dense-output weights, or None if it has none here.

Parameters:

tableauTableau: The scheme.

Returns:

numpy.ndarray or None: DOPRI5_DENSE for DOPRI5; None for any other tableau.

Examples

>>> from solprop.stepper.tableau import DOPRI5, TSIT5, continuous_extension
>>> continuous_extension(DOPRI5).shape, continuous_extension(TSIT5)
((7,), None)

solprop.stepper.tableau.extension_coefficients(tableau, start, stages, dz)[source]#

Contract an accepted DOPRI5 step into its continuous extension’s coefficients.

Parameters:

tableauTableau: DOPRI5, the scheme with a continuous extension here.
startnumpy.ndarray: The state at the step’s start, y0.
stagessequence of numpy.ndarray: The step’s stage derivatives k_i, one per stage.
dzfloat: The step h.

Returns:

list of numpy.ndarray: [y0, y1 - y0, h k1 - (y1 - y0), (y1 - y0) - h k7 - (h k1 - (y1 - y0)), h sum_i d_i k_i], with y1 = y0 + h sum_i b_i k_i.

Raises:

ValueError: If the tableau has no continuous extension here.

Examples

A constant derivative extends linearly: half a step of y' = 2 from 1 is 2.

>>> import numpy as np
>>> from solprop.stepper.tableau import (
...     DOPRI5, extension_coefficients, extension_value)
>>> stages = [np.array([2.0])] * 7
>>> coefficients = extension_coefficients(DOPRI5, np.array([1.0]), stages, 1.0)
>>> extension_value(coefficients, 0.5).round(12).tolist()
[2.0]

solprop.stepper.tableau.extension_value(coefficients, theta)[source]#

Evaluate a continuous extension at a fraction of its step.

Parameters:

coefficientssequence of numpy.ndarray: From extension_coefficients().
thetafloat: The fraction of the step, in [0, 1].

Returns:

numpy.ndarray: r1 + theta (r2 + (1 - theta) (r3 + theta (r4 + (1 - theta) r5))).