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`.
`RESTART_AFTER`	Consecutive rejections after which the history is restarted at order one.
`ORDER_HYSTERESIS`	Factor by which a neighbouring order's predicted step must win to be adopted.
`COLLAR_STABILITY`	Share of an order's negative-real stability interval the order choice allows.
`REJECTION_GROWTH`	Largest factor a rejected step's own defect may propose for its retry.
`CARRIER_PHASE_CAP`	Largest carrier phase `h \|kappa\|` a step may cover under a carrier frame, rad.
`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).
`CarrierFrame`	Which carrier phase the exponential formulation removes from its history (SEP 35).
`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.
`march_policy`(method, *, dz_min, dz_max, ...)	Return the policy of one run of method (`MarchPolicy`).
`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.
`MarchPolicy`(order, variable_order, safety, ...)	Every number the Adams march reads and never changes during a run (SEP 39).
`Multistep`([order, variable_order, ...])	A low-call Adams method, in place of a Runge-Kutta tableau.

solprop.stepper.multistep.CARRIER_PHASE_CAP = 1.0#

Largest carrier phase h |kappa| a step may cover under a carrier frame, rad.

Part of the CarrierFrame policy, not a setting (SEP 35 section 4.3): the residual rotation per step of the source components the shift does not remove – stationary ones at kappa, the third harmonic’s at 2 kappa – then stays at or below 1 and 2 rad, short of the pi at which the samples alias and the defect stops seeing the error. The proposal is shortened before a trial; a shortened step is not a landing.

solprop.stepper.multistep.COLLAR_STABILITY = 1.0#: Share of an order’s negative-real stability interval the order choice allows.

solprop.stepper.multistep.CarrierFrame#

Which carrier phase the exponential formulation removes from its history (SEP 35).

On a field grid, and on an envelope grid with the third harmonic, the operator keeps beta_0 and the source that follows the carrier rotates along z at kappa = beta(omega_0) - beta_1 omega_0. "scalar" interpolates, inside each step, G_j = exp(sigma (1 - theta_j)) F_j with sigma = i Im x at the reference bin, and integrates it against the moments at x - sigma: exact for any sigma, with the state, the stored history and the linear flow untouched, and each step capped at CARRIER_PHASE_CAP. "off", the default, is SEP 16’s arithmetic. On an envelope grid without the third harmonic the shift is the present arithmetic exactly (sigma = 0). The shift saves source calls on dense gases at matched total-field error, but the weak bands it leaves less accurate (SEP 35 phase 4, docs/guide/numerics/adams.md), which is why it is not the default. A per-band shift of the third harmonic was developed and removed by SEP 35’s gate.

alias of Literal[‘off’, ‘scalar’]

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', carrier_frame='off', 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'#

property history_key: str[source]#

Return what defines the history a march of this variant keeps.

Every field of a comparison variant may shape its history, so the whole representation: any change restarts it.

Returns:

str: The variant’s representation.

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.MarchPolicy(order, variable_order, safety, max_growth, min_factor, rejection_growth, hysteresis, restart_after, collar_stability, real_intervals, collar_rate, carrier, carrier_cap, max_calls, equal_approach, interval_cap, dz_min, dz_max, stretch, z_end, save_planes, hard_kinks, breakpoints, landing_ulps, interpolate, quiet_interval)[source]#

Bases: object

Every number the Adams march reads and never changes during a run (SEP 39).

solprop forms one per run, from the method, the run’s stepper settings, the operator’s breakpoints and kinks, the in-stage collar and this module’s constants (march_policy()). The Python march reads its controller and landing numbers from it. The native march receives it whole and holds no default and no constant of its own, so every number either execution uses has this one source (SEP 39 section 6.2).

Attributes:

orderint: Highest predictor order, q_max.
variable_orderbool: Whether the order is chosen step by step.
safety, max_growth, min_factorfloat: The method’s safety factor, largest growth and smallest factor.
rejection_growthfloat: Largest factor a rejection’s own defect may propose (REJECTION_GROWTH).
hysteresisfloat: ORDER_HYSTERESIS.
restart_afterint: RESTART_AFTER.
collar_stabilityfloat: COLLAR_STABILITY.
real_intervalstuple of float: PEC_REAL_INTERVAL, by order.
collar_ratefloat: The in-stage collar’s largest rate over two, 1/m; zero without one, outside the exponential formulation, or under the split scheme.
carrierbool: Whether a carrier frame caps each proposal (CarrierFrame).
carrier_capfloat: CARRIER_PHASE_CAP, rad.
max_callsint or None: The call budget.
equal_approach, interval_capbool: The method’s landing options.
dz_min, dz_max, stretch, z_endfloat: The stepper settings: the floor, the cap, the landing stretch and the run’s end, m (the stretch is a ratio).
save_planestuple of float: The planes to land on, m; empty when the planes are interpolated.
hard_kinkstuple of float: The operator’s hard kinks, m.
breakpointstuple of float: The edges whose intervals cap a step, m: the operator’s breakpoints with the interval cap, otherwise only its first and last.
landing_ulpsint: LANDING_ULPS.
interpolatebool: Whether the save planes are observed densely.
quiet_intervalfloat or None: The boundary monitor’s interval, m, which caps a step while the collar rests; None without a monitor. The Python march reads the same number from the monitor through the driver (StepDriver._quiet_capped).

