scipy.integrate.quad on functions returning numpy arrays

Question

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.

possible duplicate of using SciPy to integrate a function that returns a matrix or array — Saullo G. P. Castro
– Saullo G. P. Castro, Commented Mar 22, 2014 at 19:59

sergpolly · Accepted Answer · 2014-03-14 16:37:24Z

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

edited Mar 14, 2014 at 16:37

answered Mar 14, 2014 at 16:14

sergpolly

1128 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

felippe Over a year ago

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.

sergpolly Over a year ago

Are you integrating over x,y or a,b?

sergpolly Over a year ago

[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?

felippe Over a year ago

The integration is over a and b. x and y are arrays.

sergpolly Over a year ago

f = lambda a,b,x,y: np.exp(ax+by); [[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)]

Collectives™ on Stack Overflow

scipy.integrate.quad on functions returning numpy arrays

1 Answer 1

5 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related