How to make that interpolation?

Question

I have the following dataset:

SetDirectory[NotebookDirectory[]];
TableSpectrum = Import["spectrum.dat", "Table"];

I want to make the interpolation. However, when simply using

Interpolation[TableSpectrum, InterpolationOrder -> 1]

it writes that the quality of the underlying mesh is too low. The solution is to change the data in the following way:

Spectrum1[mS_, ES_] = 
 Interpolation[{Log[#[[1]]], #[[2]], #[[3]]} & /@ TableSpectrum, 
   InterpolationOrder -> 1][Log[mS], ES]
Spectrum2[mS_, ES_] = 
 Interpolation[{#[[1]], Log[#[[2]]], #[[3]]} & /@ TableSpectrum, 
   InterpolationOrder -> 1][mS, Log[ES]]

However, they each have different problems. Spectrum1 is very bad interpolation:

PlotmS[mS_] := LogLogPlot[{Spectrum1[mS, ES]}, {ES, mS, 5000}]
PlotmS[20]

While Spectrum2 is integrated much more slowly:

Timing[NIntegrate[Spectrum1[20, ES], {ES, 20, 5000},Method->"AdaptiveMonteCarlo"]]
Timing[NIntegrate[Spectrum2[20, ES], {ES, 20, 5000},Method->"AdaptiveMonteCarlo"]]]

Is there any way to make the interpolation including the best properties of Spectrum2 (smooth interpolation between grid) and Spectrum1 (fast integration)?

@MichaelE2 : I agree that this is not quite accurate to use Timing. However, the result does not change when using AbsoluteTiming. — John Taylor, Jul 29 '19 at 15:34

halirutan · Accepted Answer · 2019-07-30T15:42:43.493

If you look at your data, you'll find that the sampling is vastly different in x- and y-direction. We can show that by splitting your list of points depending on the first dimension:

data = Import["~/Downloads/spectrum.dat", "Table"];
splitted = SplitBy[data, First];

Mathematica graphics

I'm a bit unsure, what your final goal is as it seems that you try to integrate along these lines. If this is the sole purpose of the interpolation, then you might want to consider a simpler approach.

For instance, at the end of your question, you calculate the area under the curve at x=20. The sampling distance in y-direction is (almost) always 3.0, so a far simpler method is to use a Riemann sum and calculate the integral by

3.0*Total[splitted[[5, All, 3]]]
(* 5.50244*10^-6 *)

Alternatively, you can implement the Trapezoidal rule and you can quickly integrate all of your lines and also take into account when the distance between two points is not 3.0:

intC = Compile[{{pts, _Real, 2}},
   Module[{
     x1 = pts[[1, 1]],
     y1 = pts[[1, 2]],
     x2 = pts[[2, 1]],
     y2 = pts[[2, 2]]
     },
    (x2 - x1)*(y1 + y2)/2.0
    ],
   RuntimeAttributes -> Listable,
   CompilationTarget -> "C",
   Parallelization -> True];

integrateLine[pts_?MatrixQ] := Total[intC[Partition[pts[[All, {2, 3}]], 2, 1]]]

The result for the fifth line at x=20 is almost unchanged

integrateLine[splitted[[5]]]
(* 5.50241*10^-6 *)

And of course, you can simply map integrateLine over all lines in splitted and get an estimate of how the area changes along x:

intData = Transpose[
   {splitted[[All, 1, 1]], integrateLine /@ splitted}];
ListLinePlot[intData,
 PlotRange -> All,
 Mesh -> All,
 MeshStyle -> Red,
 Frame -> True,
 FrameLabel -> {"x", "Area of curve"},
 GridLines -> Automatic
 ]

Mathematica graphics

I'm not sure this is what you are really after, but maybe it helps.

How to make that interpolation?

1 Answers1