如何避免scipy solve_ivp中LSODA的Python调用开销

Great question—this is a common pain point when dealing with Python-Fortran interface overhead in ODE solvers. The issue you're hitting makes total sense: SciPy's solve_ivp LSODA wrapper is built around step-by-step integration (itask=5) to support features like event handling, adaptive output sampling, and intermediate state checks. Switching directly to itask=1 breaks things because the wrapper isn't set up to handle a single full integration call or process the resulting state correctly.

Here are your most practical solutions, ordered from easiest to most involved:

1. Use SciPy's ode Class (Simplest Fix) #

If you don't need solve_ivp's specific features (like event detection with terminal conditions or the flexible t_eval interface), SciPy's older ode class uses a different LSODA interface that defaults to full integration (itask=1) out of the box. This minimizes Python-Fortran round-trips to just one call, eliminating per-timestep overhead entirely.

Example usage:

from scipy.integrate import ode
import numpy as np

def ode_rhs(t, y):
    # Replace with your actual ODE right-hand side
    return -y  # Simple exponential decay example

# Initial conditions and integration bounds
y0 = np.array([1.0])
t_start = 0.0
t_end = 10.0

# Set up the solver with LSODA
solver = ode(ode_rhs).set_integrator('lsoda')
solver.set_initial_value(y0, t_start)

# Integrate directly to the final time (single Python-Fortran call)
final_state = solver.integrate(t_end)

print(f"Final state at t={t_end}: {final_state}")

If you need intermediate outputs, you can use the set_solout method to attach a callback that collects data during integration—this still reduces overhead drastically compared to per-timestep calls, since the callback only runs at your desired output points, not every internal LSODA timestep.

2. Modify the solve_ivp LSODA Wrapper (Advanced) #

If you rely on solve_ivp's features (like event handling or the t_eval interface), modifying the wrapper is possible but requires updating more than just the itask value. Here's what you'd need to adjust in SciPy's _ivp/lsoda.py:

• Replace the step-by-step loop with a single LSODA call using itask=1
• Update the solver state handling to process the full integration result (instead of a single step)
• Adjust output collection to pull all internal timesteps (or use dense output) from LSODA's workspace
• Add checks for integration errors, convergence, and event triggers after the single call

This is not trivial—SciPy's wrapper is tightly coupled to the step-by-step workflow. Only attempt this if you have deep familiarity with both the wrapper code and ODEPACK's LSODA documentation.

3. Call LSODA Directly (Lowest Level) #

If you want full control, you can call the compiled LSODA Fortran routine directly via SciPy's private _lsoda module. This skips all Python wrapper logic entirely, but requires adhering strictly to ODEPACK's parameter conventions.

Here's a minimal example based on SciPy's internal wrapper code:

from scipy.integrate import _lsoda
import numpy as np

def lsoda_rhs(t, y, dy):
    # LSODA's Fortran interface expects you to write results to `dy`
    dy[:] = -y  # Replace with your ODE right-hand side

n = 1  # Number of state variables
t0 = 0.0
tf = 10.0
y = np.array([1.0], dtype=np.float64)
t = np.array([t0], dtype=np.float64)

# Configure LSODA parameters (refer to ODEPACK docs for details)
itask = np.array([1], dtype=np.int32)  # Integrate to final time
istate = np.array([1], dtype=np.int32)  # Initial state flag
mf = np.array([22], dtype=np.int32)  # Auto-detect stiff/non-stiff
rtol = np.array([1e-6], dtype=np.float64)
atol = np.array([1e-9], dtype=np.float64)

# Allocate workspace arrays
lrw = 20 + 16 * n
liw = 20 + n
rw = np.zeros(lrw, dtype=np.float64)
iw = np.zeros(liw, dtype=np.int32)
iopt = np.array([0], dtype=np.int32)
opt = np.array([], dtype=np.float64)

# Call the raw LSODA routine
_lsoda.lsoda(lsoda_rhs, n, y, t, tf, itask, istate, iopt, rtol, atol,
             0, mf, rw, lrw, iw, liw, opt, 0)

print(f"Final time: {t[0]:.6f}, Final state: {y}")

Key notes for this approach:

• You must follow ODEPACK's parameter rules (e.g., workspace sizes, flag meanings) exactly
• The RHS function must write to the dy array (not return a value)
• Check the istate array after the call to detect integration errors or convergence issues

备注：内容来源于stack exchange，提问作者billbert