The brand-new NIntegrate

...and its new features

Anton Antonov

Wolfram Research, Inc.

Introduction

Integration rules: 1D and nD

Singularity handling

Monte Carlo integration

NIntegrate borrows design and implementational ideas from NDSolve

The talk is mostly examples of the new features

1D Integration Rules

Here is a list of the built-in NIntegrate one dimensional rules

{TrapezoidalRule,NewtonCotesRule,GaussKronrodRule, LobattoKronrodRule,ClenshawCurtisRule};

This is a way to use a rule

NIntegrate[(1/(x^2)),{x,0.4,1},Method->NewtonCotesRule]

... which is the shortened version of

NIntegrate[(1/(x^2)),{x,0.4,1},Method->{GlobalAdaptive,Method->NewtonCotesRule}]

Here is list of the sampling points and the order of the integrand evaluation over them for 1/(x^2)  over [0.4,1] with each one dimensional rule.

(t=Reap[NIntegrate[(1/(x^2)),{x,0.4,1},Method->#,EvaluationMonitor:>Sow[x]]][[2,1]];ListPlot[Transpose[{t,Range[Length[t]]}]])&/@{TrapezoidalRule,NewtonCotesRule,GaussKronrodRule, LobattoKronrodRule,ClenshawCurtisRule}

Single plots

nD Integration Rules

Fully-symmetric rules

NIntegrate has two new multidimensional rules based on the Invariant theory.

NIntegrate[Cos[x^2+y^2],{x,0,π},{y,0,π},Method->MultiDimensional]

Profiles

NIProfile@
NIntegrate[Cos[x^2+y^2],{x,0,π},{y,0,π},Method->MultiDimensional]

NIProfile@
NIntegrate[Cos[x^2+y^2],{x,0,π},{y,0,π},Method->{MultiDimensionalRule,Generators->9}]

2D generated rules

Graphics[{Blue,Point/@Reap[NIntegrate[x+y,{x,0,1},{y,0,1},Method->{MultiDimensionalRule,Generators->#},EvaluationMonitor:>Sow[{x,y}]]][[2,1]]},Frame->True,PlotRange->{{0,1},{0,1}}]&/@{5,6,9}

3D generated rules

Graphics3D[{Blue,Point/@Reap[NIntegrate[x+y+z,{x,0,1},{y,0,1},{z,0,1},Method->{MultiDimensionalRule,Generators->#},EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]]},PlotRange->{{0,1},{0,1},{0,1}}]&/@{5,6,9}

Cartesian rules

NIntegrate can do multi-dimensional integration with a Cartesian product of one or several 1D integration rules.

Plot3D[Exp[-(((1-x))^4) -(((1- y))^4)], {x, -2, 4}, {y, -2, 4}]

NIntegrate[Exp[-(((1-x))^4) -(((1- y))^4)], {x, -2, 4}, {y, -2, 4}, Method -> {GaussKronrodRule, NewtonCotesRule}]

NISamplingPoints@NIntegrate[(1/(Sqrt[x]))+y, {x, 0, π}, {y, 0, π}, Method -> {GaussKronrodRule, NewtonCotesRule}]/.Graphics[s_,op___]:>Graphics[{Blue,s},op]

tbl=Reap[NIntegrate[Sqrt[x]+Sqrt[y],{x,0,1},{y,0,1-x^2},Method->GaussKronrodRule,SingularityDepth->100,PrecisionGoal->5,EvaluationMonitor:>Sow[{x,y}]]][[2,1]];

dim=2;points=5;tblPart=Partition[tbl,(2*points+1)^dim];

Plot3D[Sqrt[x]+Sqrt[y],{x,(10^(-4)),1},{y,(10^(-4)),1},RegionFunction->(#2<1-#1^2&)]

Graphics3D[MapIndexed[MakePolygon[#1,#2[[1]]/Length[tblPart],2*points+1]&,tblPart],Axes->True,PlotRange->All]

Graphics3D[Point/@MapIndexed[Map[Function[{u},Append[u,#2[[1]]/Length[tblPart]]],#1]&,Take[tblPart,All]],Axes->True,PlotRange->All]

Clear[MakePolygon];
MakePolygon[points_,z_,step_:11]:=
Module[{t,t1},
t=Partition[points,step];
t1=Drop[#1,{2,-2}]&/@Take[t,{2,-2}];
t1=Transpose[t1];
t=Join[t[[1]],t1[[2]],Reverse[t[[-1]]],Reverse[t1[[1]]]];
Polygon[Append[#1,z]&/@t]
];

MakePolygon[Flatten[Table[{i,j},{i,3},{j,6}],1],0,6]

Clear[MakeRectangle];
MakeRectangle[points_,z_]:=
Module[{t,xmin,xmax,ymin,ymax},
t=#[[1]]&/@points;
{xmin,xmax}={Min[t],Max[t]};
t=#[[2]]&/@points;
{ymin,ymax}={Min[t],Max[t]};
Polygon[{{xmin,ymin,z},{xmin,ymax,z},{xmax,ymax,z},{xmax,ymin,z}}]
];

1D Singularity Handling

IMT tuning

NIProfileWithExact@
NIntegrate[(1/(Sqrt[1-x])),{x,0,1},Method->{GlobalAdaptive,SingularityHandler->{IMT,TuningParameters->1}},PrecisionGoal->12]

NIProfileWithExact@
NIntegrate[(1/(Sqrt[1-x])),{x,0,1},Method->{GlobalAdaptive,SingularityHandler->{IMT,TuningParameters->10}},PrecisionGoal->12]

DoubleExponential

NIProfileWithExact@
NIntegrate[(1/(Sqrt[1-x]))Log[(1/x)],{x,0,1},Method->{GlobalAdaptive,SingularityHandler->DoubleExponential},PrecisionGoal->12]

Automatic singularity handling

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

NISamplingPoints@
NIntegrate[(1/((10^4)((((1/32)-x))^2)+1/32)),{x,0,1},Method->{GlobalAdaptive,MaxRecursion->20,SingularityHandler->Automatic}]

nD Singularity Handling

IMT (old feature)

Duffy Coordinates

The integration strategy DuffyCoordinates gives quickly results for integrands of "spherical nature".

If dis the dimension of integration and r=((Subscript[x, 1])^2)+((Subscript[x, 2])^2)+...+((Subscript[x, d])^2), then Duffy's coordinates is a suitable technique for singularities of the type:
1. (r^α), (r^α) ln r, α>-d ;
2. ((Subscript[x, 1])^(Subscript[α, 1]))((Subscript[x, 2])^(Subscript[α, 2]))...((Subscript[x, d])^(Subscript[α, d]))(r^β), (Subscript[α, i])>-1, i∈[1,d], ∑(Subscript[α, i])+β>-d ;
3. ((((Subscript[c, 1])((Subscript[x, 1])^β)+(Subscript[c, 2])((Subscript[x, 2])^β)+...+(Subscript[c, d])((Subscript[x, d])^β)))^α), β>0, α β>-d, (Subscript[c, i])>0, i∈[1,d].

Motivation

NIntegrate[(1/(Sqrt[(x^2)+(((2y))^2)])),{x,0,1},{y,0,1},Method->{GlobalAdaptive, SingularityHandler->IMT}]//Timing

NIntegrate[(1/(Sqrt[(x^2)+(((2y))^2)])),{x,0,1},{y,0,1},Method->{GlobalAdaptive, SingularityDepth->∞,MaxRecursion->20}]//Timing

NIntegrate[(1/(Sqrt[(x^2)+(((2y))^2)])),{x,0,1},{y,0,1},Method->{DuffyCoordinates,Corners->{0,0}}]//Timing

ranges={"{x,0,1},{y,0,1}","{x,0,1},{y,0,x}","{x,0,1},{y,1-x,1}","{y,0,1},{x,0,y}"};
Manipulate[(k=0;t=NIntegrate[(1/(Sqrt[(x^2)+(((2 y))^2)])),Evaluate[Sequence@@ToExpression["{"<>r<>"}"]],EvaluationMonitor:>k++,SingularityDepth->∞,MaxRecursion->100];{k,t//InputForm}),{r,ranges}]

The moral is that we get results faster if we partition the region into simplexes

p1=Reap[NIntegrate[(1/(Sqrt[(x^2)+(((2 y))^2)])),{x,0,1},{y,0,x},Method->GaussKronrodRule,SingularityDepth->∞,EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
p2=Reap[NIntegrate[(1/(Sqrt[(x^2)+(((2 y))^2)])),{y,0,1},{x,0,y},Method->GaussKronrodRule,SingularityDepth->∞,EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
Graphics[{Blue,Point/@p1,Red,Point/@p2},Frame->True]

Explanation example

Append[#,AxesOrigin->{0,0}]&@NISamplingPoints@
NIntegrate[(1/(Sqrt[(x^2)+(y^2)])),{x,0,1},{y,0,1},Method->{DuffyCoordinates,Corners->{0,0}}]

NIntegrate`DuffyCoordinatesIntegrand[(1/(Sqrt[(x^2)+(y^2)])),{x,0,1},{y,0,1},Method->{DuffyCoordinates,Corners->{0,0}}]

Simplify[Expand[%],Assumptions->{x>0,y>0}]

Generalization

grx=GraphicsComplex[{{0,0,0},{1,0,0},{1,0,1},{1,1,1},{1,1,0}},{Polygon[{1,2,3}],Polygon[{1,3,4}],Polygon[{1,4,5}],Polygon[{1,5,2}],Polygon[{2,3,4,5}]}];
gry=MapAt[Map[RotateLeft[#]&,#]&,grx,{1}];
grz=MapAt[Map[RotateRight[#]&,#]&,grx,{1}];
Graphics3D[{Opacity[0.5],Red,grx,Cyan,gry,Yellow,grz}]

NIntegrate[(1/(Sqrt[(x^2)+(y^2)+(z^2)+(s^2)])),{x,0,1},{y,0,1},{z,0,1},{s,0,1},Method->{DuffyCoordinates,Corners->{0,0,0,0}}]//Timing

NIntegrate`DuffyCoordinatesIntegrand[(1/(Sqrt[(x^2)+(y^2)+(z^2)+(s^2)])),{x,0,1},{y,0,1},{z,0,1},{s,0,1},Method->{DuffyCoordinates,Corners->{0,0,0,0}}]

Automatic

Oscillatory Integration

NIntegrate automatically detects oscillatory (one dimensional) integrands, and automatically  decides which algorithm to use according to integrands' ranges.

The integrals detected as being of oscillatory type have the form

(∫_a^b)k[x] f[x]dx,

in which the oscillating kernel k(x) is of the form:
1. sin(ω (x^p)+c), cos(ω (x^p)+c), for (a,b) finite;
2. sin(ω (x^p)+c), cos(ω (x^p)+c), (Subscript[J, p])(x),or (Subscript[Y, p])(x),for (a,b) infinite or semi-infinite.
In these oscillating kernel forms ω, p,and c are constants, and p is an integer.

Finite range oscillatory integration

Comparison of NIntegrate and Integrate:

NIntegrate[Sin[200000*(x^2) + 50]Exp[x], {x, 1, 20}, 
  Method -> {GlobalAdaptive, MaxRecursion -> 20,SingularityDepth->∞}]//InputForm

Integrate[Sin[a*((x )^2)+ b]Exp[x], {x, 1, 20}]
%/.{a->200000,b->50}//N//InputForm

Automatic oscillatory integrand detection

Here is a piecewise oscillatory integration

NIntegrate[If[x<1,Cos[20000x](1/(Sqrt[x])),Sin[20000*x]*(1/(x^2))],{x,0,∞},Method->{GlobalAdaptive,MaxRecursion->20}]

Monte Carlo Integration

Adaptive Monte Carlo

2D example -- narrow peak

Plot3D[Exp[-((x^2)+(y^2))],{x,-5π,10π},{y,-5π,10π},PlotRange->All,PlotPoints->30]

NIntegrate[Exp[-((x^2)+(y^2))],{x,-5π,10π},{y,-5π,10π},Method->AdaptiveMonteCarlo]

exact=Integrate[Exp[-((x^2)+(y^2))],{x,-5π,10π},{y,-5π+4/5,10π}];exact//N

Mean[Table[NIntegrate[Exp[-((x^2)+(y^2))],{x,-5π,10π},{y,-5π,10π},Method->AdaptiveMonteCarlo],{n,40}]]

NISamplingPoints@NIntegrate[Exp[-((x^2)+(y^2))],{x,-5π,10π},{y,-5π,10π},Method->AdaptiveMonteCarlo]

3D example -- complicated boundary

ineqn=Abs[x]+Sqrt[Abs[y]]+Abs[z]<2/3;

ContourPlot3D[Boole[ineqn],{x,-1,1},{y,-1,1},{z,-1,1}]

pnts=Reap[NIntegrate[Boole[ineqn],{x,-1,1},{y,-1,1},{z,-1,1},SymbolicPreprocessing->None,Method->{AdaptiveMonteCarlo,BisectionDithering->1/10,Partitioning->5,MaxRecursion->20},PrecisionGoal->2,EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]];pnts//Length

Graphics3D[Point/@Take[pnts,{1,-1,4}],PlotRange->All,Axes->True,AxesLabel->{x,y,z}]

pntsIn=Select[pnts,ineqn/.Thread[{x,y,z}->#1]&];
Graphics3D[Point/@Take[pntsIn,{1,-1,1}],PlotRange->All,Axes->True,AxesLabel->{x,y,z}]

Adaptive Quasi Monte Carlo

NIntegrate[Exp[-((x^4)+(y^4))],{x,-2π+1/2,2π},{y,-2π+4/5,2π},Method->AdaptiveQuasiMonteCarlo]

NISamplingPoints@NIntegrate[Exp[-((x^4)+(y^4))],{x,-2π+1/2,2π},{y,-2π+4/5,2π},Method->AdaptiveQuasiMonteCarlo]

Symbolic Preprocessing

The capabilities of all strategies are extended through symbolic preprocessing of the integrals.The preprocessors can be seen as strategies that delegate integration to other strategies (preprocessors included).

Here are the (general) symbolic preprocessors:

1D example

(∫_(-1)^1)(Piecewise[{{(2/( (Sqrt[sin(-x)]))), x<0}, {(1/(Sqrt[x])), x≥0}}])dx.

pnts=Reap[NIntegrate[Piecewise[{{(2/(Sqrt[Sin[-x]])),x<0},{(1/(Sqrt[x])),x>= 0}}],{x,-1,1},PrecisionGoal->8,EvaluationMonitor:>Sow[x]]][[2,1]];
ListPlot[Transpose[{pnts,Range[Length[pnts]]}],PlotRange->All,AxesOrigin->{-1,0},AxesLabel->{x,ord}]

2D example

The integral (∫_(-π)^π)(∫_(-π)^π)Boole[(x^2)+(y^2)>1] sin((x^2)+(y^2))dydxhas been partitioned into(∫_(-π)^(-1))(∫_(-π)^π)sin((x^2)+(y^2))dydx   + (∫_(-1)^1)(∫_(-π)^(-(Sqrt[1-(x^2)])))sin((x^2)+(y^2))dydx + (∫_(-1)^1)(∫_(Sqrt[1-(x^2)])^π)sin((x^2)+(y^2))dydx  + (∫_1^π)(∫_(-π)^π)sin((x^2)+(y^2))dydx, that is why the sampling points form different pattern for-1≤x≤1.

gr=Plot3D[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π, π},DisplayFunction->Identity];
Off[NIntegrate::"slwcon"];grSP=NISamplingPoints[NIntegrate[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π,π},PrecisionGoal->2.5]][[1]];
On[NIntegrate::"slwcon"];
grSP=Graphics[grSP,Axes->True,AxesOrigin->{- π,- π}];
Show[GraphicsArray[{gr,grSP}]]

No even-odd handling

gr=Plot3D[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π, π},DisplayFunction->Identity];
Off[NIntegrate::"slwcon"];grSP=NISamplingPoints[NIntegrate[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π,π},Method->{SymbolicPiecewiseSubdivision,Method->GlobalAdaptive},PrecisionGoal->2.5]][[1]];
On[NIntegrate::"slwcon"];
grSP=Graphics[grSP,Axes->True,AxesOrigin->{- π,- π}];
Show[GraphicsArray[{gr,grSP}]]

Extensions

NIntegrate like NDSolve is a framework that speeds up new algorithm development.

Integration Rules

Simpson rule

NIntegrate`NewtonCotesRuleData[2,MachinePrecision]

(1/6.)-(1/2)

SimpsonRule/:NIntegrate`InitializeIntegrationRule[SimpsonRule,nfs_,ranges_,ruleOpts_,allOpts_]:=
NIntegrate`GeneralRule[{{0,(1/2),1},{(1/6),(4/6),(1/6)},{(1/6),(4/6),(1/6)}-{(1/2),0,(1/2)}}]

NIntegrate[(Sqrt[x]),{x,0,1},Method->SimpsonRule]

Multi-panel Simpson rule

Options[MultiPanelSimpsonRule]={"Panels"->5};
MultiPanelSimpsonRuleProperties=Part[Options[MultiPanelSimpsonRule],All,1];
MultiPanelSimpsonRule/:NIntegrate`InitializeIntegrationRule[MultiPanelSimpsonRule,nfs_,ranges_,ruleOpts_,allOpts_]:=
Module[{t,panels,pos,absc,weights,errweights},
t=NIntegrate`GetMethodOptionValues[MultiPanelSimpsonRule,MultiPanelSimpsonRuleProperties,ruleOpts];
If[t===$Failed,Return[$Failed]];
{panels}=t;

If[!TrueQ[NumberQ[panels]&&1≤panels<Infinity],pos=NIntegrate`OptionNamePosition[ruleOpts,"Panels"];
Message[NIntegrate::intpm,ruleOpts,{pos,2}];
Return[$Failed];
];
weights=Table[{(1/6),(4/6),(1/6)},{panels}];
weights=Fold[Join[Drop[#1,-1],{#1[[-1]]+#2[[1]]},Rest[#2]]&,First[weights],Rest[weights]]/panels;

{absc,errweights,t}=NIntegrate`TrapezoidalRuleData[((Length[weights]+1)/2),MachinePrecision];
NIntegrate`GeneralRule[{absc,weights,(weights-errweights)}]
]

NIntegrate[(Sqrt[x]),{x,0,1},Method->{MultiPanelSimpsonRule,Panels->20}]

NISamplingPoints@NIntegrate[(x^2),{x,0,1},Method->{MultiPanelSimpsonRule,Panels->100}]

Integration Strategies

Clear[DannysPlayGround];
Options[DannysPlayGround]={};
DannysPlayGroundProperties=Options[DannysPlayGround][[All,1]];
DannysPlayGround/:NIntegrate`InitializeIntegrationStrategy[DannysPlayGround,nfs_,ranges_,strOpts_,allOpts_]:=
Module[{},
DannysPlayGround[First/@ranges]
];
DannysPlayGround[vars_]["Algorithm"[regions_,opts___]]:=
Module[{ranges,funcs},
ranges=Map[Flatten/@Transpose[{vars,#@"Boundaries"}]&,regions];

ranges=Flatten/@ranges;funcs=(#@"Integrand"[])@"FunctionExpression"[]&/@regions;
Print[MapThread[Integrate[#1,#2]&,{funcs,ranges}]];
N[Total[MapThread[Integrate[#1,#2]&,{funcs,ranges}]],WorkingPrecision/.opts]
];

NIntegrate[(Sqrt[x]),{x,0,1},Method->DannysPlayGround]

NIntegrate[Piecewise[{{Sin[x],x<1},{(1/x),1<x<2}},(Sqrt[x])],{x,0,3},Method->DannysPlayGround]

Summary

Better incorporated symbolic preprocessing

New rules -- 1D and nD

New singularity handlers

Finite ranges oscillatory integration

Adaptive Monte Carlo integration

Extensibility

NIntegrate and NDSolve similarities