集成在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的值.如果有很多参数,请说你只需按顺序放入参数.那就是.aintegrandintegranddef integrand(x, a, b, c, d):argsargs = (a, b, c, d)

Answer 1

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

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

您的函数定义（至少是问题中显示的函数定义）是​​使用都支持广播的 numpy 函数实现的，因此给定x[0, 1] 上的值网格，您可以执行以下操作：

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)具有 shape (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（我认为这只是一个例子）是三次多项式，辛普森规则给出了准确的结果。一般来说，不要期望simps给出如此准确的答案。