python 积分用ODEint求解的ODE的解时的IndexError

xienkqul 于 2024-01-05 发布在 Python

关注(0)|答案(1)|浏览(149)

我用odeint来解一个关于v（t）和u（t）的耦合微分方程。

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

字符串
我想得到x，它是v（t）的积分，所以我有这样一段代码，它使用了scipy.integrate.quad：

x = lambda t: f[t,0]
delta_x=quad(x, 0,1)
print(delta_x)

型
当我运行这个时，我得到

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

型
而我的笔记本突出显示了f[t，0]。如果函数在积分范围内处处定义，为什么会有索引错误？我如何解决错误或通过另一种方法获得积分？
下面是一个最小的可重复的例子：

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)

型

python

来源：https://stackoverflow.com/questions/77736770/indexerror-when-integrating-the-solution-of-an-ode-solved-with-odeint

1条答案

按热度按时间

mzillmmw1#

正如注解中提到的，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)
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)

型

展开查看全部

赞(0）回复(0）举报 2024-01-05

我来回答

python 积分用ODEint求解的ODE的解时的IndexError

1条答案

相关问题

热门标签

最新问答