Rand Stats

Math::NIntegrate

zef:antononcube

Math::NIntegrate

MacOS Linux Win64

Introduction

This repository has the Raku code of a library for numerical integration. The design and architecture of library's main function, nintegrate, resembles that of Wolfram Language NIntegrate, [WRI1, WRI2].


Installation

From Zef ecosystem:

zef install Math::NIntegrate

From GitHub:

zef install https://github.com/antononcube/Raku-Math-NIntegrate.git

Motivation

While working for Wolfram Research Inc., during the years 2003-2006 I designed and implemented Mathematica's "new" NIntegrate and extensively documented it. (I wrote [WRI2].)

The Raku programming language has some built-in features that give the ability to:

These are also programming language features on which Wolfram Language's NIntegrate is based upon.

Hence, it seems natural to think that an implementation of a powerful numerical integration framework in Raku that has unique features is achievable, worthy, and rewarding.

I plan to write a book that describes in detail software architecture designs and decisions for the development of numerical integration frameworks. That book is going to use this library.


Design (in brief)

The NIntegrate framework is based on Object-Oriented Programming (OOP) Design Patterns (DP). NIntegrate uses the software design patterns Strategy, Composite, Decorator, and others.

Here is a description using the OOP DP language:

  • The basic objects of NIntegrate are integration regions.
  • Each region has its own integration function and integration rule.
  • Conceptually, there are two main types of algorithms: integration strategies and integration rules.
  • The integration strategies use the Template method for their "logic."
  • The integration regions use Strategy for the computation of integral and error estimates.
  • The integration regions can utilize singularity handler objects.
  • Creations of integration rules generally use Builder.
  • Symbolic preprocessing is done through Decorator.
  • User specifications are translated into creations of numerical integration algorithm objects through Interpreter.
  • The creation of the "final" integration algorithm object uses Abstract factory.

Summary of (current) algorithms

The following mind-map shows the set of algorithms in "Math::NIntegrate". The currently available ones are marked with a green-square checkmark; the planned ones are marked with a red empty circle. (More details can found in "WorkPlan.org".)

mindmap
  root(**Math::NIntegrate**)
    Integration Rules
      ✅ _General rule_
      **1D**  
        ✅ Clenshaw-Curtis
        ✅ Gauss-Kronrod 
        ⭕️ Lobatto-Kronrod rule
        ✅ Newton-Cotes rule
        ✅ Trapezoidal rule
        ⭕️ Multi-panel rule
      **nD**  
        ✅ Cartesian rule
        ✅ MultiDimensional rule
        ⭕️ DoubleExponential rule
        ✅ Monte Carlo rule
        ✅ Quasi Monte Carlo rule
        ✅ _User defined rule_

    Integration Strategies
      ✅ Global Adaptive
      ✅ Local Adaptive
      ⭕️ Trapezoidal
      ⭕️ Double Exponential
      ✅ Monte Carlo
      ✅ Quasi Monte Carlo
      ✅ Adaptive Monte Carlo
      ✅ Adaptive Quasi Monte Carlo
      ⭕️ Importance Sampling
      ✅ _User defined strategy_

    Integration Strategy Decorators
      ✅ Unit Cube Rescaling
      ⭕️ Even-Odd Subdivision
      ⭕️ Symbolic Piecewise Subdivision

The algorithms are nested depending on their type:

Here is an example in which the method specification illustrates the nesting:

use Math::NIntegrate;

nintegrate( { exp(-($^x ** 2 + $^y ** 2)) }, <x -Inf Inf>, <y 2 Inf>, 
    method => ('unit-cube-rescaling', 
               method => ('global-adaptive', 
                          method => ('cartesian-rule', 
                                     method => ('gauss-kronrod-rule', :7points)
                          )
               )
    ),
    :4precision-goal 
)
# 0.007347769817290595

The framework works with method specification "shortcuts" -- the above method specification can be written as:

method => ('unit-cube-rescaling', method => ('gauss-kronrod-rule', :7points))

That is because the global adaptive strategy is the default strategy, and for multidimensional integrals Cartesian rule is automatically used if just a one-dimensional integration rule is specified. (Like, the Gauss-Kronrod rule above.)


Usage examples

Basic usage examples:

use Math::NIntegrate;

nintegrate( -> $x { 1/sqrt($x) }, ['x', 0, 2] );
# 2.8284271225432907

With the adverb "pairs" the result shows the integral and error estimates together with the number of subregions used:

