Solvers

JNDF, JBDF, and JIDA are diff-soup's variable-step, variable-order implicit engines, and DiffraxStepper wraps an explicit Diffrax solver for identity-mass problems. See ODE/DAE solvers for how to choose.

Recommended defaults

Construct engines through the stock preset, the profile of every bare constructor, as in ds.JNDF.stock(rtol=1e-6, atol=1e-6, linear="auto"). The flash preset replaces Jacobian reuse across steps with a fresh Jacobian each corrector iteration. Leave profile to the preset.

---

Implicit engines

diff_soup.JNDF

The numerical differentiation formula (NDF) solver, the method of MATLAB's ode15s.

It solves ODEs and mass-form DAEs M(u, t) du/dt = f, including a state-dependent mass.

Attributes:

Name	Type	Description
adjoint	Any	The default differentiation strategy, where None is forward mode and CheckpointAdjoint() supports jax.grad.

adjoint class-attribute instance-attribute

adjoint: Any = eqx.field(static=True, default=None)

stock classmethod

stock(**kwargs: Any) -> '_VSVOSolver'

Build the solver with the stock profile, the reference method's step control and Jacobian reuse.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Constructor settings, which override the profile.	{}

Returns:

Type	Description
'_VSVOSolver'	The solver.

flash classmethod

flash(**kwargs: Any) -> '_VSVOSolver'

Build the solver with the flash profile, proportional-integral step control and a fresh Jacobian.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Constructor settings, which override the profile.	{}

Returns:

Type	Description
'_VSVOSolver'	The solver.

solve

solve(problem: Any, t_span: Any, y0: Any, *, saveat: Any = None, linear: Any = None, hmin: float = 0.0, hmax: float = 0.0, dt0: float | Any = 0.0, max_order: int = 5, max_steps: int = 10000, max_events: int = 16, event: Any = None, control_state0: Any = None, theta: Any = None, consistent_ic: Any = None, lambda0: Any = None, adjoint: Any = None, tolerances: Any = None, trace: bool = False, history: Any = None) -> Solution

Integrate problem over t_span from y0.

A repeat solve skips compilation for new values of t_span, y0, theta, control_state0 and dt0. Keep the problem object stable across calls, since a new one compiles again.

Parameters:

Name	Type	Description	Default
problem	Any	The built problem.	required
t_span	Any	The interval (t0, t1).	required
y0	Any	The initial state.	required
saveat	Any	The output times, or None for the end of the interval.	None
linear	Any	A linear backend for this solve, overriding the solver's.	None
hmin	float	The smallest step size, 0 to leave it to the engine.	0.0
hmax	float	The largest step size, 0 to leave it to the engine.	0.0
dt0	float | Any	The first step size, 0 to estimate it.	0.0
max_order	int	The highest integration order, at most 5.	5
max_steps	int	The step budget, which also sizes the output buffers.	10000
max_events	int	The capacity of the event log, past which event_log.overflow is set.	16
event	Any	An event or event set to monitor.	None
control_state0	Any	The initial control values, {name: value} for a Recipe with bound controls.	None
theta	Any	The parameters, or None for the problem's defaults.	None
consistent_ic	Any	True or False to decide the consistent initialization of this solve, None for the solver's default.	None
lambda0	Any	The initial values of the global unknowns, when y0 holds only the field block.	None
adjoint	Any	The differentiation strategy, or None for the solver's.	None
tolerances	Any	A ToleranceSpec, overriding the solver's scalar tolerances.	None
trace	bool	True attaches a per-step diagnostic trace.	False
history	Any	True records Solution.history, and an earlier solve's Solution.history continues from its end.	None

Returns:

Type	Description
Solution	The Solution.

Raises:

Type	Description
ValueError	If the problem needs a capability the engine does not have.
ConfigurationError	If the engine cannot honour a setting, such as algebraic_error_control="exclude" off JIDA.

diff_soup.JBDF

The backward differentiation formula (BDF) solver, with a fixed-leading-coefficient corrector.

It solves ODEs and mass-form DAEs M(t) du/dt = f with a constant or time-dependent mass.

Attributes:

Name	Type	Description
adjoint	Any	The default differentiation strategy, where None is forward mode and CheckpointAdjoint() supports jax.grad.

adjoint class-attribute instance-attribute

adjoint: Any = eqx.field(static=True, default=None)

stock classmethod

stock(**kwargs: Any) -> '_VSVOSolver'

Build the solver with the stock profile, the reference method's step control and Jacobian reuse.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Constructor settings, which override the profile.	{}

Returns:

Type	Description
'_VSVOSolver'	The solver.

flash classmethod

flash(**kwargs: Any) -> '_VSVOSolver'

