Numerical Integration

AUTHORS:

  • Josh Kantor (2007-02): first version
  • William Stein (2007-02): rewrite of docs, conventions, etc.
  • Robert Bradshaw (2008-08): fast float integration
  • Jeroen Demeyer (2011-11-23): trac ticket #12047: return 0 when the integration interval is a point; reformat documentation and add to the reference manual.
class sage.calculus.integration.PyFunctionWrapper

Bases: object

class sage.calculus.integration.compiled_integrand

Bases: object

sage.calculus.integration.monte_carlo_integral(func, xl, xu, calls, algorithm='plain', params=None)

Integrate func.

Integrate func over the dim-dimensional hypercubic region defined by the lower and upper limits in the arrays xl and xu, each of size dim. The integration uses a fixed number of function calls calls and obtains random sampling points using the default gsl’s random number generator.

ALGORITHM: Uses calls to the GSL (GNU Scientific Library) C library. Documentation can be found in [GSL] chapter “Monte Carlo Integration”.

INPUT:

  • func – The function to integrate
  • params – used to pass parameters to your function
  • xl – list of lower limits
  • xu – list of upper limits
  • calls – Number of functions calls used.
  • algorithm – valid choices are:
    • ‘plain’ – The plain Monte Carlo algorithm samples points randomly
      from the integration region to estimate the integral and its error.
    • ‘miser’ – The MISER algorithm of Press and Farrar is based on
      recursive stratified sampling
    • ‘vegas’ – The VEGAS algorithm of Lepage is based on importance
      sampling.

EXAMPLES:

sage: x, y = SR.var('x,y')
sage: monte_carlo_integral(x*y, [0,0], [2,2], 10000)   # abs tol 0.1
(4.0, 0.0)
sage: integral(integral(x*y, (x,0,2)), (y,0,2))
4

An example with a parameter:

sage: x, y, z = SR.var('x,y,z')
sage: monte_carlo_integral(x*y*z, [0,0], [2,2], 10000, params=[1.2])   # abs tol 0.1
(4.8, 0.0)

Integral of a constant:

sage: monte_carlo_integral(3, [0,0], [2,2], 10000)   # abs tol 0.1
(12, 0.0)

Test different algorithms:

sage: x, y, z = SR.var('x,y,z')
sage: f(x,y,z) = exp(z) * cos(x + sin(y))
sage: for algo in ['plain', 'miser', 'vegas']:  # abs tol 0.01
....:   monte_carlo_integral(f, [0,0,-1], [2,2,1], 10^6, algorithm=algo)
(-1.06, 0.01)
(-1.06, 0.01)
(-1.06, 0.01)

Tests with Python functions:

sage: def f(u, v): return u * v
sage: monte_carlo_integral(f, [0,0], [2,2], 10000)  # abs tol 0.1
(4.0, 0.0)
sage: monte_carlo_integral(lambda u,v: u*v, [0,0], [2,2], 10000)  # abs tol 0.1
(4.0, 0.0)
sage: def f(x1,x2,x3,x4): return x1*x2*x3*x4
sage: monte_carlo_integral(f, [0,0], [2,2], 1000, params=[0.6,2])  # abs tol 0.2
(4.8, 0.0)

AUTHORS:

  • Vincent Delecroix
  • Vincent Klein
sage.calculus.integration.numerical_integral(func, a, b=None, algorithm='qag', max_points=87, params=[], eps_abs=1e-06, eps_rel=1e-06, rule=6)

Returns the numerical integral of the function on the interval from a to b and an error bound.

INPUT:

  • a, b – The interval of integration, specified as two numbers or as a tuple/list with the first element the lower bound and the second element the upper bound. Use +Infinity and -Infinity for plus or minus infinity.

  • algorithm – valid choices are:

    • ‘qag’ – for an adaptive integration
    • ‘qags’ – for an adaptive integration with (integrable) singularities
    • ‘qng’ – for a non-adaptive Gauss-Kronrod (samples at a maximum of 87pts)

  • max_points – sets the maximum number of sample points

  • params – used to pass parameters to your function

  • eps_abs, eps_rel – sets the absolute and relative error tolerances which satisfies the relation |RESULT - I|  <= max(eps_abs, eps_rel * |I|), where I = \int_a^b f(x) d x.

  • rule – This controls the Gauss-Kronrod rule used in the adaptive integration:

    • rule=1 – 15 point rule
    • rule=2 – 21 point rule
    • rule=3 – 31 point rule
    • rule=4 – 41 point rule
    • rule=5 – 51 point rule
    • rule=6 – 61 point rule

    Higher key values are more accurate for smooth functions but lower key values deal better with discontinuities.

