NAME

Math::Integral::Romberg - Scalar numerical integration

SYNOPSIS

    use Math::Integral::Romberg 'integral';
    sub f { my ($x) = @_; 1/($x ** 2 + 1) } # ... or whatever
    $area = integral(\&f, $x1, $x2);    # Short form
    $area = integral                    # Long form
     (\&f, $x1, $x2, $rel_err, $abs_err, $max_split, $min_split);

    # an alternative way of doing the long form
    $Math::Integral::Romberg::rel_err = $rel_err;
    $Math::Integral::Romberg::abs_err = $abs_err;
    $Math::Integral::Romberg::max_split = $max_split;
    $Math::Integral::Romberg::min_split = $min_split;
    $area = integral(\&f, $x1, $x2);

DESCRIPTION

integral() numerically estimates the integral of f() using Romberg integration, a faster relative of Simpson's method.

Parameters

$f

A reference to the function to be integrated.

$x1, $x2

The limits of the integration domain. &$f(x1) and &$f(x2) must be finite.

$rel_err

Maximum acceptable relative error. Estimates of relative and absolute error are based on a comparison of the estimate computed using 2**n + 1 points with the estimate computed using 2**(n-1) + 1 points.

Once $min_split has been reached (see below), computation stops as soon as relative error drops below $rel_err, absolute error drops below $abs_err, or $max_split is reached.

If not supplied, uses the value $Math::Integral::Romberg::rel_err whose default is 10**-10. The accuracy limit of double-precision floating point is about 10**-15.

$abs_err

Maximum acceptable absolute error. If not supplied, uses $Math::Integral::Romberg::abs_err, which defaults to 10**-20.

$max_split

At most 2 ** $max_split + 1 different sample x values are used to estimate the integral of f(). If not supplied, uses the value $Math::Integral::Romberg::max_split, which defaults to 16, corresponding to 65537 sample points.

$min_split

At least 2 ** $min_split + 1 different sample x values are used to estimate the integral of f(). If not supplied, uses the value of $Math::Integral::Romberg::max_split, which defaults to 5, corresponding to 33 sample points.

$Math::Integral::Romberg::return_point_count

This value defaults to 0. If you set it to 1, then when invoked in a list context, integral() will return a two-element list, containing the estimate followed by the number of sample points used to compute the estimate.

$Math::Integral::Romberg::abort

This value is set to 1 if neither the $rel_err nor the $abs_err thresholds are reached before computation stops. Once set, this variable remains set until you reset it to 0.

Default values

Using the long form of integral() sets the convergence parameters for that call only - you must use the package-qualified variable names (e.g. $Math::Integral::Romberg::abs_tol) to change the values for all calls.

About the Algorithm

Romberg integration uses progressively higher-degree polynomial approximations each time you double the number of sample points. For example, it uses a 2nd-degree polynomial approximation (as Simpson's method does) after one split (2**1 + 1 sample points), and it uses a 10th-degree polynomial approximation after five splits (2**5 + 1 sample points). Typically, this will greatly improve accuracy (compared to simpler methods) for smooth functions, while not making much difference for badly behaved ones.

AUTHOR

Eric Boesch (ericboesch@gmail.com)