<<MaTeX`

How to integrate f (x)?

MaTeX["\\int_{x_{0}}^{x_{0}+\\Delta x} f(x) dx",FontSize->40]

f(x) is not known (we may just have some values (Subscript[y, n]) = f((Subscript[x, n])) or may just be too complicated to deal with it. Therefore we assume that is sufficiently smooth and ...

Methods

Taylor

... we Taylor expand it around the starting point of the range, (Subscript[x, 0]).

seriesLeft=Series[f[x],{x,x0,1}]

We can easily integrate any function in this way and we get a first result, that we call (Subscript[I, Left])

Ileft=Integrate[Normal[seriesLeft],{x,x0,x0+Δx}]

Note, however, that it is arbitrary around which point we expand in Taylor series. We may as well have expand it around the end-point of the integration range, (Subscript[x, 0])+Δx. In this case we define (Subscript[I, Right])

Series[f[x],{x,x0+Δx,1}]
Iright=Integrate[Normal[%],{x,x0,x0+Δx}]

The two integrals are, a priori, equally good approximations of the true integral and we can formulate a general prescription for an approximate value of our unknown integral      I≃ (Subscript[I, avg]) = (1/2) ( (Subscript[I, Left]) + (Subscript[I, Right]))

Iavg=Simplify[1/2(Ileft+Iright)]

Two f' terms are evaluated at close-by points, so we can further Taylor expand that term, e.g. expanding f' around (Subscript[x, 0])

Series[f'[x0+Δx],{Δx,0,1}]
f1=Normal[%]

Simplify[Iavg]//.f'[x0+Δx]->f1
FullSimplify[%]

The result is now depending on the the second derivative of f .

Taylor expand before Integration

Series[f[x],{x,x0+Δx,1}]
%/.f'[x0+Δx]->f1
IrightTilde=Integrate[Normal[%],{x,x0,x0+Δx}]

FullSimplify[1/2(Ileft+IrightTilde)]

Poly

GenericPoints[P_]:=Table[{x[k],y[k]},{k,1,P}]

GenericPoints[3]//TableForm

list={{a,b},{c,d},{e,f}}
Flatten[Table[{x[i]->list[[All,1]][[i]],y[i]->list[[All,2]][[i]]},{i,1,Length[list]}]]

Pol3[x_]:=Simplify[InterpolatingPolynomial[GenericPoints[3],x]]
PolN[x_,N_]:=Simplify[InterpolatingPolynomial[GenericPoints[N],x]]

PolNpoints[deg_,listofpoints_,z_]:=Module[{repl,len=Length[listofpoints]},
If[deg==len,
repl=Flatten[Table[{x[i]->listofpoints[[All,1]][[i]],y[i]->listofpoints[[All,2]][[i]]},{i,1,len}]];
(*Print[repl];*)
PolN[z,deg]//.repl,
Print["Degree of polynomial requested, deg=",deg," is too high/low for the number of points, N=",len,", to constrain it."];
PolNpoints[len,listofpoints,z]
]
]

Pol3points[listofpoints_]:=Module[{x1,y1,x2,y2,x3,y3},
{x1,y1}=listofpoints[[1]];
{x2,y2}=listofpoints[[2]];
{x3,y3}=listofpoints[[3]];
(*Print[x1,y1,x2,y2,x3,y3];*)
Pol3[z]//.{x[1]->x1,y[1]->y1,x[2]->x2,y[2]->y2,x[3]->x3,y[3]->y3}
]

Simplify[PolNpoints[3,{{0,-2},{1,1},{-1,1}},z]]
Plot[%,{z,-2,2}]

Simplify[Pol3points[{{0,-2},{1,1},{-1,1}}]]
Plot[%,{z,-2,2}]

Random Points

?RandomReal

RandomReal[{0,1},10]

Test each strategy against a function of known integral

Exact Integral

Let us study a function of known integral

g[x_]:=Exp[-x]Sin[x]

g[x_]:=Sqrt[x+1/Pi^2]Sin[(x+1/Pi^2)^(-1)]

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

Integrate[g[x],x]

We can easily obtain finite integrals in closed form, either symbolically or numerically

Integrate[g[x],{x,0,1}]

exctIntegral=NIntegrate[g[x],{x,0,1}]

Evaluation at fixed (Subscript[x, n]) points - 2-point  approximation

We can integrate it using a variable number of points according to the formula below

A quick code would return

nP=1000;
listOfPoints=Table[x,{x,0,1,1/nP}];

listOfPoints;

N[Total[Table[g[x],{x,listOfPoints}]]]

1/nP*N[Total[Table[g[x],{x,listOfPoints}]]](* 10 punti *)

1/nP*N[Total[Table[g[x],{x,listOfPoints}]]] (* 1000 punti *)

exctIntegral

1/nP*N[Total[Map[g,listOfPoints]]]

Let us define a function that does the evaluation of this integral for a number of segments that is a parameter of the function, so that I can ask for the integral for 10, (10^2), (10^4) points just specifying how many points I want.

For convenience we define the function as the difference with respect to the (in this case) known integral

Approximate[nP_]:=Module[{$listOfPoints,$approximate},
$listOfPoints=Table[x,{x,0,1,1/nP}];
$approximate=1/nP*N[Total[Table[g[x],{x,$listOfPoints}]]];
$approximate
]

ExactMinusApproximate[nP_]:=Module[{$listOfPoints,$approximate,$exctIntegral},
$listOfPoints=Table[x,{x,0,1,1/nP}];
$approximate=1/nP*N[Total[Table[g[x],{x,$listOfPoints}]]];
$exctIntegral=NIntegrate[g[x],{x,0,1}];
$exctIntegral-$approximate
]

ExactMinusApproximate[10]

ExactMinusApproximate[1000]

We tabulate the difference between the exact and approximate integral to monitor the convergence of the numerical integration.

DeviationsVsNPoints=Table[{nPoints,Abs[ExactMinusApproximate[nPoints]]},{nPoints,{10,100,1000,3000,10000,30000}}]
ListLogLogPlot[%,Frame->True,FrameLabel -> {"number of points","Exact - Approximate Integral"}]

ApproximateIntegrals=Table[{nPoints,Approximate[nPoints]},{nPoints,{10,100,1000,3000,10000,30000}}]

ApproximateIntegrals[[-1,2]]

DeviationsVsNPointsLast=Map[({#[[1]],#[[2]]-ApproximateIntegrals[[-1,2]]})&,ApproximateIntegrals[[1;;-2]]]
ListLogLogPlot[DeviationsVsNPointsLast,Frame->True,FrameLabel -> {"number of points","Last - Approximate Integral"}]

It seems to follow a power-law (it is a straight line in a log-log plot).

Evaluation at fixed (Subscript[x, n]) points - N-point  approximation

Prototyping

gXn=Table[{x,g[x]},{x,listOfPoints}]

Length[Drop[gXn,0]]

Forming several consecutive intervals each made of three points

rank=3
Flatten[Table[Partition[Drop[gXn,rank-k],rank],{k,1,rank}],1]
Unset[rank]

Sort[Out[283]]-Sort[$123]

Out[268][[-1]]-one

Merge these intervals in one set

one=Partition[gXn,3]
two=Partition[Drop[gXn,1],3]
three=Partition[Drop[gXn,2],3]
$123=Join[one,two,three]

Pol3points[one[[1]]]

$3pols=Table[{$123[[k]],Pol3points[$123[[k]]]},{k,1,Length[$123]}];

?MinMax

range$pols=Map[({MinMax[#[[1]][[All,1]]],#[[2]]})&,$3pols]

range$pols[[1]]

?ColorData

ColorData[1,"ColorList"]

ColorData[3,"ColorList"]

ColorData[13,"ColorList"]

Table[Plot[UnitStep[z-Min[range$pols[[k,1]]]]UnitStep[-z+Max[range$pols[[k,1]]]]range$pols[[k,2]],{z,0,1},PlotStyle->{Dotted,ColorData[3,"ColorList"][[k]],Thick},PlotRange->{range$pols[[k,1]],All}],{k,1,Length[range$pols]}]
Plot[g[z],{z,0,1},PlotStyle->{Dashed,Red,Thick}]
Show[Join[{%,%%}]]

Table[Plot[Piecewise[{{range$pols[[k,2]],z>=Min[range$pols[[k,1]]]&&z<=Max[range$pols[[k,1]]]},{Null,Piecewise[{{range$pols[[k,2]],z<Min[range$pols[[k,1]]]||z>Max[range$pols[[k,1]]]}}]}}],{z,0,1},PlotStyle->{Dotted,ColorData[3,"ColorList"][[k]],Thick},PlotRange->{range$pols[[k,1]],All}],{k,1,Length[range$pols]}]
Plot[g[z],{z,0,1},PlotStyle->{Dashed,Red,Thick}]
Show[Join[{%,%%}]]

Function

range$pols[[1,1]]

Clear[FunctionToPiecewise]

Options[FunctionToPiecewise]={Normalization->True}
FunctionToPiecewise[g_,nPext_,rank_,OptionsPattern[]]:=Module[{gXn,$listOfPoints,$123,$3pols,range$pols,nP=nPext},
$listOfPoints=Table[x,{x,0,1,1/nP}];
Print[Length[$listOfPoints]," points"];
gXn=Table[{x,g[x]},{x,$listOfPoints}];
If[OptionValue[Normalization] ==True
,Print["Normalization will be preserved"];$123=Partition[gXn,UpTo[rank],rank-1];
,Print["Normalization will not be preserved"];$123=Flatten[Table[Partition[Drop[gXn,rank-k],rank],{k,1,rank}],1];
];
$3pols=Table[{$123[[k]],PolNpoints[rank,$123[[k]],z]},{k,1,Length[$123]}];
range$pols=Map[({MinMax[#[[1]][[All,1]]],#[[2]]})&,$3pols];
Table[Piecewise[{{range$pols[[k,2]],z>=Min[range$pols[[k,1]]]&&z<=Max[range$pols[[k,1]]]}},Null]
,{k,1,Length[range$pols]}]
]

Partition[{a,b,c,d,e,f,g},3,2]

aprroxf=FunctionToPiecewise[g,22,3];

With the normalization option the interval are not overlapping

Plot[aprroxf,{z,0,1}]
Plot[g[z],{z,0,1},PlotStyle->{Opacity[0.13,Cyan],Thickness[0.03]}];
Show[{%,%%}]

ExactMinusApproximateInterpolation[nP_,rank_]:=Module[{aprroxf,subIntegrals},
aprroxf=FunctionToPiecewise[g,nP,rank];
subIntegrals=Table[Integrate[aprroxf[[k]],{z,0,1}]/.Null->0,{k,1,Length[aprroxf]}];
Total[subIntegrals]-exctIntegral
]

DeviationsVsNPointsInterpolation=Table[{nPoints,Abs[ExactMinusApproximateInterpolation[nPoints,2]]},{nPoints,{10,100,1000(*,3000,10000,30000*)}}]

ListLogLogPlot[DeviationsVsNPointsInterpolation,FrameLabel -> {"number of points","Exact - Approximate Integral"}]

DeviationsVsNPointsInterpolationCubic=Table[{nPoints,Abs[ExactMinusApproximateInterpolation[nPoints,3]]},{nPoints,{10,100,1000(*,3000,10000,30000*)}}]

ListLogLogPlot[DeviationsVsNPointsInterpolationCubic,FrameLabel -> {"number of points","Exact - Approximate Integral"}]

DeviationsVsNPointsInterpolationQuartic=Table[{nPoints,Abs[ExactMinusApproximateInterpolation[nPoints,4]]},{nPoints,{10,100,1000(*,3000,10000,30000*)}}]

ListLogLogPlot[DeviationsVsNPointsInterpolationQuartic,FrameLabel -> {"number of points","Exact - Approximate Integral"}]

DeviationsVsNPointsInterpolationQuintic=Table[{nPoints,Abs[ExactMinusApproximateInterpolation[nPoints,5]]},{nPoints,{(*80,120,300*)119,211,327}}]

ListLogLogPlot[DeviationsVsNPointsInterpolationQuintic,FrameLabel -> {"number of points","Exact - Approximate Integral"}]

Evaluation at random (Subscript[x, n]) points - 1-point approximation

Prototyping

MaTeX["I = \\int_0^1 f(x)dx",FontSize->40]

MaTeX["N_p \\textrm{ random values } x_i \in [0,1]",FontSize->40]

MaTeX["\\bar{x}= x_i ,  ||\\bar{x}||=N_p",FontSize->40]

MaTeX["E_I(\\bar{x})= \\frac{1}{|| \\bar{x} ||}  \\sum_{x_i \in  \\bar{x}} f(x_i)",FontSize->40]

nP=10;
listOfPoints=RandomReal[{0,1},nP]
approximate=1/nP*N[Total[Table[g[x],{x,listOfPoints}]]]
Unset[nP]

MaTeX["\\langle E_I(\\bar{x}) \\rangle_{\\bar{x}} = \\frac{1}{ N_{replicas}} \\left( E_I(\\bar{x}_1) + E_I(\\bar{x}_2) + ... + E_I(\\bar{x}_{N_{replicas}})  \\right) ",FontSize->40]

MaTeX["\\sigma_{E_I}  = Std\\left( E_I(\\bar{x}_1) + E_I(\\bar{x}_2) + ... + E_I(\\bar{x}_{N_{replicas}})  \\right) ",FontSize->40]

MaTeX["I_{MC}= \\langle E_I(\\bar{x}) \\rangle_{\\bar{x}} \\pm \\sigma_{E_I}",FontSize->40]

??NIntegrate

Options[NIntegrate]

ExactMinusApproximateRandom[nP_]:=Module[{$listOfPoints,$approximate,$exctIntegral},
$listOfPoints=RandomReal[{0,1},nP];
(*$approximate=1/nP*N[Total[Table[g[x],{x,$listOfPoints}]]];*)
$approximate=1/nP*(Map[g,$listOfPoints]//Total//N);
$exctIntegral=NIntegrate[g[x],{x,0,1}];
$exctIntegral-$approximate
]

(* 100 Replicas *)
Table[ExactMinusApproximateRandom[1000],{k,100}];  
coarser=Histogram[%,{0.001},ChartStyle->Opacity[0.2,Blue]]

Np=10000;Nr=100;Color$=Orange;
Table[ExactMinusApproximateRandom[10000],{k,100}];
finer=Histogram[%,{0.001},ChartStyle->Opacity[0.2,Color$],ChartLegends->Placed[{MaTeX["N_{points}="<>ToString[Np]<>", N_{replicas}="<>ToString[Nr]]},Right]]

Options[DispersionPlot]={Np->10000,Nr->100,Color->Orange,BinSize->0.001}
DispersionPlot[opt:OptionsPattern[]]:=Module[{$np=OptionValue[Np],$nr=OptionValue[Nr],$Color=OptionValue[Color],$binsize=OptionValue[BinSize],$diff},
$diff=Table[ExactMinusApproximateRandom[$np],{k,$nr}];
Histogram[$diff,{$binsize},ChartStyle->Opacity[0.2,$Color],ChartLegends->Placed[{MaTeX["N_{points}="<>ToString[$np]<>", N_{replicas}="<>ToString[$nr]]},Right]]
]

Show[DispersionPlot[Np->1000,Nr->100,Color->Blue],DispersionPlot[Np->10000,Nr->100,Color->Orange],DispersionPlot[Np->100000,Nr->100,Color->Green]]

Show[coarser,finer]

devs=Table[ExactMinusApproximateRandom[10],{k,100}]

Mean[devs]
Variance[devs]
StandardDeviation[devs]
%^2

Function

Map[ExactMinusApproximateRandom[11]&,Range[100]]

AbsoluteTiming[ParallelMap[{#,StandardDeviation[Table[ExactMinusApproximateRandom[#],{k,100}]]}&,{10,100,1000,3000,10000,30000}]]

AbsoluteTiming[Table[{nPoints,StandardDeviation[ParallelMap[ExactMinusApproximateRandom[nPoints]&,Range[100]]]},{nPoints,{10,100,1000,3000,10000,30000}}]]

AbsoluteTiming[Table[{nPoints,StandardDeviation[Table[ExactMinusApproximateRandom[nPoints],{k,100}]]},{nPoints,{10,100,1000,3000,10000,30000}}]]

AbsoluteTiming[Table[{nPoints,StandardDeviation[ParallelTable[ExactMinusApproximateRandom[nPoints],{k,100}]]},{nPoints,{10,100,1000,3000,10000,30000}}]]

AbsoluteTiming[ParallelTable[{nPoints,StandardDeviation[ParallelTable[ExactMinusApproximateRandom[nPoints],{k,100}]]},{nPoints,{10,100,1000,3000,10000,30000}}]]

DeviationsVsNPointsRandom=Table[{nPoints,StandardDeviation[Table[ExactMinusApproximateRandom[nPoints],{k,100}]]},{nPoints,{10,100,1000,3000,10000,30000}}]

ListLogLogPlot[DeviationsVsNPointsRandom]

Power-law fits

??NonlinearModelFit

Taylor N=2

fit=NonlinearModelFit[DeviationsVsNPoints,{const+d*1/x^r,r>0},{{const,0},d,{r,2}},x]

We want to monitor how well the fit reproduces each point of our table. We fine a function that makes the ratio between the value of the residual at each point, (Subscript[r, i]) = | fit((Subscript[x, i])) - (Subscript[y, i]) |, and the value of the table, therefore it is a fit of “quality”   (Subscript[q, i])=(Subscript[r, i])/ (Subscript[y, i])

?AppendTo

DataWithResidual[myfit_]:=Module[{$rows},
$rows=Transpose[myfit["Data"]];
Transpose[AppendTo[$rows,myfit["FitResiduals"]]]
]
QualityOfFit[fit_]:=Map[#[[3]]/#[[2]]&,DataWithResidual[fit]]

Check-out the fit results ...

fit["BestFitParameters"]

fit["ParameterErrors"]

Show[LogLogPlot[fit["Function"][x],{x,10,10^4},PlotStyle->Red],ListLogLogPlot[fit["Data"]],Frame->True, FrameLabel -> {"number of points","Exact - Approximate Integral"}]

The fit is terrible! The reason is that the points for larger number of points the fit is failing to match them.

QualityOfFit[fit]

We can look at the documentation to see that fits are done giving “weights” to each point in the data. So we can manipulate these weights to give more importance to the points where the number of points in the approximate integral is large.

betterfit=NonlinearModelFit[DeviationsVsNPoints,{const+d*1/x^r,r>0},{{const,0},d,{r,2}},x,Weights->{10,100,1000,3000,10000,30000}]
qfit=QualityOfFit[betterfit]

Show[LogLogPlot[betterfit["Function"][x],{x,10,3*10^4},PlotStyle->Red],ListLogLogPlot[betterfit["Data"]]]

Show[LogLogPlot[betterfit["Function"][x],{x,10,3*10^4},PlotStyle->Red],ListLogLogPlot[betterfit["Data"]]]

DeviationsVsNPointsLogLog=Map[{Log10[#[[1]]],Log10[#[[2]]]}&,DeviationsVsNPoints]

fit=NonlinearModelFit[Drop[DeviationsVsNPointsLogLog,{1,2}],{const-x*r,r>0},{{const,0},{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]

Show[Plot[fit["Function"][x],{x,1,5},PlotStyle->Red,PlotRange->All],ListPlot[fit["Data"]],PlotRange->All]

Histogram[qfit,{0.02}]

Newton N=2

fit=NonlinearModelFit[DeviationsVsNPointsInterpolation,{const+d*1/x^r,r>0},{{const,0},d,{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]

Show[LogLogPlot[fit["Function"][x],{x,10,10^4},PlotStyle->Red,PlotRange->All],ListLogLogPlot[fit["Data"]],PlotRange->All]

DeviationsVsNPointsInterpolationLogLog=Map[{Log10[#[[1]]],Log10[#[[2]]]}&,DeviationsVsNPointsInterpolation]

fit=NonlinearModelFit[DeviationsVsNPointsInterpolationLogLog,{const-x*r,r>0},{{const,0},{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]

Show[Plot[fit["Function"][x],{x,1,3},PlotStyle->Red,PlotRange->All],ListPlot[fit["Data"]],PlotRange->All]

Newton N>2

DeviationsVsNPointsInterpolationCubicLogLog=Map[{Log10[#[[1]]],Log10[#[[2]]]}&,DeviationsVsNPointsInterpolationCubic]

fit=NonlinearModelFit[DeviationsVsNPointsInterpolationCubicLogLog,{const-x*r,r>0},{{const,0},{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]

Show[Plot[fit["Function"][x],{x,1,3},PlotStyle->Red,PlotRange->All],ListPlot[fit["Data"]],PlotRange->All]

DeviationsVsNPointsInterpolationQuarticLogLog=Map[{Log10[#[[1]]],Log10[#[[2]]]}&,DeviationsVsNPointsInterpolationQuartic]

fit=NonlinearModelFit[DeviationsVsNPointsInterpolationQuarticLogLog,{const-x*r,r>0},{{const,0},{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]

Show[Plot[fit["Function"][x],{x,1,3},PlotStyle->Red,PlotRange->All],ListPlot[fit["Data"]],PlotRange->All]

DeviationsVsNPointsInterpolationQuinticLogLog=Map[{Log10[#[[1]]],Log10[#[[2]]]}&,DeviationsVsNPointsInterpolationQuintic]

fit=NonlinearModelFit[DeviationsVsNPointsInterpolationQuinticLogLog,{const-x*r,r>0},{{const,0},{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]

Show[Plot[fit["Function"][x],{x,1,3},PlotStyle->Red,PlotRange->All],ListPlot[fit["Data"]],PlotRange->All]

Monte Carlo

DeviationsVsNPointsRandomLogLog=Map[{Log10[#[[1]]],Log10[#[[2]]]}&,DeviationsVsNPointsRandom]

fit=NonlinearModelFit[DeviationsVsNPointsRandomLogLog,{const-x*r,r>0},{{const,0},{r,2}},x]
fit["BestFitParameters"]
QualityOfFit[fit]
Histogram[%,{0.01}]

Show[Plot[fit["Function"][x],{x,1,5},PlotStyle->Red,PlotRange->All],ListPlot[fit["Data"]],PlotRange->All]

2D

g[x_,y_]:=Sin[y]*Sqrt[x+1/Pi^2]Sin[(x+1/Pi^2)^(-1)]

ContourPlot[g[x,y],{x,0,1},{y,0,1},PlotLegends->Automatic]

NIntegrate[g[x,y],{x,0,1},{y,0,1}]

exact2D=Integrate[g[x,y],{x,0,1},{y,0,1}]

N[exact2D,12]

3D

g[x_,y_,z_]:=Cos[z/(1/Pi+y)]*Sin[y]*Sqrt[x+1/Pi^2]Sin[(x+1/Pi^2)^(-1)]

NIntegrate[g[x,y,z],{x,0,1},{y,0,1},{z,0,1}]

ContourPlot[g[x,y,0.3],{x,0,1},{y,0,1},PlotLegends->Automatic]