集成使用 ODEint 求解的 ODE 的解时出现索引错误

Question

我使用 odeint 来求解 v(t) 和 u(t) 的耦合微分方程

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

u = f[:,1]    

当我尝试通过编写这段代码来获取 x（即 v(t) 的积分）时。

x = lambda t: f[t,0]

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

我收到“IndexError：只有整数、切片 (

:

)、省略号 (

...

)、numpy.newaxis (

None

) 和整数或布尔数组是有效索引”，Jupyter Notebook 突出显示 f[t,0] 。如果函数在积分范围内的任何地方定义，为什么会出现索引错误？我该如何解决错误或通过其他方法获得积分？

编辑：最小可重现示例：

import numpy as np
from scipy.integrate import odeint
from scipy.integrate import quad
import matplotlib.pyplot as plt

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

正如评论中提到的，

quad

将在各个点评估函数并以这种方式计算积分。由于您有来自

ODEint

的离散值，因此您不能使用

quad

。相反，您可以使用

np.trapz

使用梯形规则根据离散结果执行积分。

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

或者，您可以将

scipy.integrate.solve_ivp

与

dense_output=True

一起使用。 ODE 积分结果将具有

sol

属性，这是一个可以评估的函数。然后可以将其用于

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)
t = np.linspace(*t_span, 1000)
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)