pulsesuite.libpulsesuite.integrator

ODE integrators ported from integrator.F90.

1:1 port of the Fortran integrator module. Every public Fortran subroutine has a Python callable with the same name. Fortran interfaces that dispatch on type/rank are kept as separate named functions.

Vectorised with NumPy; no Numba or guardrails required.

Implements:

  • Cash-Karp Runge-Kutta step (rkck)

  • Quality-controlled adaptive RK step (rkqs)

  • Classic 4th-order Runge-Kutta (rk4)

  • Forward Euler (idiot)

  • Adaptive ODE driver (odeint)

  • Simple fixed-step driver (simpleint)

  • Modified midpoint (mmid)

  • Polynomial extrapolation (pzextr)

  • Bulirsch-Stoer adaptive step (bsstep)

  • Semi-implicit midpoint for stiff systems (simpr)

  • Stiff Bulirsch-Stoer (stifbs)

  • LU decomposition with partial pivoting (ludcmp)

  • LU back-substitution (lubksb)

Author: Rahul R. Sah

Attributes

sp

dp

SAFETY

PGROW

PSHRNK

ERRCON

TINY

MAXSTP

KMAXX

Functions

rkck_dp(y, dydt, t, h, D)

Cash-Karp RK step for real (float64) ODEs. Returns (yout, yerr).

rkck_dpc(y, dydt, t, h, D)

Cash-Karp RK step for complex (complex128) ODEs. Returns (yout, yerr).

rkck_3D_dpc(y, dydt, t, h, D)

Cash-Karp RK step for 3D complex arrays. Returns (yout, yerr).

rk4_dp(y, dydt, t, h, D)

Classic RK4 for real ODEs. Modifies y in-place. Returns t_new.

rk4_dpc(y, dydt, t, h, D)

Classic RK4 for complex ODEs. Modifies y in-place. Returns t_new.

rk4_3D_dpc(y, dydt, t, h, D)

Classic RK4 for 3D complex arrays. Modifies y in-place. Returns t_new.

idiot_dp(y, dydt, t, h, D)

Forward Euler for real ODEs. Modifies y in-place. Returns t_new.

idiot_dpc(y, dydt, t, h, D)

Forward Euler for complex ODEs. Modifies y in-place. Returns t_new.

rkqs_dp(y, dydt, t, htry, eps, yscale, D, jacobn)

Quality-controlled RK step for real ODEs.

rkqs_dpc(y, dydt, t, htry, eps, yscale, D, jacobn)

Quality-controlled RK step for complex ODEs.

rkqs_3D_dpc(y, dydt, t, htry, eps, yscale, D)

Quality-controlled RK step for 3D complex arrays (UPPE, no Jacobian).

rkqs_3D_dpc_TOM(y, dydt, t, htry, eps, yscaleREAL, ...)

Quality-controlled RK step for 3D complex arrays with separate

odeint_dp(y, x1, x2, eps, h1, hmin, D, integ, jacobn)

Adaptive ODE driver for real ODEs.

odeint_dpc(y, x1, x2, eps, h1, hmin, D, integ, jacobn)

Adaptive ODE driver for complex ODEs.

odeint_3D_dpc(y, x1, x2, eps, h1, hmin, D, integ)

Adaptive ODE driver for 3D complex arrays (UPPE, no Jacobian).

odeint_3D_dpc_TOM(y, x1, x2, eps, h1, hmin, D, integ)

Adaptive ODE driver for 3D complex arrays with separate real/imaginary

simpleint_dp(y, t1, t2, h1, D, simple)

Simple fixed-step driver for real ODEs. Y1 = DY * DX + Y0.

simpleint_dpc(y, t1, t2, h1, D, simple)

Simple fixed-step driver for complex ODEs.

mmid(y, dydx, xs, htot, nstep, D)

Modified midpoint method.

pzextr(iest, xest, yest, yz, dy, x_save, d_save)

