cds2.integrate

Numerical integration and ODE solving on scipy.integrate.

Quadrature and ODE — convenience re-export of scipy.integrate; see cds2.sde for SDE companion.

QuadResult dataclass

Integral estimate with its absolute error bound.

Source code in src/cds2/integrate.py

@dataclass(frozen=True)
class QuadResult:
    """Integral estimate with its absolute error bound."""

    value: float
    error: float

OdeResult dataclass

Trajectory produced by an initial-value-problem solver.

Source code in src/cds2/integrate.py

@dataclass(frozen=True)
class OdeResult:
    """Trajectory produced by an initial-value-problem solver."""

    t: np.ndarray
    y: np.ndarray
    success: bool
    message: str
    t_events: tuple[np.ndarray, ...] | None = None
    y_events: tuple[np.ndarray, ...] | None = None

BVPResult dataclass

Solution of a two-point boundary value problem on the output mesh.

Source code in src/cds2/integrate.py

@dataclass(frozen=True)
class BVPResult:
    """Solution of a two-point boundary value problem on the output mesh."""

    x: np.ndarray
    y: np.ndarray
    success: bool
    message: str

solve_bvp(func, bc, initial_mesh, initial_guess, tol=0.001, max_nodes=1000)

Solve a two-point boundary value problem y' = func(t, y), bc(ya, yb) = 0.

initial_mesh is the starting node grid; initial_guess is a (n_states, len(mesh)) array of starting values for the collocation solver (4th-order collocation with residual control).

Source code in src/cds2/integrate.py

def solve_bvp(
    func: Callable[[np.ndarray, np.ndarray], np.ndarray],
    bc: Callable[[np.ndarray, np.ndarray], np.ndarray],
    initial_mesh: Sequence[float],
    initial_guess: object,
    tol: float = 1e-3,
    max_nodes: int = 1000,
) -> BVPResult:
    """Solve a two-point boundary value problem y' = func(t, y), bc(ya, yb) = 0.

    ``initial_mesh`` is the starting node grid; ``initial_guess`` is a
    (n_states, len(mesh)) array of starting values for the collocation
    solver (4th-order collocation with residual control).
    """
    res = spi.solve_bvp(
        func,
        bc,
        np.asarray(initial_mesh, dtype=float),
        np.atleast_2d(np.asarray(initial_guess, dtype=float)),
        tol=tol,
        max_nodes=max_nodes,
    )
    return BVPResult(
        x=np.asarray(res.x),
        y=np.asarray(res.y),
        success=bool(res.success),
        message=str(res.message),
    )

quad(func, a, b)

Adaptive definite integral of a scalar function over [a, b].

Source code in src/cds2/integrate.py

def quad(func: Callable[[float], float], a: float, b: float) -> QuadResult:
    """Adaptive definite integral of a scalar function over [a, b]."""
    value, error = spi.quad(func, a, b)
    return QuadResult(value=value, error=error)

integrate_2d(func, xa, xb, ya, yb)

Double integral of func(x, y) over the rectangle [xa, xb] x [ya, yb].

Source code in src/cds2/integrate.py

def integrate_2d(
    func: Callable[[float, float], float],
    xa: float,
    xb: float,
    ya: float,
    yb: float,
) -> QuadResult:
    """Double integral of ``func(x, y)`` over the rectangle [xa, xb] x [ya, yb]."""
    value, error = spi.dblquad(lambda y, x: func(x, y), xa, xb, ya, yb)
    return QuadResult(value=value, error=error)

integrate_3d(func, xa, xb, ya, yb, za, zb)

Triple integral of func(x, y, z) over the axis-aligned box.

Source code in src/cds2/integrate.py

def integrate_3d(
    func: Callable[[float, float, float], float],
    xa: float,
    xb: float,
    ya: float,
    yb: float,
    za: float,
    zb: float,
) -> QuadResult:
    """Triple integral of ``func(x, y, z)`` over the axis-aligned box."""
    value, error = spi.tplquad(lambda z, y, x: func(x, y, z), xa, xb, ya, yb, za, zb)
    return QuadResult(value=value, error=error)

solve_ivp(func, t_span, y0, t_eval=None, method='RK45', events=None, **kwargs)

Integrate an initial-value problem dy/dt = func(t, y).

Backward integration works naturally when t_span[1] < t_span[0]. Pass events (callables g(t, y) -> float; zero crossings trigger them, set terminal = True on a callable to stop integration) to receive t_events / y_events on the result. For stiff problems use method="Radau", "BDF" or "LSODA".

Source code in src/cds2/integrate.py

def solve_ivp(
    func: Callable[[float, np.ndarray], np.ndarray],
    t_span: tuple[float, float],
    y0: Sequence[float],
    t_eval: Sequence[float] | None = None,
    method: str = "RK45",
    events: Sequence[Callable[[float, np.ndarray], float]] | None = None,
    **kwargs: object,
) -> OdeResult:
    """Integrate an initial-value problem dy/dt = func(t, y).

    Backward integration works naturally when ``t_span[1] < t_span[0]``.
    Pass ``events`` (callables ``g(t, y) -> float``; zero crossings trigger
    them, set ``terminal = True`` on a callable to stop integration) to
    receive ``t_events`` / ``y_events`` on the result. For stiff problems
    use ``method="Radau"``, ``"BDF"`` or ``"LSODA"``.
    """
    res = spi.solve_ivp(
        func,
        tuple(t_span),
        np.asarray(y0, dtype=float),
        t_eval=None if t_eval is None else np.asarray(t_eval, dtype=float),
        method=method,
        events=events,
        **kwargs,
    )
    t_events = None
    y_events = None
    if events is not None:
        t_events = tuple(np.asarray(block) for block in res.t_events)
        y_events = tuple(np.asarray(block) for block in res.y_events)
    return OdeResult(
        t=np.asarray(res.t),
        y=np.asarray(res.y),
        success=bool(res.success),
        message=str(res.message),
        t_events=t_events,
        y_events=y_events,
    )

trapezoid(y, x=None, dx=1.0)

Composite trapezoidal integral of sampled values.

Source code in src/cds2/integrate.py

def trapezoid(y: Sequence[float], x: Sequence[float] | None = None, dx: float = 1.0) -> float:
    """Composite trapezoidal integral of sampled values."""
    return float(spi.trapezoid(np.asarray(y, dtype=float), x=x, dx=dx))

cumulative_trapezoid(y, x=None, dx=1.0)

Running trapezoidal integral, returned with the same length as y.

Source code in src/cds2/integrate.py

def cumulative_trapezoid(
    y: Sequence[float],
    x: Sequence[float] | None = None,
    dx: float = 1.0,
) -> np.ndarray:
    """Running trapezoidal integral, returned with the same length as ``y``."""
    values = spi.cumulative_trapezoid(np.asarray(y, dtype=float), x=x, dx=dx, initial=0.0)
    return np.asarray(values)

simpson(y, x=None, dx=1.0)

Composite Simpson-rule integral of sampled values.

Source code in src/cds2/integrate.py

def simpson(y: Sequence[float], x: Sequence[float] | None = None, dx: float = 1.0) -> float:
    """Composite Simpson-rule integral of sampled values."""
    return float(spi.simpson(np.asarray(y, dtype=float), x=x, dx=dx))