如何用NumPy/SciPy实现标量输入矩阵输出函数的积分？

我有一个输入为标量、输出为矩阵的函数（以下是2x2的最小可复现示例，实际函数的矩阵维度为问题特定，如100x100），希望对该函数在某一区间上进行积分。

尝试直接使用scipy.integrate.quad：

import numpy as np
from scipy.integrate import quad

def func(x):
    return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]])

quad(func, 0, np.pi)

运行后触发错误：

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/scipy/integrate/_quadpack_py.py", line 465, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/scipy/integrate/_quadpack_py.py", line 577, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: only size-1 arrays can be converted to Python scalars

尝试向量化quad操作仍报错：

grid_quad = np.vectorize(quad)
grid_quad(func, 0, np.pi)

仅当将被积函数拆分为函数数组时才可行：

efunc = [[np.cos, np.sin],[np.sin, np.cos]]
grid_quad(efunc, 0, np.pi)

得到预期输出：

(array([[4.92255263e-17, 2.00000000e+00],
       [2.00000000e+00, 4.92255263e-17]]), array([[2.21022394e-14, 2.22044605e-14],
       [2.22044605e-14, 2.21022394e-14]]))

需求：对输出矩阵的每个元素对应的函数逐个积分，输出积分值矩阵，求符合Pythonic风格的实现方式。

方法1：使用quad_vec（Scipy 1.10+ 推荐） #

Scipy 1.10及以上版本提供了quad_vec，专门支持向量/矩阵值函数的积分，无需额外处理，直接适配你的需求：

import numpy as np
from scipy.integrate import quad_vec

def func(x):
    return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]])

# 直接积分矩阵值函数
result_matrix, error_matrix = quad_vec(func, 0, np.pi)

print("积分结果矩阵：")
print(result_matrix)
print("\n误差矩阵：")
print(error_matrix)

输出与预期完全一致，这是最简洁的Pythonic实现。

方法2：向量化单元素积分逻辑 #

如果你的Scipy版本较低，可通过np.vectorize封装针对矩阵单个元素的积分逻辑，无需拆分原函数：

import numpy as np
from scipy.integrate import quad

def func(x):
    return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]])

# 获取函数输出的矩阵形状
sample_output = func(0)
rows, cols = sample_output.shape

# 定义针对矩阵(i,j)位置元素的积分函数
def integrate_single_element(i, j):
    return quad(lambda x: func(x)[i, j], 0, np.pi)

# 向量化该函数
vec_integrate = np.vectorize(integrate_single_element)

# 生成索引网格并计算
i_indices, j_indices = np.meshgrid(np.arange(rows), np.arange(cols), indexing='ij')
results, errors = vec_integrate(i_indices, j_indices)

print(results)
print(errors)

方法3：嵌套列表推导式（直观易读） #

对于小维度矩阵，嵌套列表推导式可读性强，也是Pythonic的实现方式：

import numpy as np
from scipy.integrate import quad

def func(x):
    return np.array([[np.cos(x), np.sin(x)], [np.sin(x), np.cos(x)]])

sample_output = func(0)
rows, cols = sample_output.shape

# 用列表推导式生成积分结果和误差
results, errors = zip(*[[quad(lambda x: func(x)[i,j], 0, np.pi) for j in range(cols)] for i in range(rows)])
result_matrix = np.array(results)
error_matrix = np.array(errors)

print(result_matrix)
print(error_matrix)

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