python - Python/Scipy - 与沿轴的四边形集成

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

使用我使用 Cython 包装的 Steven Johnson 教授的Cubature 方法，您可以立即实现集成：

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)

Answer 2

使用线性间隔的点不是很有效。一个更好的主意是转向四边形的正交方案，例如，here。例如：

import numpy
import quadpy

def f(x):
   return (x[0] - 1) / ((x[0] - 1)**2 + x[1]**2)

quad = numpy.array([
    [-1, -1],
    [+1, -1],
    [+1, +1],
    [-1, +1],
    ])

scheme = quadpy.quadrilateral.Stroud(5)

val = quadpy.quadrilateral.integrate(f, quad, scheme)