Raises:

StepperError: If an order, factor, interval, ulp count or position set is out of range.

breakpoints: tuple[float, ...]#

carrier: bool#

carrier_cap: float#

collar_rate: float#

collar_stability: float#

dz_max: float#

dz_min: float#

equal_approach: bool#

hard_kinks: tuple[float, ...]#

hysteresis: float#

interpolate: bool#

interval_cap: bool#

landing_ulps: int#

max_calls: int | None#

max_growth: float#

min_factor: float#

order: int#

quiet_interval: float | None#

real_intervals: tuple[float, ...]#

rejection_growth: float#

restart_after: int#

safety: float#

save_planes: tuple[float, ...]#

stretch: float#

variable_order: bool#

z_end: float#

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', carrier_frame='off')[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.
carrier_frame{‘off’, ‘scalar’}: The carrier phase the exponential formulation removes from its history inside each step (CarrierFrame, docs/guide/numerics/adams.md); "off" by default. The shift needs the exponential formulation and is refused by name on a backend or route that does not deliver it. The representation names it, so a march snapshot taken under another policy restarts its history.

Raises:

StepperError: If a field is outside its range, or a carrier frame is asked of the interaction formulation.

Examples

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

carrier_frame: Literal['off', 'scalar'] = 'off'#

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'#

property history_key: str[source]#

Return what defines the history a march of this method keeps.

A resume continues an Adams history only when the snapshot’s key equals the resuming method’s: the formulation, the carrier frame and the highest order define the samples a history holds and how they combine. The call budget, the safety factor and the other step-size settings choose the next step, not what the history means, so a change of them keeps it.

Returns:

str: name without the variable-order marker.

Examples

>>> from solprop.stepper.multistep import Multistep
>>> Multistep(max_calls=10).history_key == Multistep(safety=0.5).history_key
True
>>> Multistep().history_key == Multistep(order=4).history_key
False

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, an exponential- prefix for the exponential formulation and a -<policy>-frame suffix under a carrier frame.

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.ORDER_HYSTERESIS = 1.2#: Factor by which a neighbouring order’s predicted step must win to be adopted.

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.

The interval is the one at a residual phase Im x' = 0, x' the argument of the step’s moments: h D in the exponential formulation, h D - sigma under a carrier frame (CarrierFrame). In a bin whose residual phase per step is nonzero the recurrence is the same with x' in place of x; every h |lambda| up to the interval stays stable for |Im x'| up to about 5.9, 5.1, 3.5, 3.2 and 3.0 at q = 1..5 without loss, and further with it (SEP 35 section 4.4, reproduced in tests/backends/test_carrier_frame.py). Neither formulation enforces that condition; the temporal collar couples bins, so the bound is a heuristic in both.

solprop.stepper.multistep.REJECTION_GROWTH = 0.9#

Largest factor a rejected step’s own defect may propose for its retry.

Below one, so that a retry is always shorter than the step it replaces, whatever the defect’s exponent makes of a norm just above one.

solprop.stepper.multistep.RESTART_AFTER = 3#: Consecutive rejections after which the history is restarted at order one.

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='plain'))[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 plain Python forms 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='plain'))[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 plain Python forms 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.march_policy(method, *, dz_min, dz_max, stretch, z_end, save_planes, hard_kinks, breakpoints, interpolate, collar_rate=0.0, carrier=False, quiet_interval=None)[source]#

Return the policy of one run of method (MarchPolicy).

Parameters:

methodMultistep: The method.
dz_min, dz_max, stretch, z_endfloat: The run’s stepper settings.
save_planes, hard_kinks, breakpointssequence of float: The run’s save planes and the operator’s kinks and breakpoints, m, in order.
interpolatebool: Whether the planes are interpolated rather than landed on.
collar_ratefloat, optional: The in-stage collar’s largest rate over two, 1/m.
carrierbool, optional: Whether a carrier frame is active.
quiet_intervalfloat or None, optional: The boundary monitor’s interval, m.

Returns:

MarchPolicy: The record, with this module’s constants filled in.

Examples

>>> from solprop.stepper.multistep import Multistep, march_policy
>>> policy = march_policy(
...     Multistep(), dz_min=1e-9, dz_max=1e-2, stretch=1.1, z_end=1.0,
...     save_planes=[0.5, 1.0], hard_kinks=[], breakpoints=[0.0, 0.3, 1.0],
...     interpolate=False,
... )
>>> policy.order, policy.hysteresis, policy.breakpoints
(5, 1.2, (0.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