NIntegrate[(1/(Sqrt[Abs[x-1]])),{x,0,2}]

points=Reap[NIntegrate[(1/(Sqrt[Abs[x-1]])),{x,0,2},EvaluationMonitor:>Sow[x]]][[2,1]];
Graphics[{PointSize[0.006],Point/@N@Transpose[{points,Range[Length[points]]}]},PlotRange->All,AspectRatio->1,Axes->True]

(RefLink[NIntegrate, paclet:ref/NIntegrate]) integration strategies.

(RefLink[NIntegrate, paclet:ref/NIntegrate]) preprocessor strategies.

|(∫_V)fdx-(Subscript[R, V])(f)|>|(∫_V)fdx- (Subscript[R, (Subscript[V, 1])])(f)+(Subscript[R, (Subscript[V, 2])])(f)|,

NIntegrate[1/Sqrt[x],{x,0,1},Method->"GlobalAdaptive"]

"GlobalAdaptive" options.

globalError≤globalIntegral (10^(-pg))∨globalError≤(10^(-ag)),

NIntegrate[Sin[x],{x,0,2π},MaxRecursion->100]

NIntegrate[1/Sqrt[x], {x, 0, 1},MaxRecursion->100,PrecisionGoal->17]//InputForm

NIntegrate[Sin[(x^2)+x], {x, 0,80Pi},Method->{"GlobalAdaptive",Method->"GaussKronrodRule"},MaxRecursion->100]

res=NIntegrate[Sin[(x^2)+x], {x, 0,80Pi},Method->{"GlobalAdaptive","MaxErrorIncreases"->10000,Method->"GaussKronrodRule"},MaxRecursion->20]

exact=Integrate[Sin[(x^2)+x], {x, 0,80Pi}];
Abs[res-exact]

{absc,weights,errweights}=NIntegrate`GaussKronrodRuleData[5,MachinePrecision];

IRuleEstimate[f_,{a_,b_}]:=
Module[{integral,error},
{integral,error}=(b-a)Total@MapThread[{f[#1]#2,f[#1]#3}&,{Rescale[absc,{0,1},{a,b}],weights,errweights}];
{integral,Abs[error]}
]

IStrategyGlobalAdaptive[f_,{aArg_,bArg_},tol_]:=
Module[{t,integral, error,regions, r1,r2,a=aArg,b=bArg,c},

{integral,error}=IRuleEstimate[f,{a,b}];
(* boundaries, integral, error *)
regions={{{a,b},integral,error}};

While[error>= tol*integral,
(* splitting of the region with the largest error *)
{a,b}=regions[[1,1]];c=((a+b)/2);

(* integration of the left region *)
{integral,error}=IRuleEstimate[f,{a,c}];
r1={{a,c},integral,error};

(* integration of the right region *)
{integral,error}=IRuleEstimate[f,{c,b}];
r2={{c,b},integral,error};

(* sort the regions: the largest error one is the first *)
regions=Join[{r1,r2},Rest[regions]];
regions=Sort[regions,#1[[3]]>#2[[3]]&];

(* global integral and error *)
{integral,error}=Total[Map[Rest[#1]&,regions]];

];

integral
];

f[x_]:=1/Sqrt[x]

res=IStrategyGlobalAdaptive[f,{0,1},(10^(-8))]

exact=Integrate[f[x],{x,0,1}]

Abs[res-exact]/exact

NIntegrate[1/Sqrt[x],{x,0,1},Method->"LocalAdaptive"]

"LocalAdaptive" options.

initialIntegral+regionError==initialIntegral

integralEst=Min[initialIntegral (10^(-pg))/eps,(10^(-ag))/eps];
integralEst+regionError==integralEst

Plot[(1/(10((((1/2)-x))^2)+1)),{x,0,1}]

sampledPoints=Reap[NIntegrate[(1/(10((((1/2)-x))^2)+1)),{x,0,1},Method->{"LocalAdaptive"},EvaluationMonitor:>Sow[x]]][[2,1]];
ListPlot[Transpose[{sampledPoints,Range[Length[sampledPoints]]}]]

sampledPoints=Reap[NIntegrate[(1/(10((((1/2)-x))^2)+1)),{x,0,1},Method->{"LocalAdaptive","Partitioning"->3},EvaluationMonitor:>Sow[x]]][[2,1]];
ListPlot[Transpose[{sampledPoints,Range[Length[sampledPoints]]}]]

sampledPoints=Reap[NIntegrate[(((x+y))^6),{x,0,1},{y,0,1},Method->{"LocalAdaptive","Partitioning"->{3,4}},EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
sampledPoints=Partition[sampledPoints,Length[sampledPoints]/(3*4+1)];
sampledPoints=Flatten[Rest[sampledPoints],1];
ListPlot[sampledPoints,AspectRatio->1,GridLines->{Range[3]/3,Range[4]/4}]

rulePoints=5;
sampledPoints=Reap[NIntegrate[(x^6),{x,0,1},Method->{"LocalAdaptive","SymbolicProcessing"->0,Method->{"ClenshawCurtisRule","Points"->rulePoints},"SingularityHandler"->None},EvaluationMonitor:>Sow[x]]][[2,1]];
Length[sampledPoints]
ListPlot[Transpose[{sampledPoints,Range[Length[sampledPoints]]}]]

totalRulePoints=2rulePoints-1;
totalPoints=(totalRulePoints-1)totalRulePoints+totalRulePoints;
((totalPoints-Length[sampledPoints])/totalPoints)//N

{absc,weights,errweights}=NIntegrate`ClenshawCurtisRuleData[6,MachinePrecision];

