如何对函数数组积分？解决integrate.quad数组输入报错问题

错误原因 #

scipy.integrate.quad是为标量函数设计的，它期望被积函数返回单个数值。你的代码中被积函数引用了X、Y这类meshgrid生成的数组，导致被积函数返回数组，和quad的标量要求冲突，因此抛出only size-1 arrays can be converted to Python scalars错误。numpy.vectorize本质是伪向量化，底层还是循环，且无法改变quad不支持数组输入的特性，所以无效。

可行方案（除手动循环外） #

方案1：使用scipy.integrate.quad_vec（推荐） #

从scipy 1.7.0版本开始，新增了quad_vec函数，它是quad的向量化版本，直接支持被积函数返回数组，能一次性处理所有X、Y点的积分计算，效率远高于手动循环。

针对你的补充示例，修改后的代码如下：

import numpy as np
from scipy import integrate

N = 3                               # 每个方向的点数
x_start, x_end = -2.5, 2.5            # x方向边界
y_start, y_end = -3.0, 3.0            # y方向边界
x = np.linspace(x_start, x_end, N)    # x的一维数组
y = np.linspace(y_start, y_end, N)    # y的一维数组
X, Y = np.meshgrid(x, y)

def integral_x():
    def integrand(s):
        # s是标量，X、Y是数组，numpy会自动完成广播计算
        return ((2*(X - s))) / (((X - s))**2 + (Y - s)**2)
    # 用quad_vec替代quad，返回的第一个元素就是和X/Y尺寸一致的积分结果数组
    return integrate.quad_vec(integrand, 0.0, 10)[0]

result = integral_x()
print(result.shape)  # 输出 (3, 3)，与X、Y尺寸匹配

对应你最初的带参数函数，修改方式类似：

import numpy as np
import math
from scipy import integrate

def integral_x(p_i, p_j):
    def integrand(s):
        term_x = X - (p_j.xa - math.sin(p_j.beta) * s)
        term_y = Y - (p_j.ya + math.cos(p_j.beta) * s)
        return (2 * term_x) / (term_x**2 + term_y**2)
    return integrate.quad_vec(integrand, 0.0, p_j.length)[0]

方案2：数值积分+广播运算（兼容旧版scipy） #

如果你的scipy版本低于1.7.0，可以通过对积分变量s采样，结合numpy广播计算所有点的被积函数值，再沿s轴做数值积分（比如梯形法、辛普森法）。

示例代码：

import numpy as np
from scipy import integrate

N = 3
x_start, x_end = -2.5, 2.5
y_start, y_end = -3.0, 3.0
x = np.linspace(x_start, x_end, N)
y = np.linspace(y_start, y_end, N)
X, Y = np.meshgrid(x, y)

def integral_x():
    # 对积分变量s采样，采样点数越多精度越高
    s_samples = np.linspace(0.0, 10, 1000)
    # 通过广播将s_samples扩展为三维数组，与X/Y匹配计算被积函数值
    integrand_vals = (2 * (X - s_samples[:, np.newaxis, np.newaxis])) / (
        (X - s_samples[:, np.newaxis, np.newaxis])**2 + 
        (Y - s_samples[:, np.newaxis, np.newaxis])**2
    )
    # 沿s轴（axis=0）执行辛普森数值积分
    return integrate.simpson(integrand_vals, x=s_samples, axis=0)

result = integral_x()
print(result.shape)  # 输出 (3, 3)

这种方法的精度取决于采样点数量，需要在精度和计算速度之间做平衡。

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