Consider some dataset.
I use the following code to interpolate it:
SetDirectory[NotebookDirectory[]];
dataset1 = Drop[Import["dataset.txt", "Table"], 1];
yxDistr[y_, x_] =
Interpolation[dataset1, InterpolationOrder -> 1][Log10[y], Log10[x]]
xmin = 10^-1.98;
xmax = 10^4.98;
ymin = 10^-7.28;
ymax = 10^0.497;
xMax[y_] =
Piecewise[{{4000, 0 <= y < 10^-3}, {3000,
10^-3 <= y < 2.5*10^-3}, {2000,
2.5*10^-3 <= y < 3.5*10^-3}, {1500,
3.5*10^-3 <= y < 7.5*10^-3}, {1000,
7.5*10^-3 <= y < 0.025}, {500, 0.025 <= y < 0.04}, {300,
0.04 <= y < 0.07}, {200, 0.07 <= y < 0.1}, {100,
0.1 <= y < 0.25}, {50, 0.25 <= y < 0.5}, {30,
0.5 <= y < 0.8}, {20, 0.8 <= y < Pi/2}}];
Here, xMax[y]
denotes the maximal value of x for the given y (I have obtained it by hand).
Next, I want to compute an integral over x and y within some interval $ymin<y<ymax$. I introduce the following functions:
ySpectrum[method_, yminv_, ymaxv_] :=
NIntegrate[yxDistr[y, x], {y, yminv, ymaxv}, {x, xmin, 4000},
Method -> method]
ySpectrum2[method_, yminv_, ymaxv_] :=
NIntegrate[yxDistr[y, x], {y, yminv, ymaxv}, {x, xmin, xMax[y]},
Method -> method];
If using method="InterpolationPointsSubdivision", ySpectrum["InterpolationPointsSubdivision",0.2,0.3]
gives zero, while ySpectrum["AdaptiveMonteCarlo",0.2,0.3]
is non-zero (as it should be, if looking at the initial data). I may improve the accuracy of the prediction of the "AdaptiveMonteCarlo" by calling
ySpectrum2["AdaptiveMonteCarlo",0.2,0.3]
.
But how to make "InterpolationPointsSubdivision" work correctly? It does not work properly with non-rectangular regions of integration.
dat = GatherBy[data, First]; 0.05^2*(2 Total[dat[[All, All, 3]], 2] + 2 Total[dat[[2 ;; -2, 2 ;; -2, 3]], 2] - Total[dat[[{1, -1}, {1, -1}, 3]], 2])/4
$\endgroup$