IRuleEstimate[f_,{a_,b_}]:=
Module[{integral,error,scaledAbsc},
scaledAbsc=Rescale[absc,{0,1},{a,b}];
{integral,error}=(b-a)Total@MapThread[{f[#1]#2,f[#1]#3}&,{scaledAbsc,weights,errweights}];
{integral,Abs[error],scaledAbsc}
]

LocalAdaptiveRecurrence[f_,{a_,b_},integralEst_]:=
Module[{regions,integral,error,scaledAbsc},

{integral,error,scaledAbsc}=IRuleEstimate[f,{a,b}];

If[N[integralEst+error]==N[integralEst],
(* Stopping criteria is satisfied *)
integral,
(* ELSE call itself recursively *)
regions = Partition[scaledAbsc,2,1];
Total[LocalAdaptiveRecurrence[f,#1,integralEst]&/@regions]
]
];

IStrategyLocalAdaptive[f_,{aArg_,bArg_},tol_]:=
Module[{integral, error,a=aArg,b=bArg,d=1,dummy},

If[a>b,{a,b}={b,a};d=-1];

(* initial integral estimate *)
{integral,error,dummy}=IRuleEstimate[f,{a,b}];

d*LocalAdaptiveRecurrence[f,{a,b},d*integral*tol/$MachineEpsilon]
];

f[x_]:=Sqrt[x]*Sin[x]

res=IStrategyLocalAdaptive[f,{0,8 π},(10^(-8))]

exact=Integrate[f[x],{x,0,8 π}]

Abs[res-exact]/exact

"GlobalAdaptive" Is Generally Better than "LocalAdaptive"

funcs={(Sqrt[x]),(1/(Sqrt[x])),((Sin[200 x])/(Sqrt[x])),Log[x],(x^26),(1/((10^4)((((1/2)-x))^2)+1)),-(((E^(1-(1/(1-x)))) Sin[1-(1/(1-x))])/(((1-x))^2))};
precs={6,8,10,12,14};
n=10; (* number of integrations to determine the timing *)
integrationRule=Automatic;

FRangesToCube[{{x,0,∞}}]

exactvals=Integrate[#,{x,0,1}]&/@funcs;

gatimings=Map[First@Timing[Do[NIntegrate[#[[1]],{x,0,1},PrecisionGoal->#[[2]],Method->{"GlobalAdaptive","SymbolicProcessing"->0,Method->integrationRule,"SingularityHandler"->None},MaxRecursion->200],{n}]]&,Outer[List,funcs,precs,1],{2}];

latimings=Map[First@Timing[Do[NIntegrate[#[[1]],{x,0,1},PrecisionGoal->#[[2]],Method->{"LocalAdaptive","SymbolicProcessing"->0,Method->integrationRule,"SingularityHandler"->None},MaxRecursion->200],{n}]]&,Outer[List,funcs,precs,1],{2}];

ganfevals=Map[(k=0;res=NIntegrate[#[[1,1]],{x,0,1},PrecisionGoal->#[[2]],Method->{"GlobalAdaptive","SymbolicProcessing"->0,Method->integrationRule,"SingularityHandler"->None},MaxRecursion->200,EvaluationMonitor:>k++];{k,Abs[res-#[[1,2]]]/Abs[#[[1,2]]]})&,Outer[List,Transpose[{funcs,exactvals}],precs,1],{2}];

lanfevals=Map[(k=0;res=NIntegrate[#[[1,1]],{x,0,1},PrecisionGoal->#[[2]],Method->{"LocalAdaptive","SymbolicProcessing"->0,Method->integrationRule,"SingularityHandler"->None},MaxRecursion->200,EvaluationMonitor:>k++];{k,Abs[res-#[[1,2]]]/Abs[#[[1,2]]]})&,Outer[List,Transpose[{funcs,exactvals}],precs,1],{2}];

Grid[Prepend[Transpose@{funcs,ColumnForm/@Table[precs,{Length[funcs]}],ColumnForm/@(latimings/gatimings),ColumnForm[First/@#]&/@ganfevals,ColumnForm[First/@#]&/@lanfevals},Style[#,"SmallText"]&/@{"functions","precision goals","\!\(\*FractionBox[\(LocalAdaptive\\ timings\), \(GlobalAdaptive\\ timings\)]\)","GlobalAdaptive\nfunction\nevaluations","LocalAdaptive\nfunction\nevaluations"}],Frame->All]

Grid[Prepend[Transpose@{funcs,ColumnForm/@Table[precs,{Length[funcs]}],ColumnForm[#[[2]]&/@#]&/@ganfevals,ColumnForm[#[[2]]&/@#]&/@lanfevals},Style[#,"SmallText"]&/@{"functions","precision goals","GlobalAdaptive\nrelative errors","LocalAdaptive\nrelative errors"}],Frame->All]

Point Singularities

NIntegrate[(1/(Sqrt[(x-(π/6))]))(1/(Sqrt[(x-(π/3))])),{x,0,(π/6),(π/3),(π/2)}]

NIntegrate[Log[(((1-x))^2)+(((1-y))^2)],{x,0,1,2},{y,0,1,2}]

NIntegrate[(1/(Sqrt[(x-(π/6))]))(1/(Sqrt[(x-(π/3))])),{x,0,(π/2)},Exclusions->{(π/6),(π/3)}]

NIntegrate[Log[(((1-x))^2)+(((1-y))^2)],{x,0,2},{y,0,2},Exclusions->{{1,1}}]

Curve, Surface, and Hypersurface Singularities

Plot3D[Log[(1-(x^2+y^2))^2],{x,0,2},{y,0,2},Exclusions->(x^2)+(y^2)==1]

NIntegrate[Log[(1-(x^2+y^2))^2],{x,0,2},{y,0,2},Exclusions->(x^2)+(y^2)==1]//Timing

NIntegrate[Log[(1-(x^2+y^2))^2],{x,0,2},{y,0,2}]//Timing

NIntegrate[Log[(1-(x^2+y^2))^2],{x,0,1},{y,0,(Sqrt[1-(x^2)]),1}]

Example Implementation of Curve, Surface, and Hypersurface Singularity Handling

Plot3D[Log[(1-(x^2+y^2))^2],{x,0,2},{y,0,2},Exclusions->(x^2)+(y^2)==1]

NIntegrate[Log[(1-(x^2+y^2))^2]*Boole[x^2+y^2<1],{x,0,2},{y,0,2}]+NIntegrate[Log[(1-(x^2+y^2))^2]*Boole[x^2+y^2>1],{x,0,2},{y,0,2}]//Timing

Plot3D[Log[(((1-(x+(((1-y))^2))))^2)],{x,-2,2},{y,-1,3},Exclusions-> x+(((1-y))^2)==1]

NIntegrate[Log[(1-(x+(1-y)^2))^2]*Boole[x+(1-y)^2<1],{x,-2,2},{y,-1,3},PrecisionGoal->4]+NIntegrate[Log[(1-(x+(1-y)^2))^2]*Boole[x+(1-y)^2>1],{x,-2,2},{y,-1,3},PrecisionGoal->4]//Timing

gr1={Red,Point/@N@Reap[NIntegrate[Log[(1-(x+(1-y)^2))^2]*Boole[x+(1-y)^2<1],{x,-2,2},{y,-1,3},PrecisionGoal->4,EvaluationMonitor:>Sow[{x,y}]]][[2,1]]};gr2={Blue,Point/@N@Reap[NIntegrate[Log[(1-(x+(1-y)^2))^2]*Boole[x+(1-y)^2>1],{x,-2,2},{y,-1,3},PrecisionGoal->4,EvaluationMonitor:>Sow[{x,y}]]][[2,1]]};
Graphics[{PointSize[0.006],gr1,gr2},Axes->True,AxesOrigin->{-2,-1}]

SingularManifoldNIntegrate[f_,ranges___,Equal[eq_,n_?NumericQ],opts___]:=
NIntegrate[f*Boole[eq<n],ranges,opts]+NIntegrate[f*Boole[eq>n],ranges,opts]

f[x_,y_,z_]:=Log[((1-(x+(1-y)^2+(1-z)^2)))^2];

SingularManifoldNIntegrate[f[x,y,z],{x,-2,2},{y,-1,3},{z,-1,1},x+(1-y)^2+(1-z)^2==1,PrecisionGoal->3]

gr1={Red,Point[Re[#]]&/@Reap[NIntegrate[f[x,y,z]*Boole[x+(1-y)^2+(1-z)^2<1],{x,-2,2},{y,-1,3},{z,-1,1},PrecisionGoal->3,EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]]};gr2={Blue,Point[Re[#]]&/@Reap[NIntegrate[f[x,y,z]*Boole[x+(1-y)^2+(1-z)^2>1],{x,-2,2},{y,-1,3},{z,-1,1},PrecisionGoal->3,EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]]};
Graphics3D[{PointSize[0.006],gr1,gr2},Axes->True]

"GlobalAdaptive" and "LocalAdaptive" singularity handling options.

NIntegrate[(1/(Sqrt[1-x])),{x,0,1},Method->{"GlobalAdaptive","SingularityHandler"->{IMT, "TuningParameters"->{10,2}}}]

IMT singularity handler options.

φ[a_,p_,t_]:=Exp[a(1-(1/(t^p)))];
φ[t_]:=φ[1,1,t]

Limit[D[φ[a,p,t],t],t->0,Assumptions->{a>0,p>0}]

Plot[φ[t],{t,0,1},AxesOrigin->{0,-0.02},PlotRange->All,AspectRatio->Automatic]

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

f[x_]:=1/Sqrt[x]

Plot[f[x],{x,0,1}]

{a,b}={0,1/16};
f[Rescale[φ[t],{0,1},{a,b}]]D[Rescale[φ[t],{0,1},{a,b}],t]

{a,b}={0,1/16};
Plot[f[Rescale[φ[t],{0,1},{a,b}]]D[Rescale[φ[t],{0,1},{a,b}],t]//Evaluate,{t,0,1},AxesOrigin->{0,-0.02},PlotRange->All]

{absc,weights,errweight}=NIntegrate`GaussKronrodRuleData[5,MachinePrecision];

abscGK=Rescale[#1,{0,1},{0,(1/16)}]&/@absc

derivativesGK=D[Rescale[t,{0,1},{0,(1/16)}],t]/.t->#&/@absc

(f[abscGK]derivativesGK).weights

abscGKIMT=Rescale[φ[#1],{0,1},{0,(1/16)}]&/@NIntegrate`GaussKronrodRuleData[5,MachinePrecision][[1]]

derivativesGKIMT=D[Rescale[φ[t],{0,1},{0,(1/16)}],t]/.t->#&/@NIntegrate`GaussKronrodRuleData[5,MachinePrecision][[1]]

(f[abscGKIMT]derivativesGKIMT).weights

Integrate[(1/(Sqrt[x])),{x,0,(1/16)}]

NIntegrate[(1/(Sqrt[1-x])),{x,0,1},Method->{"GlobalAdaptive","SingularityHandler"->"DoubleExponential"}]

Needs["Integration`NIntegrateUtilities`"];

exact=2;
tbl=((t={"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[(1/(Sqrt[x])),{x,0,1},Method->{"GlobalAdaptive","SingularityHandler"->#1[[1]],"SingularityDepth"->#1[[2]],"SymbolicProcessing"->0},MaxRecursion->100]];{#1[[2]],Abs[t[[1,1]]-exact],t[[2]],t[[3]]})&)/@{{"IMT",Infinity},{"IMT",1},{"DoubleExponential",1},{"IMT",4},{"DoubleExponential",4},{Automatic,4}};
TableForm[tbl,TableHeadings->Map[Style[#,FontFamily->Times,FontSize->11]&,{{"No singularity handling","IMT","DoubleExponential","IMT","DoubleEponential","Automatic"},{"SingularityDepth",ColumnForm[{"Difference from","the exact integral"}],ColumnForm[{"Number of function","evaluations"}],"Time (s)"}},{-1}]]

f[x_]:=(70/((10^4) (((x-(1/32)))^2)+(1/16)));
exact=Integrate[f[x],{x,0,1}];
tbl=((t={"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[f[x],{x,0,1},Method->{"GlobalAdaptive","SingularityHandler"->#1[[1]],"SingularityDepth"->#1[[2]],"SymbolicProcessing"->0},MaxRecursion->100,PrecisionGoal->8]];{#1[[2]],Abs[t[[1,1]]-exact],t[[2]],t[[3]]})&)/@{{"IMT",Infinity},{"IMT",1},{"DoubleExponential",1},{"IMT",4},{"DoubleExponential",4},{Automatic,4}};
TableForm[tbl,TableHeadings->Map[Style[#,FontFamily->Times,FontSize->11]&,{{"No singularity handling","IMT","DoubleExponential","IMT","DoubleEponential","Automatic"},{"SingularityDepth",ColumnForm[{"Difference","from the exact integral"}],ColumnForm[{"Number of function","evaluations"}],"Time (s)"}},{-1}]]

f[x_]:=((x+(1/(-1+Log[x])))/(x Log[x]));
exact=Integrate[f[x],{x,0,1}];
tbl=((t={"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[f[x],{x,0,1},Method->{"GlobalAdaptive","SingularityHandler"->#1[[1]],"SingularityDepth"->#1[[2]],"SymbolicProcessing"->0},MaxRecursion->3000,PrecisionGoal->6]];{#1[[2]],Abs[t[[1,1]]-exact],t[[2]],t[[3]]})&)/@{{"IMT",Infinity},{"IMT",1},{"DoubleExponential",1},{"IMT",4},{"DoubleExponential",4},{Automatic,4}};
TableForm[tbl,TableHeadings->Map[Style[#,FontFamily->Times,FontSize->11]&,{{"No singularity handling","IMT","DoubleExponential","IMT","DoubleEponential","Automatic"},{"SingularityDepth",ColumnForm[{"Difference from", "the exact integral"}],ColumnForm[{"Number of function","evaluations"}],"Time (s)"}},{-1}]]

NIntegrateProfile@NIntegrate[(1/(Sqrt[x]))+y,{x,0,1},{y,0,1}]

NIntegrateProfile@NIntegrate[(1/(Sqrt[x]))+y,{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityHandler"->None},MaxRecursion->30]

NIntegrateProfile@NIntegrate[(1/(Sqrt[(x^2)+(y^2)])),{x,0,1},{y,0,1}]

NIntegrateProfile@NIntegrate[(1/(Sqrt[(x^2)+(y^2)])),{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityHandler"->"IMT"},MaxRecursion->30]

NIntegrateProfile@NIntegrate[(1/(Sqrt[(x^2)+(y^2)])),{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityHandler"->None},MaxRecursion->30]

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

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

(∫_0^1)(∫_0^x)(1/(Sqrt[4(x^2)+(y^2)]))dxdy ⇔(∫_0^1)(∫_0^1)(x/(Sqrt[4(x^2)+(((x y))^2)]))dxdy ⇔(∫_0^1)(∫_0^1)(1/(Sqrt[(y^2)+4]))dxdy.

(∫_0^1)(∫_0^1)(1/(Sqrt[4(x^2)+(y^2)]))dxdy,

(∫_0^1)(∫_0^y)(1/(Sqrt[4(x^2)+(y^2)]))dydx ⇔(∫_0^1)(∫_0^x)(1/(Sqrt[4(y^2)+(x^2)]))dxdy ⇔(∫_0^1)(∫_0^1)(x/(Sqrt[4(((x y))^2)+(x^2)]))dxdy ⇔(∫_0^1)(∫_0^1)(1/(Sqrt[1+4(y^2)]))dxdy.

Reap[NIntegrate[(1/(Sqrt[(y^2)+4]))+(1/(Sqrt[1+4 (y^2)])),{x,0,1},{y,0,1},PrecisionGoal->8,EvaluationMonitor:>Sow[{x,y}]]][[2,1]]//Length

Reap[NIntegrate[(1/(Sqrt[(y^2)+4]))+(1/(Sqrt[1+4 (x^2)])),{x,0,1},{y,0,1},PrecisionGoal->8,EvaluationMonitor:>Sow[{x,y}]]][[2,1]]//Length

exact=Integrate[1/Sqrt[x^2 + y^2],{x,0,1},{y,0,1}];tbl=((t={"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[1/Sqrt[x^2+y^2],{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityHandler"->#1[[1]],"SingularityDepth"->#1[[2]],"SymbolicProcessing"->0},MaxRecursion->12]];{#1[[2]],Abs[t[[1,1]]-exact],t[[2]],t[[3]]})&)/@{{None,Infinity},{"IMT",1},{"IMT",4},{"DuffyCoordinates",4},{"DuffyCoordinates",1}};
TableForm[tbl,TableHeadings->Map[Style[#,FontFamily->Times,FontSize->11]&,{{"No singularity handling","IMT","IMT","DuffyCoordinates","DuffyCoordinates"},{"SingularityDepth",ColumnForm[{"Difference","from the","nexact integral"}],ColumnForm[{"Number of","function","evaluations"}],"Time (s)"}},{-1}]]

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

"DuffyCoordinates" options.

NIntegrate`DuffyCoordinatesIntegrand[(1/(Sqrt[(x^3)+(((1-y))^3)+(z^3)])),{x,0,1},{y,0,1},{z,0,1},Method->{"DuffyCoordinates","Corners"->{0,1,0}}]//
Simplify[#,Assumptions->{0≤x≤1,0≤y≤ 1,0≤z≤1}]&

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

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

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

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

NIntegrate[(1/(Sqrt[x+Sin[1-y]])),{x,0,1},{y,0,1}]//Timing

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

FRangesToCube[ranges_,cubeSides:{{_,_}...}]:=
Module[{t,t1,jac,vars,rules={}},
vars=First/@ranges;
t=MapThread[(t1=Rescale[#1[[1]],#2,{#1[[2]],#1[[3]]}/.rules];
AppendTo[rules,#1[[1]]->t1];t1)&,{ranges,cubeSides}];
jac=Times@@MapThread[D[#1,#2]&,{t,vars}];
{rules,jac}
]/;Length[ranges]==Length[cubeSides];
FRangesToCube[ranges_,cubeSide:{_,_}]:=FRangesToCube[ranges,Table[cubeSide,{Length[ranges]}]];
FRangesToCube[ranges_]:=FRangesToCube[ranges,{0,1}];

FRangesToCube[{{x,0,8},{y,x,∞}}]

DuffyCoordinatesAtOrigin[F_,ranges___]:=DuffyCoordinatesBounds[F,First/@{ranges},Transpose[Rest/@{ranges}]];
DuffyCoordinatesBounds[F_,vars_,bounds_]:=
Module[{rules,jac,newF,rots,res,range},
{rules,jac}=FRangesToCube[Transpose[Prepend[bounds,vars]]];
newF=(F/.rules)*jac;
rots=NestList[RotateLeft[#1]&,vars,Length[vars]-1];
res=Prepend[Map[newF/.Thread[vars->#1]&,Rest[rots]],newF];
range=Join[{{vars[[1]],0,1}},Map[{#,0,vars[[1]]}&,Rest[vars]]];
{rules,jac}=FRangesToCube[range];
{(Total[res]/.rules)*jac,Sequence @@( {#1, 0, 1} &  /@ vars)}
];

DuffyCoordinates[F_,ranges___]:=DuffyCoordinates[F,ranges,Table[0,{Length[{ranges}]}]];
DuffyCoordinates[F_,rangesSeq__,corner_?(VectorQ[#1,IntegerQ]&)]:=
Module[{factor,ranges={rangesSeq},newrange,t},
factor=1;
newrange={};
MapIndexed[(
t=ranges[[#2[[1]]]];
If[#1==0,
newrange=Append[newrange,t],
newrange=Append[newrange,{t[[1]],t[[3]],t[[2]]}];factor=-factor])&,
corner];
DuffyCoordinatesAtOrigin[factor*F,Sequence@@newrange]
];

DuffyCoordinates[F[x,y],{x,0,1},{y,0,1}]

DuffyCoordinates[F[x,y,z],{x,0,1},{y,0,1},{z,0,1}]

NIntegrate@@DuffyCoordinates[(1/(Sqrt[(x^2)+(y^2)+(z^2)])),{x,0,4},{y,0,3},{z,0,2}]

res=NINT@@DuffyCoordinates[(1/(Sqrt[(x^2)+(((3-y))^2)+(z^2)])),{x,0,4},{y,0,3},{z,0,2},{0,1,0}];
res=Hold[Evaluate[res]]/.NINT->NIntegrate;
Timing @@res

Timing @NIntegrate[(1/(Sqrt[(x^2)+(y^2)+(z^2)])),{x,0,4},{y,0,3},{z,0,2},Method->{"GlobalAdaptive","SingularityHandler"->None}]

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

Plot[Log[(((1-x))^2)],{x,0,2},PlotRange->All]

exact=(∫_0^2)Log[(((1-x))^2)]\[DifferentialD]x;
((exact-NIntegrate[Log[(((1-x))^2)],{x,0,2}])/exact)

NIntegrate[Log[(((1-x))^2)],{x,0,2},Method->"GlobalAdaptive",MaxRecursion->100]

InputForm/@Select[#,0.9<#<1.01&]&@Reap[NIntegrate[x,{x,0,2},EvaluationMonitor:>Sow[x]]][[2,1]]

InputForm/@Select[#,0.99999<#<1.00001&]&@Reap[
NIntegrate[Log[(((1-x))^2)],{x,0,2},EvaluationMonitor:>Sow[x]]][[2,1]]

{Reap[NIntegrate[Log[(((1-x))^2)],{x,0,1,2},EvaluationMonitor:>Sow[x]]][[2,1]]//Length,Timing[NIntegrate[Log[(((1-x))^2)],{x,0,1,2}]][[1]]}

{Reap[NIntegrate[Log[(((1-x))^2)],{x,0,2},MaxRecursion->20,EvaluationMonitor:>Sow[x]]][[2,1]]//Length,Timing[NIntegrate[Log[(((1-x))^2)],{x,0,2},MaxRecursion->20]][[1]]}

NIntegrate[(1/(Sqrt[Abs[BesselJ[2,x]]])),{x,1,20},MaxRecursion->1000]//InputForm

NIntegrate[(1/(Sqrt[Abs[BesselJ[2,x]]])),{x,1,5.135622301840683`,8.417244140399848`,11.619841172149059`,14.79595178235126`,17.959819494987826`,20},PrecisionGoal->8]//InputForm

Integrate[(1/x)Sin[x],{x,1,∞}]==Integrate[(1/x)Sin[(1/x)],{x,0,1}]

NIntegrate[(1/x)Sin[x],{x,1,∞}]

NIntegrate[(1/x)Sin[(1/x)],{x,0,1},Method->"GlobalAdaptive",MaxRecursion->100]

One-Dimensional Integration

((Subscript[E, m])/(Subscript[E, n]))≥ ((e^(-((c m)/(log m))))/(e^(-((c n)/(log n)))))⇔

((log (Subscript[E, m]))/(log (Subscript[E, n])))≤((-((c m)/(log m)))/(-((c n)/(log n))))==((m log n)/(n log m)).

k=0;
f[x_]:=(70/((10^4) (((x-(1/32)))^2)+(1/16)));
gr=Reap[NIntegrate[f[x],{x,0,1},EvaluationMonitor:>Sow[Point[{N[x],k++}]]]][[2,1]];Graphics[{PointSize[0.006],gr},AspectRatio->1,Axes->True,PlotRange->All,GridLines->{None,{97,160}}]

Multidimensional Integration

pointsAutomatic=Reap[NIntegrate[(1/(Sqrt[x+y])),{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityDepth"->2,"SingularityHandler"->Automatic},EvaluationMonitor:>Sow[{x,y}]]][[2,1]];pointsDuffy=Reap[NIntegrate[(1/(Sqrt[x+y])),{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityDepth"->2,"SingularityHandler"->"DuffyCoordinates"},EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
Row[{Graphics[{PointSize[0.015], Point /@ pointsAutomatic}, Axes -> True,ImageSize->200,PlotLabel->"Sampling\ Points:\ "<>ToString[Length[pointsAutomatic]]], 
   Graphics[{PointSize[0.015], Point /@ pointsDuffy}, Axes -> True,ImageSize->200,PlotLabel->"Sampling\ Points:\ "<>ToString[Length[pointsDuffy]]]}]

TableForm[{#,Timing[NIntegrate[(1/(Sqrt[x+y])),{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityHandler"->#}]][[1]]}&/@{Automatic,"DuffyCoordinates","IMT",None},TableHeadings->{None,{StyleForm[ColumnForm[{"Singularity", "handler"}],FontFamily->Times],StyleForm["Time (s)",FontFamily->Times]}}]

TableForm[{#,Timing[NIntegrate[(1/(Sqrt[y]))+(x^10),{x,0,1},{y,0,1},Method->{"GlobalAdaptive","SingularityHandler"->#},MaxRecursion->20]][[1]]}&/@{Automatic,"DuffyCoordinates","IMT",None},TableHeadings->{None,{StyleForm[ColumnForm[{"Singularity", "handler"}],FontFamily->Times],StyleForm["Time (s)",FontFamily->Times]}}]

NIntegrate[((Sqrt[x])/(x-2)),{x,0,2,5},Method->"PrincipalValue"]

"PrincipalValue" options.

NIntegrate[(1/(x+(x^2))),{x,-(1/2),0,1},Method->PrincipalValue]

NIntegrate[1/(x+x^2), {x, -2, -1, 0, 1},Method->PrincipalValue]

NIntegrate[1/(x+x^2),{x,-2,1},Method->PrincipalValue,Exclusions->{-1,0}]

% + 2Log[2]

Sampling Points Visualization

k=0;
tbl=Reap[NIntegrate[1/Log[x],{x,0,1,2},Method->{"PrincipalValue","SingularPointIntegrationRadius"->1/4},EvaluationMonitor:>Sow[{x,++k}]]][[2,1]];
ListPlot[tbl,PlotRange->All]

Clear[f];f[x_?NumericQ]:=(AppendTo[tbl,{x,++k}];1/Log[x]);

k=0;tbl={};
NIntegrate[f[x],{x,0,1,2},Method->{"PrincipalValue","SingularPointIntegrationRadius"->1/4}];
ListPlot[tbl,PlotRange->All]

NIntegrate[1/Sqrt[x],{x,0,1},Method->"DoubleExponential"]

"DoubleExponential" options.

(∫_a^b)f(t)dx

(∫_(-∞)^(+∞))f(ϕ(t)) ϕ'(t)dx,

(∫_a^b)f(x)dx⟹x==((a+b)/2)+(1/2) (b-a) tanh((1/2) π sinh(x)),

DE(h)==h (∑_(i=-∞)^∞)f(ϕ(i h)) ϕ'(i h)

NIntegrate[(1/(Sqrt[x+y])),{x,0,1},{y,0,1},Method->"DoubleExponential",MaxRecursion->200]

tbl=Reap[NIntegrate[Sqrt[x]Sqrt[y],{x,0,1},{y,0,1},Method->"DoubleExponential",MaxRecursion->200,EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
Graphics[{PointSize[0.005],Point/@N[tbl]},Axes->True]

NIntegrate[((Log[(1/x)])/(x^(1/4))),{x,0,1},PrecisionGoal->10,Method->{"DoubleExponential","SymbolicProcessing"->0}]//Timing

NIntegrate[((Log[(1/x)])/(x^(1/4))),{x,0,1},PrecisionGoal->10,MaxRecursion->100,Method->{"GlobalAdaptive","SingularityDepth"->∞,"SymbolicProcessing"->0}]//Timing

methods={"DoubleExponential","GlobalAdaptive"};
pgoals=Range[5,15,2];
TableForm[Outer[(k=0;res=NIntegrate[(1/(Sqrt[x])),{x,0,1},Method->#1,PrecisionGoal->#2,MaxRecursion->20,EvaluationMonitor:>k++];{Abs[res-2]/2,k})&,methods,pgoals]//Transpose,TableHeadings->{pgoals,methods},TableDepth->2]

methods={"DoubleExponential","GlobalAdaptive"};
pgoals=Range[6,10,2];
TableForm[Outer[(k=0;res=NIntegrate[Exp[20(x-1)]Sin[256 x],{x,0,1},Method->{#1,"SymbolicProcessing"->0},PrecisionGoal->#2,MaxRecursion->20,EvaluationMonitor:>k++];{Abs[res-2]/2,k})&,methods,pgoals,1]//Transpose,TableHeadings->{pgoals,methods},TableDepth->2]

k=0;
{NIntegrate[Abs[Sin[3*x]],{x,0,π},Method->{"DoubleExponential","SymbolicProcessing"->0},MaxRecursion->10,EvaluationMonitor:>k++],k}

k=0;{NIntegrate[Abs[Sin[3*x]],{x,0,π},Method->{"GlobalAdaptive","SymbolicProcessing"->0},MaxRecursion->10,EvaluationMonitor:>k++],k}

NIntegrate[(1/(16 (((x-(π/4)))^2)+(1/16))),{x,0,1},Method->"DoubleExponential"]//Timing

NIntegrate[(1/(16 (((x-(π/4)))^2)+(1/16))),{x,0,1},Method->"GlobalAdaptive"]//Timing

Block[{gr,gr1},
gr=Plot[(1/(16 (((x-(π/4)))^2)+(1/16))),{x,0,1},PlotRange->All];gr1=Plot[D[(1/(16 (((x-(π/4)))^2)+(1/16))),x]//Evaluate,{x,0,1},PlotRange->All];
GraphicsGrid[{{gr,gr1}}]
]

(e^(-((c n)/(log(n))))),

DERuleEstimate[f_,{a_,b_},h_,wprec_:MachinePrecision]:=
Block[{$MaxExtraPrecision=50000,ϕ,F,i,j,temp,s1,s2},
ϕ[t_]:=Rescale[(1/2)(Tanh[(1/2)π Sinh[t]]+1),{0,1},{a,b}];
F[t_]:=Evaluate[f[ϕ[t]]*D[ϕ[t],t]];
i=1;
s1=FixedPoint[(temp=F[i*h];i++;N[N[temp,3*wprec]+#1,wprec])&,0];
j=-1;
s2=FixedPoint[(temp=F[j*h];j--;N[N[temp,3*wprec]+#1,wprec])&,0];
{i-j+1,h(s1+F[0]+s2)}
];

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

exact=Integrate[f[x],{x,0,1}]

{a,b}={0,1};
wprec=30;
range=Table[(1/i),{i,2,7}];
range=Join[range,Mean/@Partition[range,2,1]];
range=Sort[range,Greater];
err=Map[DERuleEstimate[f,{a,b},#1,wprec]&,range];
err=Map[{#1[[1]],Abs[exact-#1[[2]]]}&,err]; (* errors *)
logErr=Map[{#1[[1]],Log[#1[[2]]]}&,err]; (* logarithm of the errors *)
points=First/@err;

p[x_]:=Evaluate[Fit[logErr,{1,(x/(Log[x]))},x]]

p[x]

ListLinePlot[{logErr,{#1,p[#1]}&/@points},PlotRange->All,PlotStyle->{{Red},{Blue}}]

IRuleEstimate[F_,h_,oldSum_:None]:=
Block[{$MaxExtraPrecision=50000,step,i,temp,s1,s2},
If[oldSum===None,step=1,step=2];
i=1;
s1=FixedPoint[(temp=F[i*h];i+=step;N[N[temp,60]+#1])&,0];
i=-1;
s2=FixedPoint[(temp=F[i*h];i-=step;N[N[temp,60]+#1])&,0];
If[oldSum===None,h(s1+F[0]+s2),h(s1+s2)+(oldSum/2)]
];

Options[IStrategyDoubleExp]={"MaxRecursion"->7};
IStrategyDoubleExp[f_,{a_,b_},tol_,opts___]:=
Module[{ϕ,F,h,t,temp,k=0,maxrec},
maxrec="MaxRecursion"/.{opts}/.Options[IStrategyDoubleExp];
ϕ[t_]:=Evaluate[Rescale[(1/2)(Tanh[(1/2)π Sinh[t]]+1),{0,1},{a,b}]];
F[t_]:=Evaluate[f[ϕ[t]]*D[ϕ[t],t]];
h=1;
NestWhile[((temp=IRuleEstimate[F,h/=2,#1])&&k++<maxrec)&,IRuleEstimate[F,h,None],(Abs[#1]*tol<=Abs[#1-#2])&,2];
temp
];

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

IStrategyDoubleExp[f,{0,1},(10^(-8))]//InputForm

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

f[x_]:=Cos[64*Sin[x]]

res =IStrategyDoubleExp[f,{0,π},(10^(-8))];res//InputForm

exact=Integrate[f[x],{x,0,π}]

exact//N//InputForm

Abs[res-exact]/exact

"Trapezoidal" options.

NIntegrate[(1/(Sqrt[x])),{x,0,1},Method->{"Trapezoidal","SymbolicProcessing"->0}]

k=1;
tab=Table[{t,({"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[(1/π)Cos[t Sin[x]-k x],{x,0,π},Method->"GaussKronrodRule"],1])[[2]],({"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[(1/π)Cos[t Sin[x]-k x],{x,0,π},Method->"Trapezoidal"],1])[[2]]},{t,8,80,4}];
TableForm[tab,TableHeadings->{None,{t,"GlobalAdaptive","Trapezoidal"}}]

Example Implementation

TrapStep[f_,{a_,b_},n_?IntegerQ]:=
Module[{h,absc, is},
h=((b-a)/(n-1));
absc=Table[i,{i,a,b,h}];
is=h*Total[MapAt[#/2&,f/@absc,{{1},{-1}}]];
{is,∞,n}
];

TrapStep[f_,{a_,b_},{oldEstimate_,oldError_,oldn_}]:=
Module[{n,h,absc,is},
n=2 oldn-1;
h=((b-a)/(n-1));
absc=Table[i,{i,a+h,b-h,2 h}];
is=h*Total[f/@absc]+oldEstimate/2;
{is,Abs[is-oldEstimate],n}
];

Options[TrapezoidalIntegration]={"MaxRecursion"->7};
TrapezoidalIntegration[f_,{a_,b_},tol_,opts___]:=
Block[{maxrec,k=0,temp},
maxrec="MaxRecursion"/.{opts}/.Options[TrapezoidalIntegration];
NestWhile[((temp=TrapStep[f,{a,b},#])&&k++<maxrec)&,TrapStep[f,{a,b},5],#[[2]]>tol&][[1]];
temp[[1]]
]

f[x_]:=(1/π)Cos[80 Sin[x]- x]

res=TrapezoidalIntegration[f,{0,π},(10^(-5))]//N

exact=Integrate[f[x],{x,0,π}]

Abs[res-exact]/exact

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,-∞,∞}]

(∫_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,

k=0;NIntegrate[((Sin[200 x])/(x^2)),{x,(2/10),2},EvaluationMonitor:>k++];k

NIntegrate[((Sin[200 x])/(x^2)),{x,(2/10),2}]//Timing

k=0;NIntegrate[((Sin[20000 x])/(x^2)),{x,(2/10),2},EvaluationMonitor:>k++];k

NIntegrate[((Sin[20000 x])/(x^2)),{x,(2/10),2}]//Timing

k=0;NIntegrate[((Sin[200 x])/(x^2)),{x,(2/10),2},Method->{Automatic, "SymbolicProcessing"->0},EvaluationMonitor:>k++];k

NIntegrate[((Sin[200 x])/(x^2)),{x,(2/10),2},Method->{Automatic, "SymbolicProcessing"->0}]//Timing

k=0;NIntegrate[((Sin[20000 x])/(x^2)),{x,(2/10),2},Method->{Automatic, "SymbolicProcessing"->0},EvaluationMonitor:>k++];k

NIntegrate[((Sin[20000 x])/(x^2)),{x,(2/10),2},Method->{Automatic, "SymbolicProcessing"->0}]//Timing

NIntegrate[Sin[200 x^2+5](1/((x+1)^2)),{x,0,∞},Method->"ExtrapolatingOscillatory"]

(∫_a^∞)k(x) f(x)dx,

(∫_a^(Subscript[z, 1]))k( x) f (x)dx,(∫_a^(Subscript[z, 2]))k(x) f(x)dx,…,(∫_a^(Subscript[z, i]))k(x) f(x)dx,…

(∫_a^(Subscript[z, 1]))k(x) f(x)dx,(∫_(Subscript[z, 1])^(Subscript[z, 2]))k(x) f(x)dx,…,(∫_(Subscript[z, i-1])^(Subscript[z, i]))k(x) f(x)dx,….

Example Implementation

Clear[ω,k,f];
ω=20;
k[x_]:=Sin[ω x];
f[x_]:=(1/((x+1)^2));

Plot[k[x]f[x],{x,0,10},PlotPoints->1000,PlotRange->All]

psum[i_?NumberQ]:=NIntegrate[k[x]f[x],{x,i(1/ω) π,(i+1)(1/ω) π}]

res=NSum[psum[i],{i,0,∞},Method->"AlternatingSigns","VerifyConvergence"->False]

exact=Integrate[k[x]f[x],{x,0,∞}]//N

((Abs[exact-res])/(Abs[exact]))

resEO=NIntegrate[Sin[ω x]f[x],{x,0,∞},Method->"ExtrapolatingOscillatory"]

((Abs[exact-resEO])/(Abs[exact]))

NIntegrate[Sin[2*x]*(1/(x^2)),{x,1,∞}, Method->{"DoubleExponentialOscillatory","SymbolicProcessing"->0}]

"DoubleExponentialOscillatory" options.

(Subscript[I, s])==(∫_0^∞)f (x) sin (ω x)dx.

x==((M ϕ(t))/ω), ϕ (t)==(t/(1-(e^(-2 t-β ((e^t)-1)-α (1-(e^(-t))))))),

α==β/(Sqrt[1+((M log(M+1))/(4 π))]),β==(1/4)

(Subscript[I, s])==(∫_0^∞)f (M (ϕ t)) sin (M (ϕ t))M (ϕ' t)/ωdt.

DEO((Subscript[I, s]),h,N)==M h (∑_(n=-(SubMinus[N]))^(SubPlus[N]))f (M (ϕ (n h)/ω) sin (M ϕ (n h)) ϕ' (n h)/ω, N==(SubMinus[N])+(SubPlus[N])+1.

(Subscript[ε, M])==((|(Subscript[I, s])-DEO((Subscript[I, s]),h,N)|)/(|(Subscript[I, s])|))≃(e^(-(A/h)))≃(e^(-((A M)/π))).

ε≥((σ(((|(Subscript[I, (Subscript[M, 2])])-(Subscript[I, (Subscript[M, 1])])|)/(|(Subscript[I, (Subscript[M, 2])])|))))^2)

ε≥((σ(((|(Subscript[I, (Subscript[M, 3])])-(Subscript[I, (Subscript[M, 2])])|)/(|(Subscript[I, (Subscript[M, 3])])|))))^((Subscript[M, 3])/(Subscript[M, 2])))

ε≥((σ(((|(Subscript[I, (Subscript[M, 4])])-(Subscript[I, (Subscript[M, 3])])|)/(|(Subscript[I, (Subscript[M, 4])])|))))^((Subscript[M, 4])/(Subscript[M, 3])))

(Subscript[I, c])==(∫_0^∞)f (x) cos (ω x)dx,

exact=Limit[Integrate[Exp[-c x]Log[x]Sin[x],{x,0,∞},Assumptions->c>0],c->0]

NIntegrate[Log[x]Sin[x],{x,0,∞}]-exact

exact=Integrate[Sin[20 x]Cos[18x] (1/(Sqrt[x+1])),{x,0,∞}]

NIntegrate[Sin[20 x]Cos[18x] (1/(Sqrt[x+1])),{x,0,∞},PrecisionGoal->10]-exact

Plot[{Sin[20 x]Cos[18x] (1/(Sqrt[x+1])),Sin[20x]},{x,0,3}]

NIntegrate[(E^(-((x^4)+(y^4)))),{x,-2,2},{y,-2,2},Method->"MonteCarlo"]

NIntegrate[(E^(-((x^4)+(y^4)))),{x,-2,2},{y,-2,2},Method->"QuasiMonteCarlo"]

"MonteCarlo" options.

"QuasiMonteCarlo" options.

(∫_V)f(x)dx==vol(V) ∫(1/(vol(V)))Boole (x∈V) f(x)dx==vol(V) E(f),

NIntegrate[(1/(Sqrt[x])),{x,0.01,1},Method->{"MonteCarlo","MaxPoints"->(10^3)}]

NIntegrate[(1/(Sqrt[x])),{x,0.01,1},Method->{"MonteCarlo","RandomSeed"->12}]//InputForm

NIntegrate[(1/(Sqrt[x])),{x,0.01,1},Method->{"MonteCarlo","RandomSeed"->12}]//InputForm

pnts=Reap[NIntegrate[(1/(Sqrt[x])),{x,0.01,1},Method->{"MonteCarlo","RandomSeed"->12},EvaluationMonitor:>Sow[x]]][[2,1]];
Take[pnts,20]

SeedRandom[12];RandomReal[{0.01,1},20]

E(f)==(1/(vol(V)))(∫_v)f(x)dx.

Var(f)==(1/4) ((((Subscript[E, 1])(f)-(Subscript[E, 2])(f)))^2)+(1/2) ((Subscript[Var, 1])(f)+(Subscript[Var, 2])(f))

NIntegrate[(x^2)+(y^2),{x,0,1/3,2/3,1},{y,0,1/3,2/4,3/4,1},Method->"MonteCarlo"]

NIntegrate[(x^2)+(y^2),{x,0, 1},{y,0,1},Method->{"MonteCarlo","Partitioning"->{3,4}}]

parts={3,4};
t=Reap[NIntegrate[1,{x,0,1},{y,0,1},Method->{"MonteCarlo","Partitioning"->parts,Method->{"MonteCarloRule","Points"->3}},EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
grX=(Line[{{#1,0},{#1,1}}]&)/@Table[i,{i,0,1,(1/(parts[[1]]))}];
grY=(Line[{{0,#1},{1,#1}}]&)/@Table[i,{i,0,1,(1/(parts[[2]]))}];
grLP=Point/@t;
Graphics[{grLP,grX,grY,Red,grLP},Axes->True]

ranges={{x,0,(1/4),(3/5),1},{y,0,(2/5),(3/4),1}};
t=Reap[NIntegrate[1,Evaluate[Sequence@@ranges],Method->{"MonteCarlo",Method->{"MonteCarloRule","Points"->3}},EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
grX=Line[{{#1,0},{#1,1}}]&/@Rest@ranges[[1]];
grY=Line[{{0,#1},{1,#1}}]&/@Rest@ranges[[2]];
grLP=Point/@t;
Graphics[{grLP,grX,grY,Red,grLP},Axes->True]

n=10;res=Timing[Do[NIntegrate[(((E^x)-1)/(E-1)),{x,0,1},Method->{"MonteCarlo","Partitioning"->#1,"MaxPoints"->(10^6)},PrecisionGoal->2],{n}]][[1]]/n&/@Range[4];
ColumnForm[res]

MonteCarloEstimate[f_,strata_,n_]:=({#1[[1]],(Sqrt[#1[[2]]])}&)[Total[({((Mean[f/@#1])/strata),((Variance[f/@#1])/(((strata^2) n)/strata))}&)/@Table[Random[Real,{((i-1)/strata),(i/strata)}],{i,strata},{(n/strata)}]]]

f[x_]:=(((E^x)-1)/(E-1))
N[(∫_0^1)f[x]\[DifferentialD]x]

t=Table[MonteCarloEstimate[f,i,1000],{i,1,40}];

ratios=Transpose[t][[2]]/Transpose[t][[2,1]];

p[x_]:=Evaluate@Fit[ratios,{(1/x)},x]
p[x]

ListLinePlot[{ratios,p/@Range[Length[ratios]]},PlotRange->All,PlotStyle->{{Red},{Blue}}]

NIntegrate[(E^(-((x^4)+(y^4)))),{x,-π,π},{y,-π,π},Method->"AdaptiveMonteCarlo"]

f[x_,y_]:=(E^(- 30((((x-5))^4)+(((y-5))^4))));
Plot3D[f[x,y],{x,0,10},{y,0,10},PlotPoints->30,PlotRange->All]

Integrate[f[x,y],{x,0,10},{y,0,10}]
%//N

Mean@Table[NIntegrate[f[x,y],{x,0,10},{y,0,10},Method->{"AdaptiveMonteCarlo","BisectionDithering"->0}],{20}]

t=Reap[NIntegrate[f[x,y],{x,0,10},{y,0,10},Method->{"AdaptiveMonteCarlo","BisectionDithering"->0,"RandomSeed"->10},PrecisionGoal->2,EvaluationMonitor:>Sow[{x,y,0}]]];
Print["Integral value ",t[[1]]]
cp=Plot3D[f[x,y],{x,0,10},{y,0,10},PlotPoints->30,PlotRange->All];
Graphics3D[{cp[[1]],PointSize[0.006],Point/@t[[2,1]]},BoxRatios->{1,1,0.4},PlotRange->All,Axes->True]

Mean@Table[NIntegrate[f[x,y],{x,0,10},{y,0,10},Method->{"AdaptiveMonteCarlo", "BisectionDithering"->(1/10)}],{30}]

t=Reap[NIntegrate[f[x,y],{x,0,10},{y,0,10},Method->{"AdaptiveMonteCarlo","BisectionDithering"->(1/10),RandomSeed->10},PrecisionGoal->2,EvaluationMonitor:>Sow[{x,y,0}]]];
Print["Integral value ",t[[1]]]
cp=Plot3D[f[x,y],{x,0,10},{y,0,10},PlotPoints->30,PlotRange->All];
Graphics3D[{cp[[1]],PointSize[0.006],Point/@t[[2,1]]},BoxRatios->{1,1,0.4},PlotRange->All,Axes->True]

f[x_,y_]:=(E^(-((((x+1))^2)+(((y+1))^2))))+(E^(-((((x-1))^2)+(((y-1))^2))))

Plot3D[f[x,y],{x,-π,π},{y,-π,π}]

{k=0;
(resMC=NIntegrate[f[x,y],{x,-π,π},{y,-π,π},Method->"MonteCarlo",PrecisionGoal->2,EvaluationMonitor:>k++])//Timing//First,k}

{k=0;
(resAMC=NIntegrate[f[x,y],{x,-π,π},{y,-π,π},Method->"AdaptiveMonteCarlo",PrecisionGoal->2,EvaluationMonitor:>k++])//Timing//First,k}

exact=Integrate[f[x,y],{x,-π,π},{y,-π,π}]

tblMC=Table[NIntegrate[f[x,y],{x,-π,π},{y,-π,π},Method->"MonteCarlo",PrecisionGoal->2],{100}];//Timing

{Abs[Mean[tblMC]-exact],Mean[Abs[((tblMC-exact)/exact)]],(((tblMC-exact).(tblMC-exact))/(Length[tblMC]))}

tblAMC=Table[NIntegrate[f[x,y],{x,-π,π},{y,-π,π},Method->"AdaptiveMonteCarlo",PrecisionGoal->2],{100}];//Timing

{Abs[Mean[tblAMC]-exact],Mean[Abs[((tblAMC-exact)/exact)]],(((tblAMC-exact).(tblAMC-exact))/(Length[tblAMC]))}

NIntegrate[(E^(-((x1^4)+(x2^4)+(x3^4)))),{x1,-π,π},{x2,-π,π},{x3,-π,π},Method->"MultiPeriodic"]

tbl=(First/@#&)/@Rest[NIntegrate`MultiPeriodicDump`copyrules];
tbl=MapIndexed[#1*2^(#2[[1]]+1)&,tbl];
mlen=Max[Length/@tbl];
tbl=Map[Join[#,Table["",{mlen-Length[#]}]]&,tbl];
Style[TableForm[Transpose[tbl],TableHeadings->{Automatic,Range[2,Length[tbl]+1]}],Small]

f[x1_,x2_,x3_,x4_,x5_,x6_,x7_,x8_]:=1/(((1+0.9671190054385935` x1+0.21216802639809276` x2+0.682779542171783` x3+0.32962509624641606` x4+0.5549215440908636` x5+0.7907543870000786` x6+0.8580353669569777` x7+0.4796298578498076` x8))^9)

tbl=Table[{"IntegralEstimate","Evaluations","Timing"}/.NIntegrateProfile[NIntegrate[f[x1,x2,x3,x4,x5,x6,x7,x8],{x1,0,1},{x2,0,1},{x3,0,1},{x4,0,1},{x5,0,1},{x6,0,1},{x7,0,1},{x8,0,1},Method->meth,MaxPoints->(10^8),PrecisionGoal->pg],1],{pg,1,4},{meth,{"MultiPeriodic",{"MultidimensionalRule", "Generators"->5},{"MultidimensionalRule", "Generators"->9}}}];
TableForm[Map[#[[3]]&,tbl,{2}],TableHeadings->Map[Style[#,FontFamily->Times,FontSize->11]&,{{"Precision goal -> 1","Precision goal -> 2","Precision goal -> 3","Precision goal -> 4"},{"MultiPeriodic",ColumnForm[{"MultidimensionalRule", "with 5 generators"}],ColumnForm[{"MultidimensionalRule", "with 9 generators"}]}},{-1}],TableSpacing->3]

TableForm[Map[#[[2]]&,tbl,{2}],TableHeadings->Map[Style[#,FontFamily->Times,FontSize->11]&,{{"Precision goal -> 1","Precision goal -> 2","Precision goal -> 3","Precision goal -> 4"},{"MultiPeriodic",ColumnForm[{"MultidimensionalRule", "with 5 generators"}],ColumnForm[{"MultidimensionalRule", "with 9 generators"}]}},{-1}],TableSpacing->3]

Options of "SymbolicPiecewiseSubdivision".

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,"sampling\norder"}]

gr=Plot3D[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π, π}];
grSP=Point/@Reap[NIntegrate[Boole[(x^2)+(y^2)>1] Sin[(x^2)+(y^2)],{x,- π, π},{y,-π,π},Method->{"SymbolicPiecewiseSubdivision",Method->"GlobalAdaptive"},PrecisionGoal->3,EvaluationMonitor:>Sow[{x,y}]]][[2,1]];
grSP=Graphics[{PointSize[0.005],grSP},Axes->True,AxesOrigin->{- π,- π}];
GraphicsGrid[{{gr,grSP}}]

"ExpandSpecialPiecewise"

f[x_,y_]:=Boole[(x^2)+2(y^2)<1]Abs[(x^2)+(y^3)-2]Abs[-(x^2)+(y^2)+1]Abs[(x^2)-3(y^2)+x](1/(Sqrt[(x^2)+(y^2)+10]));
NIntegrate[f[x,y],{x,-1,1},{y,-1,1},Method->{"SymbolicPiecewiseSubdivision","ExpandSpecialPiecewise"->{Boole},Method->"MonteCarlo"}]//Timing

NIntegrate[f[x,y],{x,-1,1},{y,-1,1},Method->{"SymbolicPiecewiseSubdivision",Method->"MonteCarlo"}]//Timing

Options of "EvenOddSubdivision".

gr=Plot3D[1+Sin[(Sqrt[(x^2)+(y^2)])],{x,-π,π},{y,-π,π}];
grSP=Point[Reap[NIntegrate[1+Sin[(Sqrt[(x^2)+(y^2)])],{x,-π,π},{y,-π,π},Method->{Automatic,"SymbolicProcessing"->0},EvaluationMonitor:>Sow[{x,y}]]][[2,1]]];
grSP=Graphics[{PointSize[0.01],grSP},Axes->True,AxesOrigin->{-π,-π}];
GraphicsGrid[{{gr,grSP}}]

Graphics[{PointSize[0.01],Point/@Reap[NIntegrate[1+Sin[(Sqrt[(x^2)+(y^2)])],{x,-π,π},{y,-π,π},EvaluationMonitor:>Sow[{x,y}]]][[2,1]]},Axes->True]

Transformation Theorem

grEven=Point/@Reap[NIntegrate[x,{x,0,1},{y,-x,x},{z,2,(y^2)},Method->{"SymbolicPreprocessing","UnitCubeRescaling"->False,Method->{"LobattoKronrodRule","GaussPoints"->5}},EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]];
gr=Point/@Reap[NIntegrate[x,{x,0,1},{y,-x,x},{z,2,(y^2)},Method->{"LobattoKronrodRule","GaussPoints"->5,"SymbolicProcessing"->0},EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]];
Graphics3D[{gr,Red,grEven},PlotRange->All,Axes->True,ViewPoint->{2.813, 0.765, 1.718}]

grEven=Point/@Reap[NIntegrate[x,{x,0,1},{y,-x,x},{z,2,y},Method->{"SymbolicPreprocessing","UnitCubeRescaling"->False,Method->{"LobattoKronrodRule","GaussPoints"->5}},EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]];
Graphics3D[{Red,grEven},PlotRange->All,Axes->True,ViewPoint->{2.813, 0.765, 1.718}]

"VerifyConvergence"

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

NIntegrate[x,{x,-∞,∞},Method->{"EvenOddSubdivision","VerifyConvergence"->False}]

"OscillatorySelection" method options.

NIntegrate[((Sin[20x])/(x+1)),{x,0,∞}]

NIntegrate[((Sin[20x])/(x+1)),{x,0,∞},Method->{Automatic,"SymbolicProcessing"->0}]

NIntegrate[(Piecewise[{{((BesselJ[3,-x])/(Sqrt[-x])), x<0}, {((Cos[200 x])/(Sqrt[x])), 0<x<20}, {((Sin[2 x+3])/(x^3)), x>30}, {(1/(Log[x])), True}}]),{x,-∞,∞}]

NIntegrate[((Exp[2I (x^2)])/(x+1)),{x,0,∞}]//InputForm

NIntegrate[((Cos[2 (x^2)])/(x+1)),{x,0,∞}]+I NIntegrate[((Sin[2 (x^2)])/(x+1)),{x,0,∞}]//InputForm

NIntegrate[(Piecewise[{{((Cos[200 x])/(x^6)), 0<x<20}, {(1/(Sqrt[x-20])), True}}]),{x,1,40},Method->{"SymbolicPiecewiseSubdivision",Method->{"OscillatorySelection",Method->"DoubleExponential","FourierFiniteRangeMethod"->{"LocalAdaptive","Partitioning"->3}}}]

k=0;pointsDELA=Reap[NIntegrate[(Piecewise[{{((Cos[200 x])/(x^6)), 0<x<20}, {(1/(Sqrt[x-20 ])), True}}]),{x,1,40},Method->{"SymbolicPiecewiseSubdivision",Method->{"OscillatorySelection",Method->"DoubleExponential","FourierFiniteRangeMethod"->{"LocalAdaptive","Partitioning"->3},"FourierInfiniteRangeMethod"->"ExtrapolatingOscillatory"}},EvaluationMonitor:>Sow[{x,k++}]]][[2,1]];
k=0;points=Reap[NIntegrate[(Piecewise[{{((Cos[200 x])/(x^6)), 0<x<20}, {(1/((Sqrt[x-20]) )), True}}]),{x,1,40},EvaluationMonitor:>Sow[{x,k++}]]][[2,1]];
grDELA = Graphics[{PointSize[0.01], Point /@ pointsDELA}, AspectRatio -> 1, Axes -> True, 
    PlotRange -> {{0, 40}, All}]; 
gr = Graphics[{PointSize[0.01], Point /@ points}, AspectRatio -> 1, Axes -> True, PlotRange -> {{0, 40}, All}]; 
GraphicsGrid[{{grDELA,gr}},ImageSize->Medium]

NIntegrate[(Piecewise[{{((BesselJ[3,-x])/(Sqrt[-x])), x<0}, {((Cos[200 x])/(Sqrt[x])), 0<x<20}, {((Sin[2 x+3])/(x^2)), x>30}, {(1/(Log[x])), True}}]),{x,-∞,∞},Method->{"SymbolicPiecewiseSubdivision",Method->{"OscillatorySelection","FourierInfiniteRangeMethod"->"ExtrapolatingOscillatory"}}]

NIntegrate[((Sin[2 (x^3)+3])/(x^2)),{x,1,∞},Method->"ExtrapolatingOscillatory"]//Timing

NIntegrate[((Sin[2 (x^3)+3])/(x^2)),{x,1,∞}]//Timing

Options of "UnitCubeRescaling".

NIntegrate[Sin[(x^2)+(y^2)]/(Sqrt[(x^2)+(y^2)]),{x,0,5},{y,0,(Sqrt[x])},Method->"UnitCubeRescaling"]//Timing

NIntegrate[Sin[(x^2)+(y^2)]/(Sqrt[(x^2)+(y^2)]),{x,0,5},{y,0,(Sqrt[x])},Method->{Automatic,"SymbolicProcessing"->0}]//Timing

(∫_(Subscript[a, 1])^(Subscript[b, 1]))(∫_((Subscript[a, 2])((Subscript[x, 1])))^((Subscript[b, 2])((Subscript[x, 1]))))…(∫_((Subscript[a, d])((Subscript[x, 1]),…,(Subscript[x, d-1])))^((Subscript[b, d])((Subscript[x, 1]),…,(Subscript[x, d-1]))))f((Subscript[x, 1]),…,(Subscript[x, d]))d(Subscript[x, 1])…d(Subscript[x, d])

(Subscript[x, i])==(Subscript[a, i])((Subscript[(OverHat[x]), 1]),…,(Subscript[(OverHat[x]), i-1]))+(Subscript[(OverHat[x]), i])((Subscript[b, i])( (Subscript[(OverHat[x]), 1]),…,(Subscript[(OverHat[x]), i-1]))-(Subscript[a, i])((Subscript[(OverHat[x]), 1]),…,(Subscript[(OverHat[x]), i-1]))), i==1,…,d.

J((Subscript[(OverHat[x]), 1]),…,(Subscript[(OverHat[x]), d]))==(∏_(i=1)^d)((Subscript[b, i])((Subscript[(OverHat[x]), 1]),…,(Subscript[(OverHat[x]), d]))-(Subscript[a, i])((Subscript[(OverHat[x]), 1]),…,(Subscript[(OverHat[x]), d]))).

NIntegrate[Exp[-1/10(x+y)](x^2),{x,0,∞},{y,0,∞},Method->{"UnitCubeRescaling","FunctionalRangesOnly"->False}]//Timing

NIntegrate[Exp[-1/10(x+y)](x^2),{x,0,∞},{y,0,∞},Method->{"UnitCubeRescaling","FunctionalRangesOnly"->True}]//Timing

Example Implementation

FRangesToCube[ranges_,cubeSides:{{_,_}...}]:=
Module[{t,t1,jac,vars,rules={}},
vars=First/@ranges;
t=MapThread[(t1=Rescale[#1[[1]],#2,{#1[[2]],#1[[3]]}/.rules];
AppendTo[rules,#1[[1]]->t1];t1)&,{ranges,cubeSides}];
jac=Times@@MapThread[D[#1,#2]&,{t,vars}];
{rules,jac}
]/;Length[ranges]==Length[cubeSides];
FRangesToCube[ranges_,cubeSide:{_,_}]:=FRangesToCube[ranges,Table[cubeSide,{Length[ranges]}]];
FRangesToCube[ranges_]:=FRangesToCube[ranges,{0,1}];

Rescale[x,{0,1},{a,b}]

{transRules,jacobian}=FRangesToCube[{{x,0,1},{y,a[x],b[x]}}];
transRules
jacobian

(F[x,y]/.transRules)jacobian

{transRules,jacobian}=FRangesToCube[{{x,0,∞},{y,a[x],b[x]}}];
transRules
jacobian

{transRules,jacobian}=FRangesToCube[{{x,(Subscript[a, 0]),(Subscript[b, 0])},{y,(Subscript[a, 1])[x],(Subscript[b, 1])[x]}}];
transRules
jacobian

"SymbolicPreprocessing" options.

grarr=Map[Graphics[{PointSize[0.01],Point/@Reap[NIntegrate[((Boole[x^2+y^2<1/2])/(Sqrt[(x^2)+(y^2)])),{x,-1,1},{y,0,(x^2)},Method->{"SymbolicPreprocessing","EvenOddSubdivision"->#[[1]],"SymbolicPiecewiseSubdivision"->#[[2]],Method->{"GlobalAdaptive",Method->"GaussKronrodRule","SingularityDepth"->∞}},PrecisionGoal->3,EvaluationMonitor:>Sow[{x,y}]]][[2,1]]},PlotRange->{{-1,1},{0,1}},Axes->True]&,Outer[List,{False,True},{False,True}],{-2}];
Grid[Join[{{"",SpanFromLeft,"SymbolicPiecewiseSubdivision",SpanFromLeft},{SpanFromAbove,SpanFromBoth,False,True}},{Join[{"EvenOddSubdivision",False},grarr[[1]]]},{Join[{SpanFromAbove,True},grarr[[2]]]}],Dividers->All]

ClosedContourIntegral[fexpr_,{x_,xpareq_},{y_,ypareq_},{θ_,0,2π},opts___]:=
NIntegrate[fexpr(Sqrt[(x^2)+(y^2)])/.{x->xpareq,y->ypareq},{θ,0,2π},Evaluate[Sequence@@Append[{opts},Method->"Trapezoidal"]]]

{a,b}={2,3};
exact=Integrate[(Sqrt[(a^2) ((Cos[θ])^2)+(b^(2 ))((Sin[θ])^2)]),{θ,0,2π}]

ep=ClosedContourIntegral[1,{x,a Cos[θ]},{y,b Sin[θ]},{θ,0,2 π}]

Abs[exact-ep]

FourierAnalysis[f_,{x_,xmin_,xmax_},nterms_,integrator_:(NIntegrate[##,Method->"GlobalAdaptive",MaxRecursion->30]&)]:=Block[{a,b,funcTerms},a=(2/((xmax-xmin)))Table[integrator[Cos[((2π )/(xmax-xmin))j x] f,{x,xmin,xmax}],{j,0,nterms}];
b=(2/((xmax-xmin)))Table[integrator[Sin[((2π )/(xmax-xmin))j x] f,{x,xmin,xmax}],{j,1,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
];

func=FourierAnalysis[(x^3)+(x^2),{x,-2,2},12,Integrate]

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

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

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