Skip to main content

plotting - Joining and interpolating data points



I wonder what is the best practice for interpolating curves? Usually, I'm using BSplineCurve and adjusting SplineWeights so it would fit better (and assigning more weight around the sharp edges to drag the curve closer to it). Or if I can guess what formula describes the points, I use FindFit. But often, I can't guess the formula and adjusting weights is very tedious, so it's easier to just manually draw the curve along the points. So what is the best way to join points in Mathematica?


Consider these five sets of points describing five resonance curves for example:


 data = {{{100.434, 0.}, {102.461, 0.0909091}, {104.392, 0.318182}, {105.321, 0.545455},
{106.226, 1.}, {107.108, 0.545455}, {107.965, 0.318182}, {109.608, 0.136364},
{111.154, 0.0909091}},
{{100.434, 0.}, {102.461, 0.06}, {104.392, 0.22}, {105.321, 0.46}, {106.226, 1.},
{107.108, 0.4}, {107.965, 0.12}, {111.154, 0.02}, {113.958, 0.}},
{{100.434, 0.030303}, {102.461, 0.0505051}, {104.392, 0.0909091},
{105.321, 0.272727}, {105.867, 0.494949}, {106.226, 1.}, {106.582, 0.636364},
{107.108, 0.212121}, {107.965, 0.0505051}, {111.154, 0.}},

{{100.434, 0.}, {102.461, 0.0555556}, {104.392, 0.333333}, {105.321, 0.611111},
{105.867, 1.}, {106.226, 0.944444}, {106.405, 0.583333}, {106.582, 0.777778},
{106.933, 1.}, {107.965, 0.444444}, {109.608, 0.166667}, {111.154, 0.0555556}},
{{100.434, 0.0188679}, {102.461, 0.0566038}, {104.392, 0.}, {105.321, 0.54717},
{105.867, 0.849057}, {106.226, 1.}, {106.405, 0.433962}, {106.582, 0.886792},
{106.933, 0.924528}, {107.281, 0.660377}, {107.965, 0.320755},
{111.154, 0.0566038}}}

Answer




So, what is the best way to join points in Mathematica?




There is no one "best way" (not only in Mathematica, but in general); an interpolation scheme that behaves nicely for data set A might be a crapshoot when applied to data set B. It depends on the configuration of your points, and impositions you have on the interpolant (e.g. $C^1$/$C^2$ continuity, preservation of monotonicity/convexity, etc.), with these impositions not always being achievable all at the same time.


Having said this, if you're fine with a $C^1$ interpolant (interpolant and first derivative are continuous), one possibility is to use Akima interpolation. It is not always guaranteed to preserve shape (unless your points are specially configured), but at least for your data set, it does a decent job:


AkimaInterpolation[data_] := Module[{dy}, dy = #2/#1 & @@@ Differences[data];
Interpolation[Transpose[{List /@ data[[All, 1]], data[[All, -1]],
With[{wp = Abs[#4 - #3], wm = Abs[#2 - #1]},
If[wp + wm == 0, (#2 + #3)/2, (wp #2 + wm #3)/(wp + wm)]] &
@@@ Partition[ArrayPad[dy, 2, "Extrapolated"], 4, 1]}],
InterpolationOrder -> 3, Method -> "Hermite"]]


MapIndexed[(h[#2[[1]]] = AkimaInterpolation[#1]) &, data];

Partition[Table[Plot[{h[k][u]}, {u, 100.434, 111.154}, Axes -> None,
Epilog -> {Directive[Red, AbsolutePointSize[4]], Point[data[[k]]]},
Frame -> True, PlotRange -> All],
{k, 5}], 2, 2, 1, SpanFromBoth] // GraphicsGrid

Akima interpolants for resonance curves


Note that in the fifth plot, the Akima interpolant has a slight wiggle before hitting the third point; this, as I said, is due to the fact that Akima's scheme does not guarantee that it will respect the monotonicity of the data. If you want something that fits a bit more snugly, one scheme you can use is Steffen interpolation, which is also a $C^1$ interpolation method:


SteffenEnds[h1_, h2_, d1_, d2_] :=

With[{p = d1 + h1 (d1 - d2)/(h1 + h2)}, (Sign[p] + Sign[d1]) Min[Abs[p]/2, Abs[d1]]]

SteffenInterpolation[data_?MatrixQ] := Module[{del, h, pp, xa, ya},
{xa, ya} = Transpose[data]; del = Differences[ya]/(h = Differences[xa]);
pp = MapThread[Reverse[#1].#2 &, Map[Partition[#, 2, 1] &, {h, del}]]/
ListConvolve[{1, 1}, h];
Interpolation[Transpose[{List /@ xa, ya,
Join[{SteffenEnds[h[[1]], h[[2]], del[[1]], del[[2]]]},
ListConvolve[{1, 1}, 2 UnitStep[del] - 1] *
MapThread[Min, {Partition[Abs[del], 2, 1], Abs[pp]/2}],

{SteffenEnds[h[[-1]], h[[-2]], del[[-1]], del[[-2]]]}]}],
InterpolationOrder -> 3, Method -> "Hermite"]]

MapIndexed[(w[#2[[1]]] = SteffenInterpolation[#1]) &, data]

Partition[Table[Plot[{w[k][u]}, {u, 100.434, 111.154}, Axes -> None,
Epilog -> {Directive[Red, AbsolutePointSize[4]], Point[data[[k]]]},
Frame -> True, PlotRange -> All],
{k, 5}], 2, 2, 1, SpanFromBoth] // GraphicsGrid


Steffen interpolants for resonance curves


Note that the interpolants from Steffen's method are a lot less wiggly, though the interpolant turns more sharply at its extrema. The advantage of using Steffen is that it is guaranteed to preserve the shape of the data, which might be more important than a smooth turn in some applications.


My point, now, is that sometimes, one must try a number of interpolation schemes to see what is most suitable for the data at hand, for which plotting the interpolant along with the data it is interpolating is paramount.


Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...