NIntegrate Algorithms

Anton Antonov

Wolfram Research, Inc.

Integration strategies

({{strategies, deterministic, Monte Carlo}, {adaptive  , GlobalAdaptive,
LocalAdaptive, AdaptiveMonteCarlo,
AdaptiveQuasiMonteCarlo}, {non-adaptive , DoubleExponential,
Trapezoidal, MonteCarlo,
QuasiMonteCarlo}, {specialized, DuffyCoordinates,
Oscillatory, -}})

pnts=Reap[NIntegrate[1/(Sqrt[x]),{x,0,1},Method->"GlobalAdaptive",EvaluationMonitor:>Sow[x]]][[2,1]];ListPlot[Transpose[{pnts,Range[Length[pnts]]}]]

One-dimensional integration rules

funcs={(Sqrt[x]),Abs[x-(1/E)],(Piecewise[{{2, x≤(1/E)}, {3, x>(1/E)}}]),(1/((10^3) ((((1/π)-x))^2)+1))};

Plot[#,{x,0,1},AxesOrigin->{0,0},PlotRange->All,Frame->True,Filling->Axis]&/@funcs

fx=funcs[[1]];TableForm[{(k=0;NIntegrate[fx,{x,0,1},Method->{"GlobalAdaptive",Method->#,"SingularityDepth"->∞,MaxRecursion->100,"SymbolicProcessing"->0},EvaluationMonitor:>k++]),k}&/@{"GaussKronrodRule","GaussBerntsenEspelidRule","LobattoKronrodRule","LobattoPeanoRule","ClenshawCurtisRule","MultiPanelRule","NewtonCotesRule","TrapezoidalRule"}]

Mutli-dimensional 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->"MultiDimensionalRule"]

Profiles

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

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.

Consider the function

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

Integration with a Cartesian product rule

NIntegrate[Sqrt[x]+Sqrt[y],{x,0,1},{y,0,1-x^2},Method -> {"GaussKronrodRule", "NewtonCotesRule"}]

Integration visualization

Consider the function:

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

Sampling points of an integration with a Cartesian rule

tbl=Reap[NIntegrate[Sqrt[x]+Sqrt[y],{x,0,1},{y,0,1-x^2},Method->{"GlobalAdaptive",Method->"GaussKronrodRule","SingularityDepth"->100},PrecisionGoal->5,EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
dim=2;points=5;tblPart=Partition[tbl,(2*points+1)^dim];

Visualization of the integration process

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

Visualization of the integration process

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

Visualization code

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}}]
];

One dimensional 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}]

Multi-dimensional 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

Results and number of evaluation points for different ranges for the integrand (1/(Sqrt[(x^2)+(((2 y))^2)])):

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}"};
Table[{r,(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}]//TableForm

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

Sampling points of an example integration

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 strategies

NIntegrate's oscillatory strategies are for one dimensional integrals.

The algorithms for finite region integrals are different from the algorithms for infinite regions.

NIntegrate uses Chebyshev expansions of the integrand and the global adaptive integration strategy for finite region oscillatory integrals.

For infinite oscillatory integrals NIntegrate uses either 
-- a modification of the double exponential quadrature or
-- sequence summation acceleration over the sequence of integrals with regions between the integrand's zeroes.

Here is an example that all three algorithms:

NIntegrate[(Piecewise[{{((Cos[2000 x])/(Sqrt[x])), 0<x<2}, {((Sin[20 x])/(Sqrt[(x^2)])), x<0}, {BesselY[2,(x^3)]/x, x>2}}]),{x,-∞,∞}]

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

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), (e^(i ω (x^p))) for (a,b) finite;
2. sin(ω (x^p)+c), cos(ω (x^p)+c), (e^(i ω (x^p))), (Subscript[J, ν])(ω (x^p)+c), (Subscript[Y, ν])(ω (x^p)+c), ((Subscript[H, ν])^((1)))(ω (x^p)+c) , ((Subscript[H, ν])^((2)))(ω (x^p)+c), (Subscript[j, ν])(ω (x^p)+c), or (Subscript[y, ν])(ω (x^p)+c) for (a,b) infinite or semi-infinite.
In these oscillating kernel forms ω, c and ν are real constants, and p is a positive integer.

f(x) is assumed to be algebraic, but in many cases NIntegrate works well on non-algebraic f(x).

Finite oscillatory integration

Modified Clenshaw-Curtis quadrature is for finite region one dimensional integrals of the form

(∫_a^b)sin (ω (x^p)+c) f(x)dx  ,  (∫_a^b)cos(ω (x^p)+c) f(x)dx, or(∫_a^b)exp(i ω (x^p)+c) f(x)dx,

where a, b, k, c, p are finite real numbers.

Highly oscillatory integrand