Polynomial extrapolation step (Neville's algorithm).

bsstep(y, dydx, x, htry, eps, yscale, D, jacobn)

Bulirsch-Stoer adaptive step using modified midpoint (mmid) and

ludcmp(a, indx)

LU decomposition with partial pivoting (Crout's algorithm).

lubksb(a, indx, b)

LU back-substitution.

simpr(y, dydx, dfdx, dfdy, xs, htot, nstep, derivs)

Semi-implicit midpoint rule for stiff ODEs.

stifbs(y, dydx, x, htry, eps, yscal, derivs, jacobn)

Single step of the stiff Bulirsch-Stoer algorithm.

Module Contents

pulsesuite.libpulsesuite.integrator.sp
pulsesuite.libpulsesuite.integrator.dp
pulsesuite.libpulsesuite.integrator.SAFETY = 0.9
pulsesuite.libpulsesuite.integrator.PGROW = -0.2
pulsesuite.libpulsesuite.integrator.PSHRNK = -0.25
pulsesuite.libpulsesuite.integrator.ERRCON = 0.00018895680000000003
pulsesuite.libpulsesuite.integrator.TINY = 1e-30
pulsesuite.libpulsesuite.integrator.MAXSTP = 100000000
pulsesuite.libpulsesuite.integrator.KMAXX = 8
pulsesuite.libpulsesuite.integrator.rkck_dp(y, dydt, t, h, D)

Cash-Karp RK step for real (float64) ODEs. Returns (yout, yerr).

pulsesuite.libpulsesuite.integrator.rkck_dpc(y, dydt, t, h, D)

Cash-Karp RK step for complex (complex128) ODEs. Returns (yout, yerr).

pulsesuite.libpulsesuite.integrator.rkck_3D_dpc(y, dydt, t, h, D)

Cash-Karp RK step for 3D complex arrays. Returns (yout, yerr).

pulsesuite.libpulsesuite.integrator.rk4_dp(y, dydt, t, h, D)

Classic RK4 for real ODEs. Modifies y in-place. Returns t_new.

pulsesuite.libpulsesuite.integrator.rk4_dpc(y, dydt, t, h, D)

Classic RK4 for complex ODEs. Modifies y in-place. Returns t_new.

pulsesuite.libpulsesuite.integrator.rk4_3D_dpc(y, dydt, t, h, D)

Classic RK4 for 3D complex arrays. Modifies y in-place. Returns t_new.

pulsesuite.libpulsesuite.integrator.idiot_dp(y, dydt, t, h, D)

Forward Euler for real ODEs. Modifies y in-place. Returns t_new.

pulsesuite.libpulsesuite.integrator.idiot_dpc(y, dydt, t, h, D)

Forward Euler for complex ODEs. Modifies y in-place. Returns t_new.

pulsesuite.libpulsesuite.integrator.rkqs_dp(y, dydt, t, htry, eps, yscale, D, jacobn)

Quality-controlled RK step for real ODEs.

Modifies y in-place. Returns (t_new, hdid, hnext).

pulsesuite.libpulsesuite.integrator.rkqs_dpc(y, dydt, t, htry, eps, yscale, D, jacobn)

Quality-controlled RK step for complex ODEs.

Modifies y in-place. Returns (t_new, hdid, hnext).

pulsesuite.libpulsesuite.integrator.rkqs_3D_dpc(y, dydt, t, htry, eps, yscale, D)

Quality-controlled RK step for 3D complex arrays (UPPE, no Jacobian).

Modifies y in-place. Returns (t_new, hdid, hnext).

pulsesuite.libpulsesuite.integrator.rkqs_3D_dpc_TOM(y, dydt, t, htry, eps, yscaleREAL, yscaleIMAG, D)

Quality-controlled RK step for 3D complex arrays with separate real/imaginary error tolerances. eps is complex: real(eps) for real part tolerance, imag(eps) for imaginary part tolerance.

Modifies y in-place. Returns (t_new, hdid, hnext).

pulsesuite.libpulsesuite.integrator.odeint_dp(y, x1, x2, eps, h1, hmin, D, integ, jacobn)

Adaptive ODE driver for real ODEs.

Parameters:

  • y (ndarray, float64 — state vector (modified in-place))

  • x1 (float — integration bounds)

  • x2 (float — integration bounds)

  • eps (float — error tolerance)

  • h1 (float — initial step size)

  • hmin (float — minimum step size)

  • D (callable(t, y) -> dy — derivative function)

  • integ (callable — stepper (rkqs_dp, bsstep, stifbs))

  • jacobn (callable — Jacobian function)

Returns:

  • nok (int — number of good steps)

  • nbad (int — number of bad (retried) steps)

pulsesuite.libpulsesuite.integrator.odeint_dpc(y, x1, x2, eps, h1, hmin, D, integ, jacobn)

Adaptive ODE driver for complex ODEs.

Same interface as odeint_dp but y is complex128. Returns (nok, nbad).

pulsesuite.libpulsesuite.integrator.odeint_3D_dpc(y, x1, x2, eps, h1, hmin, D, integ)

Adaptive ODE driver for 3D complex arrays (UPPE, no Jacobian).

h1 is returned modified (remembers last used step size for the bridge between pulsesuite “Steps”).

Returns (nok, nbad, h1_out).

pulsesuite.libpulsesuite.integrator.odeint_3D_dpc_TOM(y, x1, x2, eps, h1, hmin, D, integ)

Adaptive ODE driver for 3D complex arrays with separate real/imaginary error tolerances. eps is complex.

Returns (nok, nbad, h1_out).

pulsesuite.libpulsesuite.integrator.simpleint_dp(y, t1, t2, h1, D, simple)

Simple fixed-step driver for real ODEs. Y1 = DY * DX + Y0.

Parameters:

  • y (ndarray, float64 (modified in-place))

  • t1 (float — integration bounds)

  • t2 (float — integration bounds)

  • h1 (float — step size)

  • D (callable(t, y) -> dy)

  • simple (callable — single-step method (rk4_dp, idiot_dp))

Returns:

nok

Return type:

int — number of steps taken

pulsesuite.libpulsesuite.integrator.simpleint_dpc(y, t1, t2, h1, D, simple)

Simple fixed-step driver for complex ODEs.

Returns nok.

pulsesuite.libpulsesuite.integrator.mmid(y, dydx, xs, htot, nstep, D)

Modified midpoint method.

Takes nstep sub-steps of size h = htot/nstep using the modified midpoint (Stoermer) rule. Does NOT modify y.

Returns yout.

pulsesuite.libpulsesuite.integrator.pzextr(iest, xest, yest, yz, dy, x_save, d_save)

Polynomial extrapolation step (Neville’s algorithm).

Parameters:

  • iest (int — 1-based step index)

  • xest (float — x-estimate (typically (h/nseq)**2))

  • yest (ndarray — y-estimate from midpoint/simpr)

  • yz (ndarray — output y (modified in-place))

  • dy (ndarray — output error (modified in-place))

  • x_save (ndarray, shape (KMAXX,) — persistent x storage)

  • d_save (ndarray, shape (nv, KMAXX) — persistent d storage)

pulsesuite.libpulsesuite.integrator.bsstep(y, dydx, x, htry, eps, yscale, D, jacobn)

Bulirsch-Stoer adaptive step using modified midpoint (mmid) and polynomial extrapolation (pzextr).

Modifies y in-place. Returns (t_new, hdid, hnext).

pulsesuite.libpulsesuite.integrator.ludcmp(a, indx)

LU decomposition with partial pivoting (Crout’s algorithm).

Modifies a in-place (overwritten with L and U factors). Fills indx with pivot indices (0-based). Returns d (+1.0 or -1.0 for even/odd row exchanges).

pulsesuite.libpulsesuite.integrator.lubksb(a, indx, b)

LU back-substitution.

Solves A*x = b given the LU factors in a and pivot indices in indx. Modifies b in-place with the solution.

pulsesuite.libpulsesuite.integrator.simpr(y, dydx, dfdx, dfdy, xs, htot, nstep, derivs)

Semi-implicit midpoint rule for stiff ODEs.

Uses LU decomposition of (I - h*dfdy) and performs nstep sub-steps. Does NOT modify y.

Returns yout.

pulsesuite.libpulsesuite.integrator.stifbs(y, dydx, x, htry, eps, yscal, derivs, jacobn)

Single step of the stiff Bulirsch-Stoer algorithm.

Uses semi-implicit midpoint (simpr) with polynomial extrapolation (pzextr).

Modifies y in-place. Returns (t_new, hdid, hnext).