Build the solver with the flash profile, proportional-integral step control and a fresh Jacobian.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Constructor settings, which override the profile.	{}

Returns:

Type	Description
'_VSVOSolver'	The solver.

solve

solve(problem: Any, t_span: Any, y0: Any, *, saveat: Any = None, linear: Any = None, hmin: float = 0.0, hmax: float = 0.0, dt0: float | Any = 0.0, max_order: int = 5, max_steps: int = 10000, max_events: int = 16, event: Any = None, control_state0: Any = None, theta: Any = None, consistent_ic: Any = None, lambda0: Any = None, adjoint: Any = None, tolerances: Any = None, trace: bool = False, history: Any = None) -> Solution

Integrate problem over t_span from y0.

A repeat solve skips compilation for new values of t_span, y0, theta, control_state0 and dt0. Keep the problem object stable across calls, since a new one compiles again.

Parameters:

Name	Type	Description	Default
problem	Any	The built problem.	required
t_span	Any	The interval (t0, t1).	required
y0	Any	The initial state.	required
saveat	Any	The output times, or None for the end of the interval.	None
linear	Any	A linear backend for this solve, overriding the solver's.	None
hmin	float	The smallest step size, 0 to leave it to the engine.	0.0
hmax	float	The largest step size, 0 to leave it to the engine.	0.0
dt0	float | Any	The first step size, 0 to estimate it.	0.0
max_order	int	The highest integration order, at most 5.	5
max_steps	int	The step budget, which also sizes the output buffers.	10000
max_events	int	The capacity of the event log, past which event_log.overflow is set.	16
event	Any	An event or event set to monitor.	None
control_state0	Any	The initial control values, {name: value} for a Recipe with bound controls.	None
theta	Any	The parameters, or None for the problem's defaults.	None
consistent_ic	Any	True or False to decide the consistent initialization of this solve, None for the solver's default.	None
lambda0	Any	The initial values of the global unknowns, when y0 holds only the field block.	None
adjoint	Any	The differentiation strategy, or None for the solver's.	None
tolerances	Any	A ToleranceSpec, overriding the solver's scalar tolerances.	None
trace	bool	True attaches a per-step diagnostic trace.	False
history	Any	True records Solution.history, and an earlier solve's Solution.history continues from its end.	None

Returns:

Type	Description
Solution	The Solution.

Raises:

Type	Description
ValueError	If the problem needs a capability the engine does not have.
ConfigurationError	If the engine cannot honour a setting, such as algebraic_error_control="exclude" off JIDA.

diff_soup.JIDA

The implicit DAE solver of the SUNDIALS IDA family, on the residual form F(u, du/dt, t) = 0.

It solves ODEs and fully implicit index-1 DAEs, and it is the recommended solver.

Attributes:

Name	Type	Description
consistent_ic	bool	True (default) projects the initial state before the first step, a no-op for a pure ODE.
ic_mode	str	The initialization route, "ya_ydp" (default), "ya_ydp_euler", "y_fixed_ydot" or "assume_consistent".
suppress_alg	bool	True (default) keeps the algebraic rows out of the local error test.
id	Any	The differential-row mask of the error test, by default from the problem's mass structure.
adjoint	Any	The default differentiation strategy, where None is forward mode and CheckpointAdjoint() supports jax.grad.

consistent_ic class-attribute instance-attribute

consistent_ic: bool = eqx.field(static=True, default=True)

ic_mode class-attribute instance-attribute

ic_mode: str = eqx.field(static=True, default='ya_ydp')

suppress_alg class-attribute instance-attribute

suppress_alg: bool = eqx.field(static=True, default=True)

id class-attribute instance-attribute

id: Any = eqx.field(static=True, default=None, converter=_normalize_id_mask)

adjoint class-attribute instance-attribute

adjoint: Any = eqx.field(static=True, default=None)

stock classmethod

stock(**kwargs: Any) -> '_VSVOSolver'

Build the solver with the stock profile, the reference method's step control and Jacobian reuse.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Constructor settings, which override the profile.	{}

Returns:

Type	Description
'_VSVOSolver'	The solver.

flash classmethod

flash(**kwargs: Any) -> '_VSVOSolver'

Build the solver with the flash profile, proportional-integral step control and a fresh Jacobian.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Constructor settings, which override the profile.	{}

Returns:

Type	Description
'_VSVOSolver'	The solver.

solve

solve(problem: Any, t_span: Any, y0: Any, *, saveat: Any = None, linear: Any = None, hmin: float = 0.0, hmax: float = 0.0, dt0: float | Any = 0.0, max_order: int = 5, max_steps: int = 10000, max_events: int = 16, event: Any = None, control_state0: Any = None, theta: Any = None, consistent_ic: Any = None, lambda0: Any = None, adjoint: Any = None, tolerances: Any = None, trace: bool = False, history: Any = None) -> Solution

