solprop.stepper.multistep#

Low-call Adams integrators: the method specification and its variable-step algebra.

SEP 13 (docs/dev/seps/archive/sep-0013-expensive-response-steppers.md) proposes interaction-picture Adams methods for propagation whose nonlinear evaluation dominates the cost. This module holds the backend-agnostic half: the frozen Multistep specification a run carries in place of a Runge-Kutta tableau, and the scalar algebra every backend shares – the integrated Lagrange weights of a variable-step Adams formula on the actual history positions, the leading defect constants and the Milne factor of a predictor-corrector pair (SEP 13 section 4).

Positions enter in normalised form, theta = (z - z_n) / h with z_n the anchor and h the step being taken, so a history point behind the anchor has a negative node and the predicted endpoint the node 1. The weights are exact for polynomials of degree below the node count, on any distinct nodes; constant-step Adams-Bashforth and Adams-Moulton coefficients are the special case of equally spaced nodes [1]. Variable coefficients recomputed on the actual history, rather than fixed ones, follow the variable-coefficient Adams methods of [2]; no solver code is taken from it.

References#

[1]

E. Hairer, S. P. Noersett and G. Wanner, Solving Ordinary Differential Equations I, 2nd ed., Springer (1993), sections III.1 and III.5.

[2]

P. N. Brown, G. D. Byrne and A. C. Hindmarsh, “VODE: A Variable-Coefficient ODE Solver”, SIAM J. Sci. Stat. Comput. 10, 1038 (1989). doi:10.1137/0910062

Module Attributes

`MAX_ADAMS_ORDER`	Highest admitted order; SEP 13 section 4.4 bounds the initial study at six.
`PEC_REAL_INTERVAL`	Negative-real-axis stability interval of PEC `q/(q+1)` at constant step, by `q`.
`Interpolant`	How a comparison method interpolates save planes when they are interpolated.
`Formulation`	How an Adams method integrates its source history (SEP 16 section 4.2).
`Corrector`	The execution contract of a step (SEP 13 section 4.2).

Functions

`defect_integral`(nodes[, upper])	Return the integral from 0 to upper of the node polynomial.
`integrated_weights`(nodes[, upper, check, ...])	Return the weights of the integral from 0 to upper of an interpolant.
`lagrange_coefficients`(nodes, *[, check, helpers])	Return the monomial coefficients of the Lagrange basis on distinct nodes.
`milne_factor`(predictor_nodes, corrector_nodes)	Return the factor turning a predictor-corrector difference into the defect.

Classes

`ComparisonMultistep`([order, variable_order, ...])	An Adams variant SEP 13 compared and did not select, kept for benchmarks.
`Multistep`([order, variable_order, ...])	A low-call Adams method, in place of a Runge-Kutta tableau.

class solprop.stepper.multistep.ComparisonMultistep(order=5, variable_order=False, max_calls=None, interval_cap=False, equal_approach=True, max_growth=2.0, min_factor=0.2, safety=0.8, formulation='interaction', corrector='pece', corrector_order=None, interpolant='adams')[source]#

Bases: Multistep

An Adams variant SEP 13 compared and did not select, kept for benchmarks.

PECE, explicit Adams-Bashforth and the matched predictor-corrector pairs controlled by Milne’s estimate each lost to Multistep on every case measured (SEP 13 Part II). They stay executable so that the comparison can be rerun, and are not part of the documented interface.

Attributes:

corrector{‘none’, ‘pec’, ‘pece’}: The execution contract (Corrector).
corrector_orderint or None: Derivatives of the corrector, the predicted one included: order for a matched pair or order + 1, the selected family. None is order. Ignored for corrector='none'.
interpolant{‘adams’, ‘hermite’}: How the planes are interpolated when the run interpolates them (interpolate_planes, ADR-036, ADR-039): the accepted step’s own derivative polynomial, at the step’s order and no call, or the generic cubic Hermite interpolant of the step’s end states and derivatives, a reference whose error is one order below what the step controls (Interpolant).
formulation{‘interaction’}: The variants are stepped in the interaction picture only.

