How to solve an ODE initial value problem with SciPy

Time-dependent models often start from a known state and a rule for how that state changes. SciPy handles that form as an initial value problem, where scipy.integrate.solve_ivp() advances the state from the starting time to a chosen ending time.

The solve_ivp() call takes a derivative function with signature fun(t, y), an integration interval, an initial state vector, and optional sample times. The returned object stores solver status, the time grid, state values, and counters such as the number of derivative evaluations.

A scalar exponential decay equation keeps the check easy to audit because its analytic solution is known. Matching the numerical values against 2 * exp(-0.4t) confirms that the solver reached the requested time points and that the state array is being read in the right orientation.

1. Create a Python script named ode_initial_value_demo.py.

   ode_initial_value_demo.py

   import numpy as np
   from scipy.integrate import solve_ivp
    
    
   decay_rate = 0.4
   initial_value = 2.0
   sample_times = np.array([0.0, 1.0, 2.5, 5.0, 10.0])
    
    
   def exponential_decay(t, y):
       return [-decay_rate * y[0]]
    
    
   solution = solve_ivp(
       exponential_decay,
       t_span=(sample_times[0], sample_times[-1]),
       y0=[initial_value],
       t_eval=sample_times,
       method="RK45",
       rtol=1e-7,
       atol=1e-10,
   )
    
   expected = initial_value * np.exp(-decay_rate * sample_times)
   max_abs_error = np.max(np.abs(solution.y[0] - expected))
    
   print(f"success: {solution.success}")
   print(f"message: {solution.message}")
   print(f"t: {np.array2string(solution.t, precision=2)}")
   print(f"y: {np.array2string(solution.y[0], precision=8)}")
   print(f"expected: {np.array2string(expected, precision=8)}")
   print(f"max_abs_error: {max_abs_error:.2e}")
   print(f"function_calls: {solution.nfev}")

   t_eval requests solution values at specific times. Without it, solve_ivp() returns the solver's internally chosen time grid.

2. Run the ODE initial value script.

   $ python ode_initial_value_demo.py
   success: True
   message: The solver successfully reached the end of the integration interval.
   t: [ 0.   1.   2.5  5.  10. ]
   y: [2.         1.34064009 0.73575889 0.27067058 0.03663128]
   expected: [2.         1.34064009 0.73575888 0.27067057 0.03663128]
   max_abs_error: 1.27e-08
   function_calls: 176

3. Confirm that the solver reached the final time.

   success should be True, and message should say the integration interval was reached before the state values are used downstream.

4. Read the sampled time and state arrays.

   solution.t contains the requested sample times. solution.y stores one row per state variable and one column per time value, so the scalar state is read as solution.y[0].

5. Compare the numerical values with the known solution.

   The printed max_abs_error is below 1e-6 for this demonstration, which confirms the returned values match 2 * exp(-0.4t) at the requested sample times.

6. Replace the derivative function and initial state with the model being simulated.

   RK45 is a common first choice for non-stiff equations. Compare Radau, BDF, or LSODA when small tolerance changes cause many more function calls or the solution changes unexpectedly.

7. Remove the demo script when it was only created for the check.

   $ rm ode_initial_value_demo.py