3

I have a function returning a numpy array that I need to integrate. I would like to use the scipy.integrate.dblquad, but it requires the function to return a float.

I tried to use the numpy.vectorize on dblquad, it doesn't work. After some thought, I realized the problem with this approach is that the vectorized dblquad requires a vector of functions, not a function returning vectors.

I want to solve this by any means, except splitting the arrays and integrating case by case (already did that and the code is a mess). The first thing that came in my mind was to transform the function returning the array in a array returning equivalent functions, but I don't know if this is possible.

I'll put some illustrative code. Below, x and y are numpy arrays, a and b are floats.

import numpy as np
from scipy.integrate import dblquad

def foo(a, b, x, y):
    return np.exp(a*x + b*y)

I would like to do 

x = np.random.rand((3,3))
y = np.random.rand((3,3))
a = 1.5
b = 3.0
I = dblquad(foo, 0, 1, lambda x: 0, lambda x: 1, args=(x,y))

to get an array I with shape (3,3) with entry I[i][j] being the integral of
foo(a,b,x[i][j],y[i][j]), over 0 and 1 for the first two arguments.

I believe there is a smart way to do this with np.vectorize, but I couldn't think of any.

felippe
  • 353
  • 3
  • 7
  • possible duplicate of [using SciPy to integrate a function that returns a matrix or array](http://stackoverflow.com/questions/17070333/using-scipy-to-integrate-a-function-that-returns-a-matrix-or-array) – Saullo G. P. Castro Mar 22 '14 at 19:59

1 Answers1

2

This what worked for me:

f = lambda x,y,a,b: np.exp(a*x+b*y)
dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(1.0,1.0), epsabs=1.49e-08, epsrel=1.49e-08)

f - 2D function with x,y being variables a,b - just parameters. So we integrate x from zero to one - second and third arg in dblquad function, then y from zero to one - 4th and 5th args in dblquad, but in a form of functions; then you pass your a,b in a form of args=(1.0,1.0). The rest is the default stuff

If you would like to vary you parameters over some range of values you can do something like this:

[dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(a,b)) for a in range(2) for b in range(2)]

this will get you a flat array with the results (integrals) corresponding to (a,b):

 [(0, 0), (0, 1), (1, 0), (1, 1)]

You also can do:

  [[dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(a,b)) for a in range(2)] for b in range(2)]

Which will yield a 2D structured list of results corresponding (a,b):

  [[(0, 0), (0, 1)], [(1, 0), (1, 1)]]
sergpolly
  • 112
  • 8
  • I meant x and y as arrays, so foo value is also an array. Try your suggestion with x = np.ones((2,2)) and y = np.ones((2,2)). That way quad doesn't work and that is what I want. – felippe Mar 14 '14 at 16:23
  • Are you integrating over x,y or a,b? – sergpolly Mar 14 '14 at 16:26
  • [dblquad(f, 0, 1, lambda x:0, lambda x:1, args=(a,b)) for a in range(2) for b in range(2)] - list comrehension for varying your parameters... isn't that what you want? – sergpolly Mar 14 '14 at 16:29
  • The integration is over a and b. x and y are arrays. – felippe Mar 14 '14 at 16:59
  • f = lambda a,b,x,y: np.exp(a*x+b*y); [[dblquad(f, 0, 1, lambda var:0, lambda var:1, args=(x,y)) for x in np.arange(3)] for y in np.arange(3)] – sergpolly Mar 14 '14 at 18:40