Raises:

StepperError: If a field is outside its range, variable order is asked of a matched pair, or the exponential formulation is asked of a comparison variant (SEP 16 delivers it for the selected family only).

Examples

>>> from solprop.stepper.multistep import ComparisonMultistep
>>> ComparisonMultistep(corrector="pece", order=4).name
'adams-pece4'

corrector: Literal['none', 'pec', 'pece'] = 'pece'#

corrector_nodes(order)[source]#

Return the corrector’s derivative count at a current predictor order.

Parameters:

orderint: The predictor order in force.

Returns:

int: order for a matched pair, order + 1 for the higher-order corrector.

corrector_order: int | None = None#

formulation: Literal['interaction', 'exponential'] = 'interaction'#

interpolant: Literal['adams', 'hermite'] = 'adams'#

property name: str[source]#

Return the variant’s record name.

Returns:

str: adams-<corrector><order>, with /<corrector order> for the higher-order corrector, ab standing for none.

variable_order: bool = False#

solprop.stepper.multistep.Corrector#

The execution contract of a step (SEP 13 section 4.2).

"pec" corrects once and stores the derivative of the predicted state, the selected method’s contract. "pece" pays a second call for the derivative of the corrected, accepted state, and "none" is explicit Adams-Bashforth; both remain only as comparators (ComparisonMultistep).

alias of Literal[‘none’, ‘pec’, ‘pece’]

solprop.stepper.multistep.Formulation#

How an Adams method integrates its source history (SEP 16 section 4.2).

"interaction" interpolates each source sample carried into the step’s anchor frame by the exact linear flow (Lawson’s interaction picture, SEP 13). "exponential", Multistep’s default, interpolates the untransported samples and integrates the polynomial exactly against the linear flow (exponential Adams).

alias of Literal[‘interaction’, ‘exponential’]

solprop.stepper.multistep.Interpolant#

How a comparison method interpolates save planes when they are interpolated.

"adams" integrates the accepted step’s own derivative polynomial, the selected method’s interpolant; "hermite" is the generic cubic Hermite reference (SEP 13 section 6.4). Whether planes are interpolated at all is the run’s interpolate_planes.

alias of Literal[‘adams’, ‘hermite’]

solprop.stepper.multistep.MAX_ADAMS_ORDER = 6#: Highest admitted order; SEP 13 section 4.4 bounds the initial study at six.

class solprop.stepper.multistep.Multistep(order=5, variable_order=True, max_calls=None, interval_cap=False, equal_approach=True, max_growth=2.0, min_factor=0.2, safety=0.8, formulation='exponential')[source]#

Bases: object

A low-call Adams method, in place of a Runge-Kutta tableau.

SEP 13’s selected method: one source call per accepted step, an Adams-Bashforth predictor of order q and an Adams-Moulton corrector one order higher through the new source (PEC), controlled on the predictor-corrector difference. Its predicted states advance by Adams-Bashforth of order q + 1, so it has that formula’s stability, and the accepted value is the corrector’s (docs/guide/numerics/adams.md). Two formulations share that recurrence, its history, controller and lifecycle (SEP 16): exponential Adams, the default, and the interaction picture (Formulation).

Attributes:

