The London Perl and Raku Workshop takes place on 26th Oct 2024. If your company depends on Perl, please consider sponsoring and/or attending.


Math::Integral::Romberg - Scalar numerical integration


    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);


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



A reference to the function to be integrated.

$x1, $x2

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


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.


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


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.


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.


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.


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.


Eric Boesch (