NIntegrate[((Sin[2000000 (x^2)+50000])/(x^2)),{x,(2/10),2}]//Timing

Oscillatory integration with singularity

NIntegrate[((Sin[20000 x])/(Sqrt[x])),{x,0,2},MaxRecursion->20]//Timing

Sampling points of oscillatory integration with singularity

Manipulate[
pnts=Reap[NIntegrate[((Sin[20000 x])/(Sqrt[x])),{x,a,2},MaxRecursion->20,EvaluationMonitor:>Sow[x]]][[2,1]];
ListPlot[Transpose[{pnts,Range[Length[pnts]]}],PlotRange->{{0,2},{0,700}}]
,{a,0,0.1,0.01}]

Infinite oscillatory integration

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:
sin(ω (x^p)+c), cos(ω (x^p)+c), (e^(i ω (x^p))), (Subscript[J, ν])(ω (x^p)+c), (Subscript[Y, ν])(ω (x^p)+c), ((Subscript[H, ν])^((1)))(ω (x^p)+c) , ((Subscript[H, ν])^((2)))(ω (x^p)+c), (Subscript[j, ν])(ω (x^p)+c), or (Subscript[y, ν])(ω (x^p)+c), for (a,b) infinite or semi-infinite.
In these oscillating kernel forms ω, c and ν are real constants, and p is a positive integer.

{Sin[w x^p+c],Cos[w x^p+c],Exp[I w x^p],BesselJ[n,w x^p+c],BesselY[n,w x^p+c],SphericalBesselJ[n,w x^p+c],SphericalBesselY[n,w x^p+c],HankelH1[n,w x^p+c],HankelH2[n,w x^p+c]}

NIntegrate[BesselJ[2,(x^3)+20],{x,1,∞},WorkingPrecision->30,PrecisionGoal->20]

Using ExtrapolatingOscillatory with different methods

Manipulate[
pnts=Reap[NIntegrate[BesselJ[2,(x^2)+20],{x,0,∞},MaxRecursion->20,EvaluationMonitor:>Sow[x],Method->{"ExtrapolatingOscillatory",Method->method}]][[2,1]];
ListPlot[Transpose[{pnts,Range[Length[pnts]]}],PlotRange->All]
,{method,{"GlobalAdaptive","LocalAdaptive","DoubleExponential"}}]

Parsing sums of oscillatory functions

In some cases it might better if sums of trigonometric functions are treated with the oscillatory alforithms

SphericalBesselJ[2,x]/x//FunctionExpand//Simplify

