pythonscipynumerical-integrationquad

Python: Evaluating Integral for array


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?


Solution

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