7

I'm looking for a super duper numerical quadrature function. It should have the following three properties:

  • Adaptive - it automatically adjusts the density of sampling points to fit the integrand. This is absolutely necessary because my integrand is very nonuniform and expensive to compute.
  • Vectorized - it calls the integrand on lists of sample points rather than one point at a time, for efficiency.
  • Able to handle vector-valued functions - all components of the vector-valued integrand are computed at the same time for no additional cost, so it makes no sense to integrate all the components separately.

In addition, it should be:

  • 2D - the integral I want to compute is a double integral over a planar region, and I want to be able to specify an overall (relative) tolerance for the whole integral and have it manage the error budget appropriately.

Does anybody know of a library that has such a function? Even two or three of the four properties would be better than nothing.

I'm using Python and SciPy, so if it already works with Python that's a bonus. (But I'm also able to write glue code to let it call my integrand if necessary.)

Andrew Mao
  • 35,740
  • 23
  • 143
  • 224
Keenan Pepper
  • 321
  • 1
  • 3
  • 10
  • 1
    Alas, no answers yet! I am writing my own numerical integration algorithm in C#. It is adaptive and handles N-dimensions, but is not vectorized, however. Getting the terminating conditions (tolerance) right is proving difficult. – Paul Chernoch Oct 01 '12 at 21:38
  • @Keenan Pepper Maybe the process described [in this question can give you some insight](http://stackoverflow.com/questions/14441541/performance-behaviour-of-vectorized-functions-in-numpy/16551313#16551313) – Saullo G. P. Castro May 18 '13 at 12:36

3 Answers3

5

I just implemented vectorized adaptive quadrature for 1D and 2D domains in quadpy. All you need to provide is a triangulation of your domain and the function you want to integrate. It may be vector-valued.

Install quadpy with

pip3 install quadpy

and run

import numpy
import quadpy


triangles = numpy.array(
    [[[0.0, 0.0], [1.0, 0.0]], [[1.0, 0.0], [1.0, 1.0]], [[0.0, 1.0], [0.0, 1.0]]]
)

val, error_estimate = quadpy.triangle.integrate_adaptive(
    lambda x: [numpy.sin(x[0]), numpy.exp(x[0])], triangles, 1.0e-10
)

print(val)
print(error_estimate)

This gives

[ 0.45969769  1.71828183]
[  7.10494337e-12   3.68776277e-11]
Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
  • 1
    I just did a copy and paste of this code snippet and it did not work with version 0.13.2. Turns out that `adaptive_integrate` should be `integrate_adaptive`. – Pedro H. N. Vieira Aug 16 '19 at 11:45
1

I used this library, it does everything you want, except it is written in C. But it also has an R interface, so maybe you can call R from Python (that is possible).

http://ab-initio.mit.edu/wiki/index.php/Cubature_(Multi-dimensional_integration)

Or, you can call the library using ctypes (it is not straight forward, but it is doable).

mher
  • 369
  • 3
  • 7
1

The quadrature function in scipy.integrate satisfies the first two requirements of what you are looking for. The similar romberg function uses a different method.

Other functions only satisfy one of the requirements:

  • The similarly-named quad function does adaptive quadrature, but only supports a function with a scalar argument. You can pass it a ctypes function for increased performance, but normal Python functions will be very slow.
  • The simps function and related sampling methods can be passed a vector of (typically evenly-spaced) samples, but aren't adaptive.

The third requirement you listed (simultaneous integral of a vector-valued function) is a bit esoteric and conflicts with the ability to accept a vectorized function in the first place (the function argument would have to take a matrix!) Similarly, the ability to compute a double integral would complicate the specification of the function significantly.

In most cases, the quadrature function would be the way to go.

Andrew Mao
  • 35,740
  • 23
  • 143
  • 224