在 Python 中集成返回数组的函数答案

【问题标题】：Integrating functions that return an array in Python在 Python 中集成返回数组的函数
【发布时间】：2019-11-05 17:19:01
【问题描述】：

我有很多数据要整合，我想找到一种只用矩阵来完成这一切的方法，并且愿意在准确性上做出妥协以提高性能。我的想法是这样的：

import numpy
import scipy

a = np.array([1,2,3])

def func(x):
    return x**2 + x

def func2(x):
    global a
    return a*x

def integrand(x):
    return func(x)*func2(x)

integrated = quad(integrand, 0, 1)

所以我试图整合来自integrand的数组中的每个元素。

我知道有可能像这样使用numpy.vectorize()：

integrated = numpy.vectorize(scipy.integrate.quad)(integrand, 0, 1)

但我无法让它工作。有没有办法在python中做到这一点？

解决方案

现在我学到了更多的python，如果有人碰巧稳定它并且有同样的问题，我可以回答这个问题。这样做的方法是编写函数，就好像它们将采用标量值而不是向量作为输入一样。所以按照我上面的代码，我们会得到类似

import numpy as np
import scipy.integrate.quad

a = np.array([1, 2, 3]) # arbitrary array, can be any size

def func(x):
    return x**2 + x

def func2(x, a):
    return a*x

def integrand(x, a):
    return func(x)*func2(x, a)

def integrated(a):
    integrated, tmp = scipy.integrate.quad(integrand, 0, 1, args = (a))
    return integrated

def vectorizeInt():
    global a
    integrateArray = []
    for i in range(len(a)):
        integrate = integrated(a[i])
        integrateArray.append(integrate)
    return integrateArray

并不是您要积分的变量必须是函数的第一个输入。这是 scipy.integrate.quad 所必需的。如果您正在集成一个方法，它是典型的self 之后的第二个参数（即x 集成在def integrand(self, x, a): 中）。此外，args = (a) 是告诉quad 函数integrand 中a 的值所必需的。如果integrand 有很多参数，比如def integrand(x, a, b, c, d):，您只需将参数按顺序排列在args 中。那就是args = (a, b, c, d)。

【问题讨论】：

标签： python numpy scipy

【解决方案1】：

vectorize 不会帮助提高使用quad 的代码的性能。要使用quad，您必须为integrate 返回的值的每个组件单独调用它。

对于矢量化但不太准确的近似值，您可以使用numpy.trapz 或scipy.integrate.simps。

您的函数定义（至少问题中显示的那个）是使用所有支持广播的 numpy 函数实现的，因此给定 [0, 1] 上的 x 值网格，您可以这样做：

In [270]: x = np.linspace(0.0, 1.0, 9).reshape(-1,1)

In [271]: x
Out[271]: 
array([[ 0.   ],
       [ 0.125],
       [ 0.25 ],
       [ 0.375],
       [ 0.5  ],
       [ 0.625],
       [ 0.75 ],
       [ 0.875],
       [ 1.   ]])

In [272]: integrand(x)
Out[272]: 
array([[ 0.        ,  0.        ,  0.        ],
       [ 0.01757812,  0.03515625,  0.05273438],
       [ 0.078125  ,  0.15625   ,  0.234375  ],
       [ 0.19335938,  0.38671875,  0.58007812],
       [ 0.375     ,  0.75      ,  1.125     ],
       [ 0.63476562,  1.26953125,  1.90429688],
       [ 0.984375  ,  1.96875   ,  2.953125  ],
       [ 1.43554688,  2.87109375,  4.30664062],
       [ 2.        ,  4.        ,  6.        ]])

也就是说，通过将x 设为形状为(n, 1) 的数组，integrand(x) 返回的值的形状为(n, 3)。 a 中的每个值对应一列。

您可以使用axis=0 将该值传递给numpy.trapz() 或scipy.integrate.simps()，以获得积分的三个近似值。您可能需要更精细的网格：

In [292]: x = np.linspace(0.0, 1.0, 101).reshape(-1,1)

In [293]: np.trapz(integrand(x), x, axis=0)
Out[293]: array([ 0.583375,  1.16675 ,  1.750125])

In [294]: simps(integrand(x), x, axis=0)
Out[294]: array([ 0.58333333,  1.16666667,  1.75      ])

与重复调用quad相比：

In [296]: np.array([quad(lambda t: integrand(t)[k], 0, 1)[0] for k in range(len(a))])
Out[296]: array([ 0.58333333,  1.16666667,  1.75      ])

你的函数integrate（我假设这只是一个例子）是一个三次多项式，Simpson's rule 给出了精确的结果。一般来说，不要指望simps 给出如此准确的答案。

【讨论】：

为了完整起见，我会提到scipy.integrate.romb，它可以从相同的数据中挤出更准确的结果（假设它有2**k+1 样本）。
我知道这个答案已经有将近五年的历史了，但它今天确实帮助了我。谢谢！ :)

【解决方案2】：

quadpy（我的一个项目）是完全矢量化的。安装

pip install quadpy

然后做

import numpy
import quadpy


def integrand(x):
    return [numpy.sin(x), numpy.exp(x)]  # ,...


res, err = quadpy.quad(integrand, 0, 1)
print(res)
print(err)

[0.45969769 1.71828183]
[1.30995437e-20 1.14828375e-19]

【讨论】：