python - 在 Python 中集成返回数组的函数

Question

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

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)。

Answer 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 中的每个值都有一列。

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

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 会给出如此准确的答案。