Python: Evaluating Integral for array

Question

import numpy

import matplotlib.pyplot as plt

from scipy import integrate

def f(x,y):
    return x*y + x**2 

def integral(x,y):
    I = integrate.quad(f, 0, x, args=(y,))[0]
    return I

def gau(x,y):
    return (1+x)*integral(x,y)


xlist = numpy.linspace(-3.0, 3.0, 100)
ylist = numpy.linspace(-3.0, 3.0, 100)
X, Y = numpy.meshgrid(xlist, ylist)
Z = gau(2, Y)

print(Z)

I keep on getting the error message "Supplied function does not return a valid float." , I think the problem is that I try to pass an array to the quad function. I thought about evaluating the integral for every entry of the array with something like that:

yi=numpy.linspace(-3.0,3.0,100)
for i, item in enumerate(yi):
    return integral[i]=integrate.quad(f,0,x,args=(yi,))[0]

It doesn't work but is it the right way? Any other/better suggestions?

Answer 1

You could use a universal function (see https://docs.scipy.org/doc/numpy/reference/ufuncs.html) which operates on arrays element-by-element. You can create these universal functions from any function using the frompyfunc function (https://docs.scipy.org/doc/numpy/reference/generated/numpy.frompyfunc.html):

ugau = numpy.frompyfunc(gau,2,1)
Z=ugau(X,Y)

Answer 2

It if your f() that does not provide a valid float when passed an array, not the scipy.integral itself;

why do you pass an array to your f() ?

Answer 3

You can use quadpy (one of my projects). quadpy is fully vectorized with respect to the dimensionality of the function range and the domains, so you can plug in a function that returns a vector and integrate that function over many intervals at once. You just have to make sure that the input function deals with vectorized input correctly. In your case, that would be

import numpy
import quadpy


def f(x, y):
    return numpy.multiply.outer(y, x) + numpy.multiply.outer(numpy.ones_like(y), x ** 2)


def integral(x, y):
    scheme = quadpy.line_segment.gauss_legendre(5)
    intervals = numpy.array([numpy.zeros_like(x), x])
    out = scheme.integrate(lambda t: f(t, y), intervals)
    return out


def gau(x, y):
    return (1 + x) * integral(x, y)


xlist = numpy.linspace(-3.0, 3.0, 100)
ylist = numpy.linspace(-3.0, 3.0, 100)
Z = gau(2, ylist)

print(Z)

You also insert xlist instead of 2 here to compute it all at once.