Python/Scipy - 与 Quad Along Axis 集成

Question

我有一个二维数组。 “xy”平面是从 (-1,-1) 到 (1,1) 的网格。我想在函数取决于点坐标的每个点处计算和积分。

我知道对于离散数据，我可以使用 simps 或 trapz 并指定一个轴进行积分(参见示例)。 scipy.integrate.quad 是否可以在不使用如下所示的丑陋循环的情况下实现？

import numpy as np
import scipy as sp
import scipy.integrate

x = np.linspace(-1, 1, 10)
y = np.linspace(-1, 1, 10)
X,Y = np.meshgrid(x, y)

z = np.linspace(1, 10, 100)

# Integrate discrete values using simps
def func(z):
   return (X - z) / ((X - z)**2. + Y**2)
res1 = sp.integrate.simps(func(z.reshape(-1, 1, 1)), z, axis=0)
print(res1)

# Integrate the function using quad at each point in the xy plane
res2 = np.zeros(X.shape)
for i in range(res2.shape[0]):
   for j in range(res2.shape[1]):
      def func2(z):
         return (X[i,j] - z) / ((X[i,j] - z)**2. + Y[i,j]**2)
      res2[i,j] = sp.integrate.quad(func2, 1, 10)[0]
print(res2)

Answer 1

使用 Cubature method来自 Steven Johnson 教授，我 wrapped using Cython ，您可以立即实现集成:

import numpy as np

from cubature import cubature

x = np.linspace(-1, 1, 10)
y = np.linspace(-1, 1, 10)
X,Y = np.meshgrid(x, y)

z = np.linspace(1, 10, 100)

def func(z):
   return (X.ravel() - z) / ((X.ravel() - z)**2. + Y.ravel()**2)

res = cubature(1, func, np.array([1.]), np.array([10.]))[0].reshape(X.shape)