orderint: Highest predictor order q in [1, 6]; the corrector uses q + 1 derivatives. The method starts at order one and raises it as history accrues.
variable_orderbool: Choose the order between one and order step by step from the neighbouring orders’ estimates, at no extra call (SEP 13 section 4.4). On by default: it settled at the best fixed order on every case SEP 13 measured.
max_callsint or None: Budget of charged source calls for the run: startup, accepted and rejected trials and restarts alike. A trial the remaining budget cannot pay for is not started; the run stops with solprop.exceptions.ResponseBudgetError at its last accepted step. None for no budget.
interval_capbool: Cap each step at the width of the Chebyshev linear-table interval containing its start, as the Runge-Kutta march does. Off by default: the cumulative operator is exact at any position, and on SEP 13’s canonical capillary the cap – dozens of micrometre intervals at a choked vacuum end – cost 20-40 % of the calls for no accuracy.
equal_approachbool: Divide the distance to a forced landing into two equal steps when it lies between one and two proposals, instead of a full step and a sliver.
max_growthfloat: Largest ratio of one step to the step before, above one. Variable-step Adams stability needs it bounded (SEP 13 section 4.4).
min_factorfloat: Smallest step-size factor after a rejection, in (0, 1).
safetyfloat: Safety factor on the step proposed from an error estimate, in (0, 1].
formulation{‘interaction’, ‘exponential’}: How the source history is integrated (Formulation); exponential by default. A run’s record hashes the method’s representation, which names it (output schema 15); a record made before then that names none was stepped in the interaction picture.

Raises:

StepperError: If a field is outside its range.

Examples

>>> from solprop.stepper.multistep import Multistep
>>> Multistep().name
'exponential-adams-pec5/6-variable'
>>> Multistep(order=3, variable_order=False, formulation="interaction").name
'adams-pec3/4'

property corrector: Literal['none', 'pec', 'pece'][source]#

Return the execution contract: one call per step, the predicted source kept.

Returns:

str: "pec".

corrector_nodes(order)[source]#

Return the corrector’s derivative count at a current predictor order.

Parameters:

orderint: The predictor order in force, which is below order during startup.

Returns:

int: order + 1: the corrector is one order higher than the predictor.

equal_approach: bool = True#

formulation: Literal['interaction', 'exponential'] = 'exponential'#

interval_cap: bool = False#

max_calls: int | None = None#

max_growth: float = 2.0#

min_factor: float = 0.2#

property name: str[source]#

Return the method’s record name.

Returns:

str: adams-pec<q>/<q + 1>, with -variable when the order is chosen and an exponential- prefix for the exponential formulation.

order: int = 5#

property real_stability_interval: float[source]#

Return the negative-real stability interval of the highest order run.

Returns:

float: PEC_REAL_INTERVAL at order: the largest h |lambda| for a damping rate lambda integrated explicitly.

safety: float = 0.8#

variable_order: bool = True#

solprop.stepper.multistep.PEC_REAL_INTERVAL = (1.0, 0.545, 0.3, 0.163, 0.0877, 0.0465)#

Negative-real-axis stability interval of PEC q/(q+1) at constant step, by q.

The predicted states of the one-order-higher PEC family advance by Adams-Bashforth of order q + 1 (SEP 13), so the interval is that formula’s: h |lambda| up to 1, 6/11, 0.3, 0.1633, 0.0877 and 0.0465 for q = 1..6 (rounded down), against 3.3 for DOPRI5. Computed from the spectral radius of the PEC recurrence for u' = lambda u in tests/stepper/test_multistep.py. It binds a damping term integrated explicitly – the in-stage collar – in either formulation: the exponential one integrates the diagonal operator exactly, but the collar is in its source.

solprop.stepper.multistep.defect_integral(nodes, upper=1.0)[source]#

Return the integral from 0 to upper of the node polynomial.

integral_0^upper prod_j (theta - theta_j) dtheta. Divided by q! for q nodes it is the leading defect constant C of the formula interpolating on those nodes, exact - formula = C h^{q+1} g^{(q)} (SEP 13 section 4.3): the interpolation remainder, integrated.

Parameters:

nodessequence of float: Normalised positions, as for integrated_weights().
upperfloat, optional: Upper limit of the integral.

Returns:

float: The integral, dimensionless.

Examples

The constant-step AB4 and AM4 constants, 251/720 and -19/720, times 4!:

>>> from solprop.stepper.multistep import defect_integral
>>> round(defect_integral([0, -1, -2, -3]) * 30, 9)
251.0
>>> round(defect_integral([1, 0, -1, -2]) * 30, 9)
-19.0

