1
$\begingroup$

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]

enter image description here

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

enter image description here

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

$\endgroup$
  • $\begingroup$ BTW, why use Timing[]? $\endgroup$ – Michael E2 Jul 29 at 14:55
  • $\begingroup$ @MichaelE2 : I agree that this is not quite accurate to use Timing. However, the result does not change when using AbsoluteTiming. $\endgroup$ – John Taylor Jul 29 at 15:34
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.