IndexError when integrating the solution of an ODE solved with ODEint

Question

I used odeint to solve a coupled differential equation in terms of v(t) and u(t).

f = odeint(ODEs, f0, t)
    
v = f[:,0]    
u = f[:,1]    

I want to get x, which is the integral of v(t), so I have this code that uses scipy.integrate.quad:

x = lambda t: f[t,0]

delta_x=quad(x, 0,1)
print(delta_x)

When I run this, I get

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

and my Jupyter Notebook highlights f[t,0]. Why is there an index error if the function is defined everywhere within the bounds of integration? How can I solve the error or get the integral via another method?

Here is a minimal reproducible example:

import numpy as np
from scipy.integrate import odeint, quad 

def ODEs(f,t):
    v = f[0]
    u = f[1]
    
    dvdt = u-v
    dudt = (v-u)/3
    
    return [dvdt, dudt]

f0 = [0.3, 0.2]

t = np.linspace(0,4,10000)
f = odeint(ODEs, f0, t)
v = f[:,0]    
u = f[:,1]  

x = lambda t: f[t,0]

delta_x=quad(x, 0,1)
print(delta_x)

Answer 1

As mentioned in the comments, quad will evaluate a function at various points and compute the integral that way. Since you have discrete values from odeint, you cannot use quad. Instead, you can use np.trapz to perform the integration using the trapezoid rule from the discrete results.

import numpy as np
from scipy.integrate import odeint

def ODEs(f,t):
    v, u = f
    dvdt = u - v
    dudt = -dvdt/3
    
    return [dvdt, dudt]

f0 = [0.3, 0.2]

# changed t to the range you want the later integral over
t = np.linspace(0, 1, 1000)
f = odeint(ODEs, f0, t)
v, u = f.T

integral = np.trapz(v, t)
print(integral) # 0.266422665699249

Alternatively, you can use scipy.integrate.solve_ivp with dense_output=True. The ODE integration result will then have a sol attribute which is a function that can be evaluated. This can then be used for quad.

import numpy as np
from scipy.integrate import solve_ivp, quad

# Note solve_ivp by default expect 't' to be first
def ODEs(t, f):
    v, u = f
    dvdt = u - v
    dudt = -dvdt/3
    
    return [dvdt, dudt]

f0 = [0.3, 0.2]

t_span = [0, 4]
f = solve_ivp(ODEs, t_span, f0, dense_output=True)
integral = quad(lambda t: f.sol(t)[0], 0, 1)
# the first term is the integral result and the second is the error
print(integral)  # (0.2663479453362735, 4.315387022498807e-09)