solprop.stepper.multistep.integrated_weights(nodes, upper=1.0, *, check=True, helpers=Helpers(name='numba'))[source]#

Return the weights of the integral from 0 to upper of an interpolant.

sum_j w_j f(theta_j) equals integral_0^upper p(theta) dtheta for the polynomial p of degree len(nodes) - 1 interpolating f at the nodes: each weight is the integral of the j-th Lagrange basis polynomial. With the nodes in units of the step, h * sum_j w_j g_j is the Adams increment.

Parameters:

nodessequence of float: Distinct normalised positions (z_j - z_n) / h.
upperfloat, optional: Upper limit of the integral in the same units; 1 for a whole step, a fraction of it for dense output.
checkbool, optional: Whether to refuse empty, non-finite or repeated nodes. The Adams march, whose nodes are accepted positions and so distinct, skips it on every step.
helperssolprop._helpers.Helpers, optional: The evaluating backend’s implementation of the basis (SEP 29 section 6.6); the shared Numba one by default, for callers outside a backend.

Returns:

numpy.ndarray: One weight per node.

Raises:

StepperError: If the nodes are empty or not distinct.

Examples

Constant-step Adams-Bashforth 4 and Adams-Moulton 4, times 24:

>>> from solprop.stepper.multistep import integrated_weights
>>> (24 * integrated_weights([0, -1, -2, -3])).round(12).tolist()
[55.0, -59.0, 37.0, -9.0]
>>> (24 * integrated_weights([1, 0, -1, -2])).round(12).tolist()
[9.0, 19.0, -5.0, 1.0]

solprop.stepper.multistep.lagrange_coefficients(nodes, *, check=True, helpers=Helpers(name='numba'))[source]#

Return the monomial coefficients of the Lagrange basis on distinct nodes.

Row j holds c_jk with ell_j(theta) = sum_k c_jk theta**k, ascending in k, so that ell_j(theta_i) is one at i = j and zero elsewhere. The exponential Adams weights integrate each power against the linear flow (w_j(x) = sum_k c_jk k! phi_{k+1}(x), SEP 16 section 4.2); at x = 0 they reduce to integrated_weights().

Parameters:

nodessequence of float: Distinct normalised positions (z_j - z_n) / h.
checkbool, optional: Whether to refuse empty, non-finite or repeated nodes.
helperssolprop._helpers.Helpers, optional: The evaluating backend’s implementation of the basis; the shared Numba one by default.

Returns:

numpy.ndarray: The coefficients, shape (n_nodes, n_nodes).

Raises:

StepperError: If the nodes are empty or not distinct.

Examples

The linear basis on the nodes 1 and 0: theta and 1 - theta.

>>> from solprop.stepper.multistep import lagrange_coefficients
>>> lagrange_coefficients([1.0, 0.0]).tolist()
[[0.0, 1.0], [1.0, -1.0]]

solprop.stepper.multistep.milne_factor(predictor_nodes, corrector_nodes)[source]#

Return the factor turning a predictor-corrector difference into the defect.

With the predictor’s defect C_B h^{q+1} g^{(q)} and the corrector’s C_M h^{q+1} g^{(q)} on the same number of nodes, v_C - v_P is (C_B - C_M) h^{q+1} g^{(q)}, so the corrector’s local defect is C_M / (C_B - C_M) times the difference (Milne’s device, SEP 13 section 4.3). The common 1/q! cancels.

Parameters:

predictor_nodes, corrector_nodessequence of float: The two formulas’ normalised nodes, equal in number; the corrector’s include 1.

Returns:

float: C_M / (C_B - C_M).

Raises:

StepperError: If the node counts differ, or the two constants coincide so the difference carries no information.

Examples

>>> from solprop.stepper.multistep import milne_factor
>>> round(milne_factor([0, -1, -2, -3], [1, 0, -1, -2]) * 270, 9)
-19.0