OUTPUT:

A tuple whose first component is the answer and whose second component is an error estimate.

REMARK:

There is also a method nintegral on symbolic expressions that implements numerical integration using Maxima. It is potentially very useful for symbolic expressions.

EXAMPLES:

To integrate the function \(x^2\) from 0 to 1, we do

sage: numerical_integral(x^2, 0, 1, max_points=100)
(0.3333333333333333, 3.700743415417188e-15)

To integrate the function \(\sin(x)^3 + \sin(x)\) we do

sage: numerical_integral(sin(x)^3 + sin(x),  0, pi)
(3.333333333333333, 3.700743415417188e-14)

The input can be any callable:

sage: numerical_integral(lambda x: sin(x)^3 + sin(x),  0, pi)
(3.333333333333333, 3.700743415417188e-14)

We check this with a symbolic integration:

sage: (sin(x)^3+sin(x)).integral(x,0,pi)
10/3

If we want to change the error tolerances and gauss rule used:

sage: f = x^2
sage: numerical_integral(f, 0, 1, max_points=200, eps_abs=1e-7, eps_rel=1e-7, rule=4)
(0.3333333333333333, 3.700743415417188e-15)

For a Python function with parameters:

sage: f(x,a) = 1/(a+x^2)
sage: [numerical_integral(f, 1, 2, max_points=100, params=[n]) for n in range(10)]  # random output (architecture and os dependent)
[(0.49999999999998657, 5.5511151231256336e-15),
 (0.32175055439664557, 3.5721487367706477e-15),
 (0.24030098317249229, 2.6678768435816325e-15),
 (0.19253082576711697, 2.1375215571674764e-15),
 (0.16087527719832367, 1.7860743683853337e-15),
 (0.13827545676349412, 1.5351659583939151e-15),
 (0.12129975935702741, 1.3466978571966261e-15),
 (0.10806674191683065, 1.1997818507228991e-15),
 (0.09745444625548845, 1.0819617008493815e-15),
 (0.088750683050217577, 9.8533051773561173e-16)]
sage: y = var('y')
sage: numerical_integral(x*y, 0, 1)
Traceback (most recent call last):
...
ValueError: The function to be integrated depends on 2 variables (x, y),
and so cannot be integrated in one dimension. Please fix additional
variables with the 'params' argument

Note the parameters are always a tuple even if they have one component.

It is possible to integrate on infinite intervals as well by using +Infinity or -Infinity in the interval argument. For example:

sage: f = exp(-x)
sage: numerical_integral(f, 0, +Infinity)  # random output
(0.99999999999957279, 1.8429811298996553e-07)

Note the coercion to the real field RR, which prevents underflow:

sage: f = exp(-x**2)
sage: numerical_integral(f, -Infinity, +Infinity)  # random output
(1.7724538509060035, 3.4295192165889879e-08)

One can integrate any real-valued callable function:

sage: numerical_integral(lambda x: abs(zeta(x)), [1.1,1.5])  # random output
(1.8488570602160455, 2.052643677492633e-14)

We can also numerically integrate symbolic expressions using either this function (which uses GSL) or the native integration (which uses Maxima):

sage: exp(-1/x).nintegral(x, 1, 2)  # via maxima
(0.50479221787318..., 5.60431942934407...e-15, 21, 0)
sage: numerical_integral(exp(-1/x), 1, 2)
(0.50479221787318..., 5.60431942934407...e-15)

We can also integrate constant expressions:

sage: numerical_integral(2, 1, 7)
(12.0, 0.0)

If the interval of integration is a point, then the result is always zero (this makes sense within the Lebesgue theory of integration), see trac ticket #12047:

sage: numerical_integral(log, 0, 0)
(0.0, 0.0)
sage: numerical_integral(lambda x: sqrt(x), (-2.0, -2.0) )
(0.0, 0.0)

In the presence of integrable singularity, the default adaptive method might fail and it is advised to use 'qags':

sage: b = 1.81759643554688
sage: F(x) = sqrt((-x + b)/((x - 1.0)*x))
sage: numerical_integral(F, 1, b)
(inf, nan)
sage: numerical_integral(F, 1, b, algorithm='qags')    # abs tol 1e-10
(1.1817104238446596, 3.387268288079781e-07)

AUTHORS:

  • Josh Kantor
  • William Stein
  • Robert Bradshaw
  • Jeroen Demeyer

ALGORITHM: Uses calls to the GSL (GNU Scientific Library) C library. Documentation can be found in [GSL] chapter “Numerical Integration”.