Integrate problem over t_span from y0.

A repeat solve skips compilation for new values of t_span, y0, theta, control_state0 and dt0. Keep the problem object stable across calls, since a new one compiles again.

Parameters:

Name	Type	Description	Default
problem	Any	The built problem.	required
t_span	Any	The interval (t0, t1).	required
y0	Any	The initial state.	required
saveat	Any	The output times, or None for the end of the interval.	None
linear	Any	A linear backend for this solve, overriding the solver's.	None
hmin	float	The smallest step size, 0 to leave it to the engine.	0.0
hmax	float	The largest step size, 0 to leave it to the engine.	0.0
dt0	float | Any	The first step size, 0 to estimate it.	0.0
max_order	int	The highest integration order, at most 5.	5
max_steps	int	The step budget, which also sizes the output buffers.	10000
max_events	int	The capacity of the event log, past which event_log.overflow is set.	16
event	Any	An event or event set to monitor.	None
control_state0	Any	The initial control values, {name: value} for a Recipe with bound controls.	None
theta	Any	The parameters, or None for the problem's defaults.	None
consistent_ic	Any	True or False to decide the consistent initialization of this solve, None for the solver's default.	None
lambda0	Any	The initial values of the global unknowns, when y0 holds only the field block.	None
adjoint	Any	The differentiation strategy, or None for the solver's.	None
tolerances	Any	A ToleranceSpec, overriding the solver's scalar tolerances.	None
trace	bool	True attaches a per-step diagnostic trace.	False
history	Any	True records Solution.history, and an earlier solve's Solution.history continues from its end.	None

Returns:

Type	Description
Solution	The Solution.

Raises:

Type	Description
ValueError	If the problem needs a capability the engine does not have.
ConfigurationError	If the engine cannot honour a setting, such as algebraic_error_control="exclude" off JIDA.

Settings shared by JNDF, JBDF, and JIDA

The three engines take the same keyword settings.

The keyword settings of the solver constructors and of stock() and flash().

Attributes:

Name	Type	Description
rtol	float	The relative integration tolerance, default 1e-6.
atol	float	The absolute integration tolerance, default 1e-6.
newton_tol_factor	float	The scale of the Newton acceptance threshold, default 1.0, while the step error test keeps the unscaled tolerances.
ic_max_newton_iters	int	The Newton iteration cap of the consistent initialization, default 10, on JNDF and JBDF only.
linear	Any	The linear backend, a name ("auto", "cudss", "klu", "superlu", "dense") or a ds.linear object.
precision	Precision	The Precision, used only when the problem declares none, so leave it unset.
jacobian_values_dtype	Any	Leave it unset, and set Precision.jacobian_values_dtype on the problem instead.
profile	ExecutionProfile	The ExecutionProfile, by default the class's stock profile.

profile class-attribute instance-attribute

profile: ExecutionProfile = _stock_profile('qndf')

ExecutionProfile dataclass

The step control, corrector and Jacobian policy a solver runs.

stock() and flash() on each solver class build the two presets.

Attributes:

Name	Type	Description
method_source	Literal['qndf', 'bdf', 'ida']	The engine, "qndf", "bdf" or "ida".
controller_mode	Literal['source', 'flash']	"source" runs the reference method's step and order control, and "flash" resizes the step on every accepted step.
nonlinear_mode	Literal['adaptive', 'fixed']	"adaptive" runs the adaptive corrector loop, and "fixed" a bounded number of iterations.
jacobian_policy	Literal['source', 'fresh_each_iteration', 'fresh_each_step']	"source" reuses the Jacobian as the reference method does, and "fresh_each_iteration" or "fresh_each_step" refreshes it.
controller_law	Literal['pi', 'h211b']	The flash step-size law, "pi" (default) or "h211b", read only with controller_mode="flash".

---

Explicit Diffrax solvers

DiffraxStepper

An explicit Diffrax solver run on a problem with identity mass.

Every differential row needs unit mass, and the only other rows allowed are Dirichlet pins. Events can only terminate the solve. It needs the optional diffrax package.

Attributes:

Name	Type	Description
inner	Optional[Any]	Diffrax solver instance, diffrax.Tsit5() when None.
rtol	float	Relative tolerance of the default diffrax.PIDController.
atol	float	Absolute tolerance of the default diffrax.PIDController.
dt0	Optional[float]	Initial step size, chosen by the step size controller when None.
max_steps	int	Step budget of a solve.
adjoint	Optional[AbstractAdjoint]	Differentiation strategy, mapped to the matching Diffrax adjoint, with None meaning CheckpointAdjoint.

