Choosing an Algorithm

CuBIE ships several integration algorithm families. This page helps you pick the right one. Pass the algorithm name as method= to solve_ivp() or algorithm= to Solver; names are case-insensitive.

Decision Guide

Problem type	Recommended family	Notes
Non-stiff	ERK	Fast per step; `dormand-prince-54` is a good default.
Mildly stiff	DIRK or Rosenbrock-W	DIRK is robust; Rosenbrock-W avoids Newton iteration.
Very stiff	FIRK	`radau_iia_5` handles extreme stiffness well.
Fixed step required	`euler` or any non-adaptive tableau	Forward Euler for explicit, `backwards_euler` for implicit.

The bare family names erk, dirk, firk, and rosenbrock are also accepted and select each family’s default tableau (dormand-prince-54, l_stable_dirk_3, firk_gauss_legendre_2, and ros3p respectively).

Adaptive or fixed?

“Adaptive: Yes” below means the method produces an embedded error estimate, so it can drive an adaptive step-size controller. Methods without an error estimate always run at a fixed step: if you pair one with an adaptive controller, CuBIE issues a UserWarning and silently swaps in the fixed-step controller, because there is no error signal to adapt on.

Available Algorithms

Explicit Runge–Kutta (ERK)

Name	Order	Adaptive	Notes
`heun-21`	2	No	Heun’s method.
`ralston-33`	3	No	Ralston’s third-order method.
`bogacki-shampine-32` / `rk23` / `ode23`	3(2)	Yes	Low-order, cheap.
`classical-rk4` / `rk4`	4	No	The classical Runge–Kutta method.
`fehlberg-45`	5(4)	Yes	Fehlberg’s method.
`cash-karp-54`	5(4)	Yes	Cash–Karp coefficients; defaults to Gustafsson control.
`dormand-prince-54` / `dopri54` / `rk45` / `ode45`	5(4)	Yes	Industry standard; good default and the ERK family default.
`tsit5`	5(4)	Yes	Tsitouras 5(4); often slightly more efficient than Dormand–Prince.
`vern7`	7(6)	Yes	Verner’s high-order method; defaults to Gustafsson control.
`dormand-prince-853` / `dop853`	8(5,3)	Yes	High order for smooth problems; defaults to Gustafsson control.

Diagonally Implicit Runge–Kutta (DIRK)

Name	Order	Adaptive	Notes
`implicit_midpoint`	2	No	Symmetric, energy-preserving.
`trapezoidal_dirk` / `ode23t`	2	No	Trapezoidal rule.
`sdirk_2_2`	2	No	L-stable SDIRK; has no embedded error estimate.
`l_stable_dirk_3`	3	No	Default DIRK tableau; L-stable, stiffly accurate, 3 stages.
`kvaerno3`	3	Yes	A-L stable ESDIRK, 4 stages.
`kvaerno5`	5	Yes	A-L stable ESDIRK, 7 stages; defaults to exact Newton.
`l_stable_sdirk_4`	4	Yes	L-stable, 5 stages; defaults to exact Newton.
`eldirk32_euler`	2	Yes	ELDIRK: implicit Euler, two explicit stages, third-order estimate.
`eldirk32_trapezoidal`	2	Yes	ELDIRK: trapezoidal rule, one explicit stage, third-order estimate.
`eldirk32_ellsiepen`	2	Yes	ELDIRK: `sdirk_2_2`, one explicit stage, third-order estimate.

Fully Implicit Runge–Kutta (FIRK)

Name	Order	Adaptive	Notes
`firk_gauss_legendre_2`	4	No	2-stage Gauss–Legendre; default FIRK.
`firk_gauss_legendre_4`	8	Yes	4-stage Gauss–Legendre; conservative step control from a second-order error estimate.
`radau_iia_3`	3	Yes	2-stage Radau IIA; the cheapest of the family.
`radau_iia_5` / `radau`	5	Yes	3-stage Radau IIA; excellent for stiff problems.
`radau_iia_9`	9	Yes	5-stage Radau IIA; for tight tolerances, five coupled stages per step. Defaults to BiCGSTAB with exact Newton.

