在 Python 中集成返回数组的函数 [英] Integrating functions that return an array in Python

问题描述

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

I have a lot of data to integrate over and would like to find a way of doing it all with just matrices, and would be willing to compromise on accuracy for a performance boost. What I have in mind is something like this:

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 的数组中的每个元素.

So I am trying to integrate each element in the array that comes out of integrand.

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

I'm aware that there is a possibility of using numpy.vectorize() like this:

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

但我不能让它工作.有没有办法在python中做到这一点?

but I can't get that working. Is there a way to do this in python?

解决方案

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

Well now that I learnt a bit more python I can answer this question if anyone happens to stable upon it and has the same question. The way to do it is to write the functions as though they are going to take scalar values, and not vectors as inputs. So follow from my code above, what we would have is something like

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

Not that the variable which you are integrating over must be the first input to the function. This is required for scipy.integrate.quad. If you are integrating over a method, it is the second argument after the typical self (i.e. x is integrated in def integrand(self, x, a):). Also the args = (a) is necessary to tell quad the value of a in the function integrand. If integrand has many arguments, say def integrand(x, a, b, c, d): you simply put the arguments in order in args. So that would be args = (a, b, c, d).

推荐答案

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

vectorize won't provide help with improving the performance of code that uses quad. To use quad, you'll have to call it separately for each component of the value returned by integrate.

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

For a vectorized but less accurate approximation, you can use numpy.trapz or scipy.integrate.simps.

您的函数定义(至少是问题中显示的函数定义)是​​使用都支持广播的 numpy 函数实现的，因此给定 [0, 1] 上 x 值的网格，您可以执行这个:

Your function definition (at least the one shown in the question) is implemented using numpy functions that all support broadcasting, so given a grid of x values on [0, 1], you can do this:

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 中的每个值对应一列.

That is, by making x an array with shape (n, 1), the value returned by integrand(x) has shape (n, 3). There is one column for each value in a.

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

You can pass that value to numpy.trapz() or scipy.integrate.simps(), using axis=0, to get the three approximations of the integrals. You'll probably want a finer grid:

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 的重复调用进行比较:

Compare that to repeated calls to 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(我认为这只是一个例子)是一个三次多项式，对于它 辛普森法则给出了准确的结果.一般来说，不要指望 simps 给出如此准确的答案.

Your function integrate (which I assume is just an example) is a cubic polynomial, for which Simpson's rule gives the exact result. In general, don't expect simps to give such an accurate answer.

这篇关于在 Python 中集成返回数组的函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！