On NIntegrate's algorithms

New development and little known old features

Anton Antonov

Wolfram Research, Inc.

What is this talk about

The talk will present the improvements made in NIntegrate to handle piecewise integrands,and even or odd integrands over symmetric domains.The singularity flattening functionality will be explained. NIntegrate's DoubleExponent algorithm will be outlined. The usage of functional boundaries and Cartesian rules will be demonstrated for multi dimensional integrals.

How Global Adaptive algorithms NIntegrate work

New algorithms in 5.1

Piecewise integrands

Even or odd integrands

Singularity flattening

Double Exponent

Multi-dimensional integration topics

Functional boundaries

Cartesian rules

How Global Adaptive algorithms NIntegrate work

How it number crunches

Function

f[x_] := 
        Sech[100*(-4/10 + x)]^4 + 
          Sech[10*(-2/10 + x)]^2;
Plot[f[x], {x, 0, 1}]

Animation

Final interval configuration

\(grSin = 
      Import["/home/antonov/pics/KillSmthAndEatIt_NIntegrateSin.jpg", "JPEG"]\)

New algorithms in 5.1

NIntegrate was enhanced with handling of the piecewise, even, and odd functions.

Handling the piecewise functions makes NIntegrate faster since its rules are constructed for smooth integrands.

The even and odd integrand handling reduces the work of NIntegrate by half in each even/odd dimension, and simplifies the integration process.

Several new messages were added, the Method option takes nested option values

Piecewise Integrands

Algorithm

The PiecewiseNIntegrate algorithm works in the following way

1. After checking for correctness of the options given to NIntegrate, the variables are localized and the integrand is evaluated (if the global option "EvaluateNumericalFunctionArgument"  is set to True).

2. If the integrand has any of the piecewise functions

Piecewise,  Max, Min, UnitStep, Boole, Abs, Sign, If, Which, Switch, Floor, Ceiling, Round, IntegerPart, FractionalPart, Mod, Quotient

and the global option "EvaluateNumericalFunctionArgument"  is set to True, then the integration proceeds with the piecewise handling routine.

3. The piecewise handling routine calls NIntegrate for the integrants derived from the integrand f, over their corresponding disjoint regions derived from integration ranges limits.

4. Since PiecewiseNIntegrate might invoke quite a few NIntegrations only one message is given for each tag.

5. If PiecewiseNIntegrate does not return a number,  the integration proceeds using the NIntegrate method derived from the method option.

Method option

The option Method->Automatic is equivalent to Method->{SymbolicPiecewiseSubdivision, Method->GaussKronrod}:

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}]

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}, 
      Method -> Automatic]

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}, 
      Method -> {"SymbolicPiecewiseSubdivision", 
          Method -> GaussKronrod}]

The following two commands

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}, 
      Method -> SymbolicPiecewiseSubdivision]

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}, 
      Method -> {SymbolicPiecewiseSubdivision}]

are equivalent to

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}, 
      Method -> {"SymbolicPiecewiseSubdivision", 
          Method -> Automatic}]

Things we couldn't do (or did them hard)

GaussKronrod

We can do this integral with PiecewiseNIntegrate in V5.1

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}]

but we can't in V5.0.1:

NIntegrate[1/Sqrt[Abs[x]], {x, -1, 1}]

As the message says, the reason V5.0.1 to fail is that the GausKronrod rule has an abscissa on 0. In V6.0 we essentially have the sum of the two integrals:

NIntegrate[
        1/Sqrt[-x], {x, -1, 0}] + 
NIntegrate[
        1/Sqrt[x], {x, 0, 1}]

MonteCarlo

We can do this integral with the MonteCarlo method and PiecewiseNIntegrate

NIntegrate[
      If[x^2 + y^2 ≤ 1, 1, 0], {x, -1000, 
        1000}, {y, -1000, 1000}, 
      Method -> {"SymbolicPiecewiseSubdivision", 
          Method -> MonteCarlo[1000]}]

but we cannot with just MonteCarlo