Rosenbrock-W

Name	Order	Adaptive	Notes
`ros3p`	3	Yes	Default Rosenbrock; linearly implicit, no Newton iteration.
`rodas3p`	3	Yes	Stiffly accurate variant.
`rosenbrock23_sciml` / `rosenbrock23` / `ode23s`	2(3)	Yes	MATLAB `ode23s`-equivalent / SciML-compatible tableau.

Simple methods

Name	Order	Adaptive	Notes
`euler`	1	No	Forward Euler; explicit.
`backwards_euler`	1	No	Backward Euler; implicit, L-stable.
`backwards_euler_pc`	1	No	Predictor-corrector backward Euler.
`crank_nicolson`	2	Yes	Implicit trapezoidal rule; a second backward-Euler solve provides the embedded error estimate.

Choosing a Controller

Adaptive algorithms use an error controller to adjust the step size. Select one by name with step_controller inside step_control_settings; the registered names are "fixed", "i", "pi", "pid", and "gustafsson".

fixed

Constant step size; no error control. Used automatically for methods without an error estimate.

i, pi, pid

The i, pi, and pid controllers are all the same PID controller, but the i and pi device functions compile to more efficient code when you don’t need the additional gain terms. After every step they compare the estimated error with the tolerance (atol + rtol * max(|this_state|, |last_state|)) and propose a next-step size to bring the ratio of error/tolerance to 1. Each gain term acts on the current error and a subset of previous errors:

integral_gain multiplies by the error/tolerance ratio of the current step. It sets how much this error adds to the running correction.
proportional_gain multiplies by how much the ratio changed since the previous step. It sets how much the running correction is adjusted by the change in error.
derivative_gain multiplies by how much the change in the ratio changed over two steps. It acts against large changes in the rate of change of error, damping the controller.

Passing gains without any step_controller selected will select the fastest of the i, pi, pid controllers that can provide your gains, or you can provide a string to select a common literature controller using the filter_coefficients argument.

gustafsson

Predictive controller that accounts for the previous step’s error ratio. Widely used with implicit methods; useful when step rejections are frequent.

Each algorithm family picks a default controller (i for ERK, pi for DIRK, gustafsson for FIRK and Rosenbrock, fixed without an error estimate); some tableaus override it. To override:

solver = qb.Solver(
    LV,
    algorithm="dormand-prince-54",
    step_control_settings={"step_controller": "gustafsson"},
)

Automatic calibration

The fastest algorithm order, linear solver, preconditioner, and Newton variant depend on the system. Given a representative batch, the solver can measure them directly:

solver = qb.Solver(LV, algorithm="ros3p", atol=1e-6, rtol=1e-6)
report = solver.calibrate(
    {"x": x0_values},
    {"alpha": alpha_values},
    duration=10.0,
)
print(report.summary())

Solver.calibrate races different solver configurations on your problem: it compares a few orders of each algorithm family and, for the implicit families, the preconditioner, linear-solver, Newton-variant, smoothed-error, and dense-predictor settings. Candidates that fail to integrate the grid are dropped before timing; survivors are ranked on a few full-length solves, and the leaders are returned with solve times and failure counts. By default the winning configuration is applied to the solver in place; pass apply=False to get the race results without modifying the solver.

With the algorithm settled, Solver.optimize tries a few different configurations of where buffers sit in memory, which loops get unrolled, and how many threads run at once on the GPU to find the fastest configuration for your system:

report = solver.optimize(
    {"x": x0_values},
    {"alpha": alpha_values},
    duration=10.0,
)
print(report.summary())

Settings you gave, or an earlier optimize applied, stay fixed unless force=True. The fastest configuration is applied in place; apply=False only reports.

By default optimize times a small batch at a short duration. Pass auto_size=False only when you will run one grid size at one duration repeatedly; it spends longer optimising that exact solve.

For the mathematical background behind these algorithms, see Numerical Integration.