4

I am trying to compute a double integral (over a triangle with nodes at (0,0), (0,1), (1,0)) using Gaussian quadrature of order n. However, running

import scipy.integrate as integrate
f = lambda x,y: x+y
inside = lambda x: integrate.fixed_quad(f, 0, 1-x, args=(x,), n=5)[0]
outside = integrate.fixed_quad(inside, 0, 1, n=5)

gives

Traceback (most recent call last): File "", line 1, in File "/Users/username/anaconda/lib/python3.5/site-packages/scipy/integrate/quadrature.py", line 82, in fixed_quad return (b-a)/2.0 * np.sum(w*func(y, *args), axis=0), None File "", line 1, in File "/Users/username/anaconda/lib/python3.5/site-packages/scipy/integrate/quadrature.py", line 78, in fixed_quad if np.isinf(a) or np.isinf(b):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

This is the second part of the question Can scipy.integrate.fixed_quad compute integral with functional boundaries?.

7
  • 2
    are you looking for double integrals?
    – MB-F
    Feb 16, 2017 at 11:48
  • @kazemakase I'm trying to implement double integrals with parameter n, which is called Gaussian quadrature integration order. Double integrals can be computed using scipy.integrate.dblquad() but this is not what I'm looking for.
    – Jan
    Feb 16, 2017 at 11:52
  • Could you write down the integral you try to solve?
    – Cleb
    Feb 16, 2017 at 11:54
  • @Cleb int_0^1 int_0^(x-1) f(x,y) dy dx (integration of function f(x,y) over the triangle with nodes (0,0), (0,1), (1,0))
    – Jan
    Feb 16, 2017 at 12:06
  • using dblquad: result = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1-x)[0]
    – Jan
    Feb 16, 2017 at 12:13

2 Answers 2

2

The answer to your question is, yes, under certain conditions.

For demonstration purposes, I first choose different bounds than you (11 instead of 1 - x).

Generally, one can solve these types of integrals using dblquad:

area_dblquad = integrate.dblquad(lambda x, y: x + y, 0, 1, lambda x: 0, lambda x: 11)[0]

which here returns 66. That is not an option as you mentioned in the comments.

One can now do this integration stepwise and it works fine for quad as well as fixed_quad:

def integrand(x, y):
    return x + y

def fint_quad(x):
    return integrate.quad(integrand, 0, 11, args=(x, ))[0]

def fint_fixed_quad(x):
    return integrate.fixed_quad(integrand, 0, 11, args=(x, ), n=5)[0]

res_quad = integrate.quad(lambda x: fint_quad(x), 0, 1)
res_fixed_quad = integrate.fixed_quad(lambda x: fint_fixed_quad(x), 0, 1, n=5)

They both return 66 as well, as expected. That shows that it can work to compute double integrals using scipy.integrate.fixed_quad.

However, when one now changes the upper bound back to the one you had (from 11 to 1 - x), it still works for quad but crashes for fixed_quad:

area_dblquad = integrate.dblquad(lambda x, y: x + y, 0, 1, lambda x: 0, lambda x: 1 - x)[0]
res_quad = integrate.quad(lambda x: fint_quad(x), 0, 1)

both return 0.333333..., the call with fixed_quad results in the error you received. One can understand the error by looking on the source code:

x, w = _cached_roots_legendre(n)
x = np.real(x)
if np.isinf(a) or np.isinf(b):
    raise ValueError("Gaussian quadrature is only available for "
                         "finite limits.")
y = (b-a)*(x+1)/2.0 + a
return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None

When one prints a and b one gets:

a:  0
b:  1
a:  0
b:  [ 0.95308992  0.76923466  0.5         0.23076534  0.04691008]

So for the call with 1-x, b is actually a numpy array with n elements and one cannot compare an array to infinity, that's why it crashes. Whether that is an intended behavior or a bug, I can't answer; might be worth opening an issue on github.

4
  • That's a NumPy array, not a list. Feb 17, 2017 at 1:11
  • As for whether it's a bug, no. fixed_quad documents that the function you integrate "must accept vector inputs", which inside does not. Feb 17, 2017 at 1:13
  • @user2357112: Could very well be; but would it then not be printed as array(...)? But no matter whether it is a numpy array or a list, it still leads to the same issue when compared to infinity. Ah, very well spotted, thanks! I overlooked this part in the documentation.
    – Cleb
    Feb 17, 2017 at 1:16
  • It wouldn't be printed as array(...) with just print. You'd have to do print(repr(...)) to get the repr instead of str output. Feb 17, 2017 at 1:56
0

fixed_quad requires f to accept vector inputs. And the result should be the mapped values for the inputs (i.e. something like map(f, xs)). With this in mind, just ensure your inside function returns mapped values, and you're ready to go.

import scipy.integrate as integrate
f = lambda y,x: x+y
inside = lambda xs, n: np.array([integrate.fixed_quad(f, 0, 1-x, args=(x,), n=n)[0]
                                 for x in np.array(xs)])
order = 5
outside = integrate.fixed_quad(inside, 0, 1, n=order, args=(order,))

Also, be careful with the order of arguments for your integrand. Judging from your code arg=(x,), it seems that you want the inner integral to be done along y dimension. The first argument of the integrand is the dimension along which it is integrated. So it should be lambda y,x instead (note that this is also the shape of integrand expected by dblquad).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.