NIntegrate[
      If[x^2 + y^2 ≤ 1, 1, 0], {x, -1000, 
        1000}, {y, -1000, 1000}, 
      Method -> MonteCarlo[1000]]

Even or Odd Integrands

Algorithm

1. The range of each variable (Subscript[x, i]) is tested is it symmetrical around the origin or not. If the range is symmetric then the integrand is tested is it an even or an odd function wrt (Subscript[x, i]).

2. If the integrand is an even or odd function wrt (Subscript[x, i]), and if the ranges of (Subscript[x, j]), j>i are even wrt (Subscript[x, i]), then the integration will be done only on the positive (right) half of (Subscript[x, i])' s range.

3. In case of odd integrand detected, if the numerical integration indicates lack of convergence, a message is printed that NIntegrate's result is zero, but the integral might be divergent.

4. Otherwise, the result is multiplied by 2^(number of even variables).for even integrands, or zero for odd ones

Method Options

The following integral is computed much faster using the EvenOddSubdivision option value:

NIntegrate[Sin[x^2] + Sin[y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
        Method -> {SymbolicPiecewiseSubdivision}] // Timing

NIntegrate[Sin[x^2] + Sin[y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
        Method -> {EvenOddSubdivision, 
            Method -> {SymbolicPiecewiseSubdivision}}] // 
      Timing

When it is very fast (Steinhaus tests)

Timing[
Print[
        NIntegrate[
          Sin[x^2] + Sin[y^2] + Sin[z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}]]; 
Print[
        NIntegrate[
          2^Sin[x] + 3^Cos[y] + 4^Sin[z], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}]]; 
Print[
        NIntegrate[
          2^Log[x^2 + y^2 + z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}]]; 
Print[
        NIntegrate[
          Sin[x^2 - y^2 - z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}]]; 
Print[
        NIntegrate[
          Cos[x^2 + y^2 + z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}]]; 
Print[
        NIntegrate[
          Sin[x^2 + y^2 + z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}]];
]

The integration is very fast using Cartesian product rules.

Timing[
      Print[NIntegrate[
          Sin[x^2] + Sin[y^2] + Sin[z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}, 
          Method -> {EvenOddSubdivision, 
              Method -> GaussKronrod}]]; 
Print[
        NIntegrate[
          2^Sin[x] + 3^Cos[y] + 4^Sin[z], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}, 
          Method -> {EvenOddSubdivision, 
              Method -> GaussKronrod}]]; 
Print[
        NIntegrate[
          2^Log[x^2 + y^2 + z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}, 
          Method -> {EvenOddSubdivision, 
              Method -> GaussKronrod}]]; 
Print[
        NIntegrate[
          Sin[x^2 - y^2 - z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}, 
          Method -> {EvenOddSubdivision, 
              Method -> GaussKronrod}]]; 
Print[
        NIntegrate[
          Cos[x^2 + y^2 + z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}, 
          Method -> {EvenOddSubdivision, 
              Method -> GaussKronrod}]]; 
Print[
        NIntegrate[
          Sin[x^2 + y^2 + z^2], {x, -Pi, Pi}, {y, -Pi, 
            Pi}, {z, -Pi, Pi}, 
          Method -> {EvenOddSubdivision, 
              Method -> GaussKronrod}]];
]

When it should be suspected

NIntegrate[x, {x, -∞, ∞}]

Example

NIntegrate[
      x, {x, -1, 1}, {y, -x, x}, {z, -1, 1}]

Example Plots

t = Reap[
            NIntegrate[x^2, {x, -1, 1}, {y, 0, 1 + x}, 
              MinRecursion -> 3, MaxRecursion -> 3, 
              EvaluationMonitor :> Sow[{x, y}]]][[2, 
            1]];

Show[Graphics[Point /@ t], Axes -> True, 
      AxesLabel -> {"x", "y"}]

Singularity flattening

The IMT Rule (the Iri, Moriguti, Takasawa rule) changes the integration variable

Clear[ϕ]
ϕ[{xmin_, xmax_}, t_] := 
        xmin + (xmax - 
                xmin) (E^(1 - (1/t)));
