2D numerical integration in Python

Question

I know how to perform a double integral in python

import numpy as np
import scipy.integrate as integrate

integrate.dblquad(x*y, 0, 1, lambda x: -np.sqrt(1-x**2), lambda x: np.sqrt(1-x**2))

where x and y are, say, (200,) numpy arrays.

However, what if the integrand (x*y) above is a 2D array rather than a function? In my case, I have an array Z which has a value at every coordinate (x,y), i.e. it has shape (200,200). However, I do not know in advance the continuous function(s) that it would correspond to.

How would I perform this integral? Thanks.

Integration is over a continuous domain by definition; over a discrete set of values it's just a summation. So can't you just sum up all the values in the array multiplied by the area of each grid cell? — Thomas
– Thomas, Commented Jul 12, 2019 at 9:43
@Thomas thanks for your comment. You are absolutely right. I am stupid and also exhausted from work, which is a dangerous mix. Thanks again. — martin
– martin, Commented Jul 12, 2019 at 9:59
There is a similar question, for which nested 1D integration using Simpson's method is proposed (same idea as the direct sum, but accuracy is improved by assuming that the function is locally quadratic). However, the integration domain is still rectangular... — xdze2
– xdze2, Commented Jul 12, 2019 at 12:28

Nico Schlömer · Accepted Answer · 2019-07-21 13:29:32Z

2

So you're trying to integrate a function over the unit disk? You can use quadpy (a project of mine) to do that. It does support vector-valued integrands, too:

import numpy
import quadpy

scheme = quadpy.disk.lether(6)  # there are many other schemes available
# scheme.show()
val = scheme.integrate(
    lambda x: [numpy.exp(x[0] + x[1]), numpy.cos(x[0])],
    [0.0, 0.0],  # disk center
    1.0  # disk radius
)
print(val)

[3.99523707 2.76491937]

edited Jul 21, 2019 at 13:29

answered Jul 16, 2019 at 8:22

Nico Schlömer

59.7k35 gold badges216 silver badges291 bronze badges

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

4 Comments

martin Over a year ago

Following what was pointed out in the comments, I just did a double summation and it worked fine. However, quadpy seems interesting. Is there anywhere I can find an explanation on the syntax used? I am unsure what, for example, [0.0, 0.0], 1.0 or x[0] and x[1] stand for in the above code. Thanks

Nico Schlömer Over a year ago

center, radius of the disk. x[0] and x[1] are simply the coordinates (x and y if you want).

martin Over a year ago

What if the function I want to integrate is not actually a function, but an array of values corresponding to various points within the domain of integration? Thanks.

Nico Schlömer Over a year ago

If the points where the function is prescribed are fixed, you're outside of the domain of classical numerical integration. There not too much you can do there. One sensible thing would be to construct the Voronoi tesselation of your domain with the points inside and sum them up weighted with the Voronoi volumes.

Collectives™ on Stack Overflow

2D numerical integration in Python

1 Answer 1

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related