Math::NIntegrate

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:
- A strategy decorator has a strategy
- A strategy that is not a decorator, has a rule
- A composite rule (like Cartesian Rule) has one or more component rules
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