inner class-attribute instance-attribute

inner: Optional[Any] = eqx.field(default=None)

rtol class-attribute instance-attribute

rtol: float = eqx.field(default=1e-06, static=True)

atol class-attribute instance-attribute

atol: float = eqx.field(default=1e-06, static=True)

dt0 class-attribute instance-attribute

dt0: Optional[float] = eqx.field(default=None, static=True)

max_steps class-attribute instance-attribute

max_steps: int = eqx.field(default=10000, static=True)

adjoint class-attribute instance-attribute

adjoint: Optional[AbstractAdjoint] = eqx.field(default=None)

solve

solve(problem: Any, t_span: tuple[float, float], y0: Array, *, saveat: Optional[SaveAt] = None, adjoint: Any = None, theta: Any = None, lambda0: Any = None, trace: bool = False, event: Any = None, control_state0: Any = None, consistent_ic: bool = True, dt0: Optional[float] = None, max_steps: Optional[int] = None, tolerances: Any = None, linear: Any = None, **kwargs: Any) -> Solution

Solve problem over t_span from y0 with the wrapped Diffrax solver.

Options that need JNDF, JBDF, or JIDA raise ConfigurationError, among them SaveAt(dense=True) and an integer SaveAt(steps=...).

Parameters:

Name	Type	Description	Default
problem	Any	Problem as normalized by pdesolve.	required
t_span	tuple[float, float]	Pair (t0, t1).	required
y0	Array	Initial state.	required
saveat	Optional[SaveAt]	Output plan, a SaveAt.	None
adjoint	Any	Differentiation strategy, overriding the stepper's own.	None
theta	Any	Parameter pytree, the problem's defaults when None.	None
lambda0	Any	Not supported, and must stay None.	None
trace	bool	Not supported, and must stay False.	False
event	Any	Terminating Event or EventSet, whose continuous events need direction="any".	None
control_state0	Any	Not supported, and must stay None.	None
consistent_ic	bool	If True, move Dirichlet-pinned rows of y0 onto their boundary values first.	True
dt0	Optional[float]	Initial step size, overriding the stepper's.	None
max_steps	Optional[int]	Step budget, overriding the stepper's.	None
tolerances	Any	Not supported, since the stepper's rtol and atol apply.	None
linear	Any	Not supported, and must stay None.	None
**kwargs	Any	Passed to diffrax.diffeqsolve, for example stepsize_controller=.	{}

Returns:

Type	Description
Solution	A Solution, whose result reports a failed solve instead of raising.

---

Tolerances

ToleranceSpec is the per-solve tolerance request of JNDF, JBDF, and JIDA, and ResolvedTolerances is the record a solve returns on Solution.tolerances.

diff_soup.ToleranceSpec dataclass

Per-solve tolerance request for JNDF, JBDF, and JIDA, resolved row by row against the problem.

Pass it as tolerances= to a solve, and read the resolved record from Solution.tolerances. In an atol mapping the most specific key wins, and every row must be covered.

Parameters:

Name	Type	Description	Default
rtol	float	Positive relative tolerance.	1e-06
atol	Any	Positive absolute tolerance, a scalar, one value per row, or a mapping keyed by field, global, role, partition, or "default".	1e-06
ic_rtol	Optional[float]	Relative tolerance of consistent initialization, rtol when None.	None
ic_atol	Any	Absolute tolerance of consistent initialization, in the forms of atol, which it copies when None.	None
newton_tol_factor	float	Positive scale on the corrector's acceptance threshold, where values above one loosen it.	1.0
event_tol	Optional[float]	Positive tolerance stored on the resolved record and not read by the engines.	None
algebraic_error_control	str	"include" keeps algebraic rows in the local error norm, and "exclude", on JIDA only, drops them.	'include'

diff_soup.ResolvedTolerances

Tolerances resolved to one value per row of a problem's state, as carried on Solution.tolerances.

It can be passed back as tolerances= for a problem with the same number of rows. to_dict and from_dict convert it to and from plain data.

Attributes:

Name	Type	Description
rtol	float	Relative tolerance.
atol	ndarray	Absolute tolerance per row.
ic_rtol	float	Relative tolerance of consistent initialization.
ic_atol	ndarray	Absolute tolerance per row for consistent initialization.
newton_tol_factor	float	Scale on the corrector's acceptance threshold.
event_tol	Optional[float]	Event tolerance from the request, or None.
algebraic_error_control	str	"include" or "exclude".
layout_fingerprint	str	Hash of the problem's row layout.
provenance	Tuple[str, ...]	For each field and global, the rule that set its atol.