nintegrate( -> $x { 1 / $x ** 2 }, <x 1 2>, working-precision => Rat ):pairs;
# {error => 1.7652322865675113e-07, integral => 0.5000000000000238, region-count => 1}

Compute a two-dimensional integral:

nintegrate( -> $x, $y { $x + $y² }, <x 0 2>, <y 0 12>);
# 1176.0000000000007

Adaptive Monte-Carlo integration (not implemented yet):

nintegrate( -> $x, $y, $z { $x + $y² + 1/$z⁻¹ }, x => (0, 2), y => (0, 12), z => (1, 4), method => 'adaptive-monte-carlo' );

Integration with functional boundaries (using a 2D Cartesian integration rule specified with a 1D component):

nintegrate( -> *@_ { 1 }, ['x', 0, 1], ['y', 0, {$^x.sqrt}], method => 'clenshaw-curtis-rule'):pairs
# {error => 6.46800288251937e-07, integral => 0.6666666120967495, region-count => 9}

Compute a two-dimensional integral with using a full specification of a Cartesian integration rule which based on a one-dimensional Gauss-Kronrod rule with 5 Gauss points and a one-dimensional Trapezoidal rule with 10 base-points:

nintegrate( { ($^x + $^y ** 2).sqrt }, <x 0 1>, <y 0 1>, 
        method => ('global-adaptive', 
                   method => ('cartesian-rule',
                              method => (('gauss-kronrod-rule', points => 5), ('trapezoidal-rule', points => 10))
                   )
        )):pairs
# {error => 4.208261570325412e-07, integral => 0.8786346480633286, region-count => 8}

Compare with the Wolfram Language results:

wolframscript -code 'Through[{Integrate, N@*Integrate, NIntegrate}[Sqrt[x+y^2], {x, 0, 1}, {y, 0, 1}]]'
# {(-2 + 7*Sqrt[2] + 3*ArcSinh[1])/12, 0.8786346414725245, 0.8786346981741846}

Using natural language specifications

Natural language commands for computational integration workflows can be translated to Raku commands via:

Here is full computational pipeline:

use DSL::English::IntegrationWorkflows;
use Math::NIntegrate;

my $cmd = 'integrate 1/x^2 over the range [1,2] with a local adaptive strategy and precision goal 5';

ToIntegrationWorkflowCode($cmd, 'Raku')
        ==> { say "interpreted: $_"; $_ }()
        ==> {
  use MONKEY-SEE-NO-EVAL;
  my $res = $_.EVAL;
  say "result: {$! ?? 'N/A' !! $res}"
}()
# interpreted: nintegrate({1/$^x **2}, ['x',1,2], method => "local-adaptive", precision-goal => 5)
# result: 0.5000000000000238

CLI

The package provides a Command Line Interface (CLI) script. Here is its usage message:

nintegrate --help
# Usage: nintegrate INTEGRAND RANGE-SPEC ... [OPTIONS]
# 
# The first positional argument is Raku code describing the integrand.
# 
# The remaining positional arguments are Raku code range specifications.
# 
# Examples of integrands:
#   '{$^x ** 2}'
#   'sub ($x, $y) { $x + $y }'
# 
# Examples of range specifications:
#   '<x 0 1>'
#   '["x", -2, 20]'
# 
# Named options are passed to Math::NIntegrate, including:
#   --precision-goal=VALUE
#   --method=VALUE
#   --max-recursion=VALUE
#   --OPTION=VALUE
# 
# Examples:
#   program '{$^x ** 2}' '<x 0 1>' --precision-goal=10
#   program 'sub ($x, $y) { $x + $y }' '["x", 0, 1]' '["y", 0, 2]' --method=adaptive
# 
# Use --help to display this message.

Here is an example invocation:

nintegrate '{$^x + $^y}' '<x 0 1>' '<y 0 10>' --pairs
# {error => 1.918654139689546e-14, integral => 54.99999999999999, region-count => 1}

Status

The implementation status is tracked in the Org-mode file Work-plan.org.


References

[WRI1] Wolfram Research (1988), NIntegrate, Wolfram Language function, https://reference.wolfram.com/language/ref/NIntegrate.html (updated 2014).

[WRI2] Wolfram Research (2006), Advanced Numerical Integration in the Wolfram Language, Wolfram Monograph, https://reference.wolfram.com/language/tutorial/NIntegrateOverview.html.

[AAr1] Anton Antonov, NIntegrate - The Missing Manual, (2019), GitHub/antononcube.


Anton Antonov
Windermere, Florida, USA
2021-04-05
2026-09-28