ϕ[t_] := ϕ[{0, 1}, t]

Plot[ϕ[t], {t, 0, 1}, PlotRange -> All, 
      AxesOrigin -> {-0.05, -0.05}]

Plot[Evaluate[D[ϕ[t], t]], {t, 0, 1}, 
      PlotRange -> All]

NIntegrate puts the singularities at the left side of the interval

f[x_] := 1/(x^2)

Plot[f[x], {x, 0.001, 0.05}, PlotRange -> All];

{a, b} = {0.001, 0.05};
Plot[
        f[ϕ[{a, b}, t]] D[ϕ[{a, b}, t], t] // Evaluate, {t, 0, 1}, 
        PlotRange -> All];

The singularity is smashed! But ...

This comes with a price to pay -- the evaluation points become too close to the singular end.

Double Exponent: Outline

The DoubleExponent algorithm was proposed by Mori and Takashi in 1974. It is inspired from the IMT and TANH rules.

Transformation

Can we transform the integrant in such a way that f[ϕ[t]] ϕ'[
            t] ≈ Exp[-c Exp[|t|]] as \(TraditionalForm\`\(\(|\)\(t\)\(|\)\) -> \(±∞\)\)?

ϕ[
        t_] := (1/2) (Tanh[(1/2) π Sinh[t]] + 1)

Plot[ϕ[t], {t, -3, 3}, PlotRange -> All];

Plot[Evaluate[D[ϕ[t], t]], {t, -3, 3}, 
        PlotRange -> All];

Implementation

\(f[x_] := 1\/\@x\%4\)

Integrate[f[x], {x, 0, 1}] // N // InputForm

{a, b} = {-4, 4};
n = 100;
h = ((b - a)/(n - 1));
s = 
        h*Total@Table[
              f[ϕ[t]]*D[ϕ[t], t] /. t -> a + i*h, {i, 1, 
                n - 2}];

N[s, 30] // InputForm

N[Abs[s - 4/3], 30]

Transformed Integrand

res = 
        Table[{x, f[ϕ[t]]*D[ϕ[t], t] /. t -> x}, {x, -3, 
            3, 0.1}];
ListPlot[res, PlotRange -> All, 
        PlotJoined -> True];

Convergence rate

<< Graphics`MultipleListPlot`

$MaxExtraPrecision  =  200. ;

err = 
Map[({a, b} = {-5, 
                5}; 
n = #1; 
h = ((b - a)/(n - 1)); 
s = 
              h*Total@Table[
                    f[ϕ[t]]*D[ϕ[t], t] /. t -> a + i*h, {i, 1, 
                      n - 2}]; 
N[Abs[s - 4/3], 
              40]) &, 
Range[10, 100, 10]]

c = 5; b = 20;
terr = 
        Map[Exp[-c #(1/Log)[#1] + b] &, 
          Range[10, 100, 10]];

MultipleListPlot[{Log@err, Log@terr}, PlotJoined -> True, 
      PlotStyle -> {Hue[0], Hue[0.7]}, 
      PlotLegend -> {"Computed\nError", "Model\nError"}]

Double Exponent: Examples

\(TraditionalForm\`∫\_0\%1\( 1\/\@x\%4\ \(log(
            1\/x)\)\) dx\)

Plot[(x^(-1/4)) Log[(1/x)], {x, 0, 1}]

exact = Integrate[(x^(-1/4)) Log[(1/x)], {x, 0, 1}] // N // 
        InputForm

NIProfile[NIntegrate[(x^(-1/4)) Log[(1/x)], {x, 0, 1}]]

NIProfile[
      NIntegrate[(x^(-1/4)) Log[(1/x)], {x, 0, 1}, 
        Method -> DoubleExponential]]

(∫_0^π) cos(
          64 sin(x)) dx

Plot[Cos[64*Sin[x]], {x, 0, Pi}]

Integrate[Cos[64*Sin[x]], {x, 0, Pi}]//N//InputForm

NIProfile[NIntegrate[Cos[64*Sin[x]], {x, 0, Pi}]]

NIProfile[
      NIntegrate[Cos[64*Sin[x]], {x, 0, Pi}, 
        Method -> DoubleExponential]]

\(TraditionalForm\`∫\_0\%1\( e\^\(20\ \((x - 1)\)\)\ \(sin(256\ x)\)\) dx\)

Plot[Exp[20 (x - 1)] Sin[256 x], {x, 0, 1}, 
      PlotRange -> All]

Integrate[Exp[20 (x - 1)] Sin[256 x], {x, 0, 1}] // N // 
      InputForm

NIProfile[
      NIntegrate[Exp[20 (x - 1)] Sin[256 x], {x, 0, 1}]]

NIProfile[
      NIntegrate[Exp[20 (x - 1)] Sin[256 x], {x, 0, 1}, 
        Method -> DoubleExponential]]

(∫_0^π) sin(
          3 x) dx

Plot[Abs[Sin[3*x]], {x, 0, π}]

Integrate[Abs[Sin[3*x]], {x, 0, π}]

NIProfile[
      NIntegrate[Abs[Sin[3*x]], {x, 0, π}, MaxRecursion -> 10, 
        Method -> GaussKronrod], 100]

NIProfile[
      NIntegrate[Abs[Sin[3*x]], {x, 0, π}, MaxRecursion -> 10, 
        Method -> DoubleExponential], 1]

Evaluantion points plots

t = {};
NIntegrate[Abs[Sin[3*x]], {x, 0, π}, MaxRecursion -> 20, 
      Method -> GaussKronrod, 
      EvaluationMonitor :> AppendTo[t, x]]

ListPlot[Transpose[{t, Range[Length[t]]}], 
      PlotRange -> All]

t = {};
NIntegrate[Abs[Sin[3*x]], {x, 0, π}, MaxRecursion -> 20, 
      Method -> DoubleExponential, 
      EvaluationMonitor :> (t = {t, 
              x})]
t = Flatten[t];

ListPlot[Transpose[{t, Range[Length[t]]}], 
      PlotRange -> All]

(∫_0^1) 
            (x^30) dx

Plot[x^30, {x, 0, 1}, PlotRange -> All]

NIProfile[
      NIntegrate[x^30, {x, 0, 1}, MaxRecursion -> 10, 
        Method -> GaussKronrod], 100]

NIProfile[
      NIntegrate[x^30, {x, 0, 1}, MaxRecursion -> 10, 
        Method -> DoubleExponential], 100]

(∫_0^1) 
            (1/(16 (((x - (π/4)))^2) + 
              (1/16))) dx

Plot[(1/(16 (((x - (π/4)))^2) + (1/16))), {x, 0, 1}, 
      PlotRange -> All]

Integrate[f[x], {x, 0, 1}] // N // InputForm

NIProfile[
      NIntegrate[f[x], {x, 0, 1}, Method -> GaussKronrod]]

de = 
      NIProfile[
        NIntegrate[f[x], {x, 0, 1}, Method -> DoubleExponential], 
        1]

t = {};
NIntegrate[f[x], {x, 0, 1}, MaxRecursion -> 20, 
      Method -> GaussKronrod, 
      EvaluationMonitor :> AppendTo[t, x]]

Evaluation plots points

ListPlot[Transpose[{t, Range[Length[t]]}], 
      PlotRange -> All]

t = {};
NIntegrate[f[x], {x, 0, 1}, MaxRecursion -> 20, 
      Method -> DoubleExponential, 
      EvaluationMonitor :> (t = {t, 
              x})]
t = Flatten[t];

ListPlot[Transpose[{t, Range[Length[t]]}], 
      PlotRange -> All]

NIProfile Definition

Clear[NIProfile]
SetAttributes[NIProfile, HoldFirst];
NIProfile[expr_, 
        n_:  100] := 
Module[{k, res, 
          t}, 
k = 0; 
res = 
          Hold[expr] /. {HoldPattern[NIntegrate[s___]] :> 
                NIntegrate[s, EvaluationMonitor :> k++], 
              HoldPattern[NewNIntegrate[s___]] :> 
                NewNIntegrate[s, 
                  EvaluationMonitor :> k++]}; 
res =  ReleaseHold[res]; 
t = 
          Hold[Timing[Do[expr, {i, 1, n}]]]; 
t = 
          ReleaseHold[t]; 
{FullForm[res], k,  
          t[[1]]/n // N}
]

Functional Boundaries for Multi-Dimensional Integration

t = Reap[
          NIntegrate[
            1, {y, 0, 4}, {z, ((y^2)/4), 4}, {x, 0, (Sqrt[4 z - (y^2)])}, 
            Method -> MultiDimensional, 
            EvaluationMonitor :> 
              Sow[{x, y, z}]]][[2, 
        1]];

Show[Graphics3D[Point /@ t]]

t = Reap[
            TimeConstrained[
              NIntegrate[
                1, {y, 0, 4}, {z, ((y^2)/4), 4}, {x, 0, (Sqrt[4 z - (y^2)])}, 
                Method -> GaussKronrod, 
                EvaluationMonitor :> Sow[{x, y, z}]], 
              0.2]][[2, 
          1]];

Show[Graphics3D[Point /@ t]]

Cartesian Rules for Multi-Dimensional Integration

Steinhaus related

k = 0;
NIntegrate[Cos[x^2 + y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
        Method -> {EvenOddSubdivision, Method -> GaussKronrod}, 
        EvaluationMonitor :> k++];
k

k = 0;
NIntegrate[Cos[x^2 + y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
        Method -> GaussKronrod, 
        EvaluationMonitor :> k++];
k

MultiDimensional

k = 0;
NIntegrate[Cos[x^2 + y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
        Method -> {EvenOddSubdivision, Method -> MultiDimensional}, 
        EvaluationMonitor :> k++];
k

k = 0;
NIntegrate[Cos[x^2 + y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
        Method -> MultiDimensional, 
        EvaluationMonitor :> k++];
k

Plot3D[Cos[x^2 + y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
      PlotPoints -> 100, PlotRange -> All]

Plot3D[Cos[x^2 + y^2], {x, -Pi, Pi}, {y, -Pi, Pi}, 
      PlotPoints -> 100, PlotRange -> All, 
      ViewPoint -> {0.000,  0.000,  3.260}]

t = {};
NIntegrate[Cos[x^2 + y^2], {x, 0, Pi}, {y, 0, Pi}, 
        Method -> MultiDimensional, 
        EvaluationMonitor :> (t = {t, {x, 
                  y}})];
tMD = Partition[Flatten[t], 2];

ListPlot[tMD, AspectRatio -> Automatic]

t = {};
NIntegrate[Cos[x^2 + y^2], {x, 0, Pi}, {y, 0, Pi}, 
        Method -> GaussKronrod, 
        EvaluationMonitor :> (t = {t, {x, 
                  y}})];
tGK = Partition[Flatten[t], 2];

ListPlot[tGK, AspectRatio -> Automatic]

Number of evaluations for x^p + y^p + z^p

<< Graphics`MultipleListPlot`

n = 15;

tGK = Table[(k = 0; NIntegrate[x^p + y^p + z^p, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Method -> GaussKronrod, EvaluationMonitor :> k++]; k), {p, 0, n}];

tMD = Table[(k = 0; NIntegrate[x^p + y^p + z^p, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Method -> MultiDimensional,EvaluationMonitor :> k++]; k), {p, 0, n}];

tGK

tMD

TableForm[Transpose[{tGK, tMD}], TableHeadings -> {Range[0, n], None}]

MultipleListPlot[Transpose[{Range[0, Length[tGK] - 1], tGK}], 
      Transpose[{Range[0, Length[tMD] - 1], tMD}], PlotJoined -> True, 
      PlotRange -> All]

Summary

PiecewiseNIntegrate an EvenOddNIntegrate leverage the symbolic capabilities of Mathematica

DoubleExponential is suitable for integrals with end singularities

GaussKronrod Cartesian rules a valuable alternative of the MultiDimensional rule for low dimensions (2 or 3)