如何简化梯形法则积分的Numpy实现？求解数组索引报错

一、简化梯形法则积分代码 #

你的原始循环版本已经正确实现了梯形法则计算$\int_a^b \frac{\sin x}{x} dx$，现在用Numpy向量化改写其实很直接，核心是利用Numpy的数组运算替代循环，同时紧扣梯形法则的权重逻辑：

梯形法则的公式可以拆解为：
$$\Delta x \times \left( \frac{f(x_0)}{2} + \sum_{i=1}^{n-1} f(x_i) + \frac{f(x_n)}{2} \right)$$

对应Numpy实现的话，推荐用np.linspace生成采样点（比np.arange更精准，能确保最后一个点刚好是b），然后向量化计算所有点的$f(x)$值，再按权重求和：

import numpy as np

def integral(a, b, n):
    x = np.linspace(a, b, n+1)  # 生成n+1个点，对应n个梯形区间
    f = np.sin(x) / x
    delta = (b - a) / n
    # 按梯形法则加权求和
    return delta * (0.5*f[0] + np.sum(f[1:-1]) + 0.5*f[-1])

如果想写成你期望的合并形式，也可以调整为：

return delta * (f[0] + 2*np.sum(f[1:-1]) + f[-1]) / 2

这和上面的代码逻辑完全一致，只是把0.5提出来统一除以2，更贴合你想要的写法。

二、解决too many indices for array错误 #

你这段代码里的问题出在数组索引的误用：

def integral(a,b,n):
    d = (b-a)/float(n)
    x = np.arange(a,b,d)  # x是一维数组，shape为(n,)
    J = np.where(x[:,1] < np.sin(x[:,0])/x[:,0])[0]

np.arange生成的x是一维数组，而x[:,1]、x[:,0]是二维数组的索引方式（表示取所有行的第1列/第0列），一维数组只有一个轴，自然会报“索引过多”的错误。

如果你是想对每个x元素判断$\sin(x)/x$和某个值的大小，直接用一维数组的运算即可，比如：

# 示例：找所有满足 x < sin(x)/x 的元素索引
f = np.sin(x)/x
J = np.where(x < f)[0]

这样就不会有索引错误了。

内容的提问来源于stack exchange，提问作者JKM