NIntegrate[SphericalBesselJ[2,x]//FunctionExpand,{x,1,∞},Method->{"OscillatorySelection","TermwiseOscillatory"->True}]

NIntegrate[SphericalBesselJ[2,x]//FunctionExpand,{x,1,∞}]

Application: Fourier series calculation

This is a Mathematica function that calculates a truncated Fourier series approximation of a function.

FourierAnalysis[f_,{x_,xmin_,xmax_},nterms_,integrator_:(NIntegrate[##,Method->{"GlobalAdaptive","MaxRecursion"->30}]&)]:=
FourierAnalysis[f,{x,xmin,xmax},nterms,True,True,integrator];

FourierSinAnalysis[f_,{x_,xmin_,xmax_},nterms_,integrator_:(NIntegrate[##,Method->{"GlobalAdaptive","MaxRecursion"->30}]&)]:=
FourierAnalysis[f,{x,xmin,xmax},nterms,True,False,integrator];

FourierCosAnalysis[f_,{x_,xmin_,xmax_},nterms_,integrator_:(NIntegrate[##,Method->{"GlobalAdaptive","MaxRecursion"->30}]&)]:=
FourierAnalysis[f,{x,xmin,xmax},nterms,False,True,integrator];

FourierAnalysis[f_,{x_,xmin_,xmax_},nterms_,ds:True|False,dc:True|False,integrator_]:=Block[{a,b,funcTerms},a=(2/((xmax-xmin)))Table[integrator[Cos[((2π )/(xmax-xmin))j x] f,{x,xmin,xmax}],{j,0,If[dc,nterms,0]}];
If[!dc,a=Join[a,Table[0,{j,0,nterms}]]];
b=If[ds,(2/((xmax-xmin)))Table[integrator[Sin[((2π )/(xmax-xmin))j x] f,{x,xmin,xmax}],{j,1,nterms}],Table[0,{j,0,nterms}]];
funcTerms=((a[[1]])/2)+Total[Table[Cos[((2π )/(xmax-xmin))j x] a[[j+1]]+Sin[((2π )/(xmax-xmin))j x] b[[j]],{j,1,nterms}]];
funcTerms
];

This calculates a 60 terms Fourier approximation of Sin[(x^3)+(1/2)] over [-4,4] using NIntegrate. If Integrate is used the calculation will be very slow.

func=FourierAnalysis[Sin[(x^3)+(1/2)],{x,-4,4},60];//Timing

This a plot of Sin[(x^3)+(1/2)] and the Fourier series approximation calculated above.

Plot[{Tooltip[Sin[(x^3)+(1/2)],"Original\nfunction"],Tooltip[func,"Fourier\napproximation"]},{x,-4,4},PlotRange->All]

This calculates a 40 terms Fourier approximation of x-⌊x⌋ over [0,1] using NIntegrate. If Integrate is used the calculation will be very slow.

func=FourierSinAnalysis[x-Floor[x],{x,0,1},4000,NIntegrate];//Timing

This a plot of x-⌊x⌋ and the Fourier series approximation calculated above.

Plot[{Tooltip[x-Floor[x],"Original\nfunction"],Tooltip[func,"Fourier\napproximation"]},{x,0,2},PlotRange->All,Filling->{1->Axis},PlotStyle->{Thickness[0.003],Thickness[0.0015]}]

Monte Carlo integration

X=Array[x,{30}];

Timing[NIntegrate[Evaluate[1/(1+Norm[X])],Evaluate[Sequence@@({#1,-1,1}&)/@X],Method->"AdaptiveQuasiMonteCarlo"]]

X=Array[x,{10^3}];
Timing[NIntegrate[Evaluate[1/(1+Norm[X])],Evaluate[Sequence@@({#1,-1,1}&)/@X],Method->"AdaptiveQuasiMonteCarlo"]]

2D example -- narrow peak

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

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

Timing[NIntegrate[(E^(-((x^2)+(y^2)))),{x,-5 π,10 π},{y,-5 π,10 π},Method->"MonteCarlo",MaxPoints->(10^6)]]

exact=(∫_(-5 π)^(10 π))(∫_(-5 π)^(10 π))(E^(-((x^2)+(y^2))))\[DifferentialD]y\[DifferentialD]x;
N[exact]

t=Table[NIntegrate[(E^(-((x^2)+(y^2)))),{x,-5 π,10 π},{y,-5 π,10 π},Method->"AdaptiveMonteCarlo"],{n,40}];
Mean[t]
Variance[t]

t=Table[NIntegrate[(E^(-((x^2)+(y^2)))),{x,-5 π,10 π},{y,-5 π,10 π},Method->"MonteCarlo"],{n,40}];Mean[t]
Variance[t]

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

3D example -- complicated boundary

Clear[x,y,z]

ineqn=Abs[x]+Sqrt[Abs[y]]+Abs[z]<4/5;

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

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

t=Reap[NIntegrate[Boole[ineqn],{x,-1,1},{y,-1,1},{z,-1,1},Method->{"AdaptiveMonteCarlo","BisectionDithering"->1/10,"Partitioning"->6,"MaxRecursion"->20,"SymbolicProcessing"->0},PrecisionGoal->2,EvaluationMonitor:>Sow[{x,y,z}]]];
t[[1]]
pnts=t[[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->{"AdaptiveMonteCarlo",RandomSeed->20,Method->{"MonteCarloRule","PointGenerator"->"Hammersley"}}]

NIntegrate[Exp[-((x^4)+(y^4))],{x,-2π+1/2,2π},{y,-2π+4/5,2π},Method->{"AdaptiveMonteCarlo",RandomSeed->20,Method->{"MonteCarloRule","PointGenerator"->"Sobol"}}]

NIntegrate[Exp[-((x^4)+(y^4))],{x,-2π+1/2,2π},{y,-2π+4/5,2π},Method->{"AdaptiveMonteCarlo",RandomSeed->20,Method->{"MonteCarloRule","PointGenerator"->"Niederreiter","Points"->100}}]

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.

Plot of the function

gr=Plot3D[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π, π}]

Even-odd preprocessing sampling points

grSP=NISamplingPoints[NIntegrate[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π,π},PrecisionGoal->2.5]][[1]];
grSP=Graphics[grSP,Axes->True,AxesOrigin->{- π,- π}]

No even-odd preprocessing sampling points

grSP=NISamplingPoints[NIntegrate[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π,π},Method->{SymbolicPiecewiseSubdivision,Method->GlobalAdaptive},PrecisionGoal->2.5]][[1]];grSP=Graphics[grSP,Axes->True,AxesOrigin->{- π,- π}]

Documentation

Integration strategies

Integration rules

Extenistion

Summary

Better incorporated symbolic preprocessing

Rules -- 1D and nD

Singularity handlers

Finite ranges oscillatory integration

Infinite ranges oscillatory integration

Adaptive Monte Carlo integration

Extensibility

NIntegrate and NDSolve similarities