1
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Maybe a manual trapezoidal rule: 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$
    – Michael E2
    Commented Nov 26, 2020 at 18:07

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.