5
$\begingroup$

I have a data, for example

data = {{0, 0.443395}, {0.0171663, 0.443395}, {0.0465206, 
      0.440663}, {0.0727851, 0.432468}, {0.0959596, 0.412555}, {0.116044, 
      0.37872}, {0.133039, 0.332491}, {0.146943, 0.276475}, {0.157758, 
      0.213267}, {0.165483, 0.145097}, {0.170118, 0.0741668}, {0.171663, 
      3.00302*10^-9}};

This is not a very smooth dataset (kinda few points, not terribly few, but it could use more mid points), so I wanted to interpolate it. When I did, something strange happened

Plot[Interpolation[data,InterpolationOrder->2][x],{x,0,0.171663}]

bump in the function There is a noticeable bump in the plot around x = 0.165. It looks like the interpolating function is not differentiable around that point (that is strange, as the differentiability is on of the basic assumptions when dealing with interpolations of order 2 and higher). What causes this behaviour? Dataset itself looks completely "normal" everywhere

ListPlot[f]

ListPlot of dataset

What causes this behaviour? What can I do to fix this problem?

$\endgroup$
  • 3
    $\begingroup$ It loos somewhat better if you use Method -> "Spline". $\endgroup$ – user46676 Apr 9 '17 at 20:42
  • $\begingroup$ But how can I interpolate it with this algorithm: for argument in some range between two points (x1, x2), the function is ax^2+bx+c. Determine the coefficients from three equations: value of interpolation must match function values in points x1, x2 and derivative of interpolation must be continuous. So to have interpolation of order two, one has to provide at least three points. Derivative at the ending points can be set by hand. For me it's obvious f'(0) = 0, f'(end of the interval) = infinity (like for function sqrt(1-x)). $\endgroup$ – user16320 Apr 9 '17 at 21:00
  • $\begingroup$ With Hermite interpolation of order $n$, the derivative of order $n-1$ is often not continuous, I believe. (It's certainly true for $n = 3$.) Each segment is the unique quadratic polynomial through three adjacent points. You get an obvious kink because of the data. There are less obvious kinks at the other points. $\endgroup$ – Michael E2 Apr 9 '17 at 21:01
  • $\begingroup$ You can specify the derivatives at the points -- see the documentation. $\endgroup$ – Michael E2 Apr 9 '17 at 21:02
  • 1
    $\begingroup$ The reason I asked about fitting is it is a natural procedure for data with error or noise present. If your data are like that it is ill to constrain the function to pass through points exactly - it will result in those bumps. Usually what people do they find an analytic model (you can do better than that formula). Besides fitting the points it also gives you insight into the law behind the data. So it is conceptual question: do you need to fit or interpolate? $\endgroup$ – Vitaliy Kaurov Apr 9 '17 at 21:08
2
$\begingroup$

If your data comes from numerical simulation, then at a fine scale it can be quite noisy, due to iterative solvers and stopping criteria, etc. This is a common issue, for example, with gradient based solvers that are solving a problem with simulations as their kernel. The first derivative of the simulation results can be noisy, and so you fit the data in proximity.

I'd use the spline fit to ensure a smooth curve. That's my answer. But...

MMa used to have a radial basis function capability, which is an excellent interpolation method, but it seems to be gone (although it still pops up in an on-line search of MMa documentation).

I thought the machine learning Predictor[ ] function would get it, as it has a "GaussianProcess" option, but I couldn't figure out how to tell the predictor to assume a deterministic model. Still, here are the results... start by massaging the data into the form the predictor needs:

    pdata = Map[#[[1]] -> #[[2]] &, data];

    {0 -> 0.443395, 0.0171663 -> 0.443395, 0.0465206 -> 0.440663, 
     0.0727851 -> 0.432468, 0.0959596 -> 0.412555, 0.116044 -> 0.37872, 
     0.133039 -> 0.332491, 0.146943 -> 0.276475, 0.157758 -> 0.213267, 
     0.165483 -> 0.145097, 0.170118 -> 0.0741668, 
     0.171663 -> 3.00302*10^-9}

Create predictor

    p = Predict[pdata, Method -> "GaussianProcess", PerformanceGoal -> "Quality"];

Get info

    PredictorInformation[p]

Note the "AssumeDeterministic" is False. Tried to set it to "True" but can't.

Output of info on predictor

Plot it...

    Show[ListPlot[data], Plot[p[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]]

enter image description here

* EDIT *

If you have some physics or other knowledge to tie to the data, don't throw it away. For example, your curve looks like a really good fit to

    fn[x_]:=0.443*Sqrt[1 - (x/.17167)^3.3]

enter image description here

So use the interpolation to model the difference between fn[x] and the data. I did, got this plot:

    del = Map[{#[[1]], (fn[#[[1]]] - #[[2]]) // Re} &, data];
    int = Interpolation[del, Method -> "Spline"];
    Show[ListPlot[data], Plot[fn[x] - int[x], {x, 0, .171663}]]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Haha, yes, there is some physics tied to it, however the physics is really just a mean-field theory (that means critical exponent is 1/2 for the order parameter). Now I realized that this doesn't mean that there should be term T/Tc under the square root! You've lead me to this realization, thank you. $\endgroup$ – user16320 Apr 11 '17 at 18:43
2
$\begingroup$

The data looks smooth enough. The glitch is due to the parabola with a vertical axis that fits the last three points. What if the data was rotated? The glitch should go away. But, rotated by what angle and around which point?

Suppose we have the abscissa $x_0$ and we want the interpolant $y_0$. Then the point we will rotate the data about will be $(x_0,0)$ and the angle will be 45 degrees, CCW. This choice makes the math simple enough to demonstrate feasibility of this approach.

Here's the code for a rotating interpolation function and the code for two test plots.

rifn[x0_, pts_] := 
 Block[{\[Theta] = \[Pi]/4, fwd, rot, ifn, pt0, x1, x, y1},
  fwd = {{Cos[\[Theta]], -Sin[\[Theta]]}, {Sin[\[Theta]], 
     Cos[\[Theta]]}};
  rot = Dot[fwd, # - {x0, 0}] & /@ data;
  ifn = Interpolation[rot, InterpolationOrder -> 2];
  x1 = x /. First[FindRoot[ifn[-x] == x, {x, 0}] // Quiet];
  y1 = ifn[-x1];
  pt0 = fwd.{x1, y1} + {x0, 0};
  Last[pt0]
  ]
Plot[rifn[x, data], {x, 0, 0.171663},
 Epilog -> {Red, Point@data}, AspectRatio -> GoldenRatio]

Plot[{rifn[x, data],
  Interpolation[data, InterpolationOrder -> 2][x]},
 {x, 0.15, 0.171663},
 Epilog -> {Red, Point@data}, AspectRatio -> GoldenRatio]

The first plot shows the smoothness of interpolation over the entire range of the data. The second plot zooms in on an interval that contains the glitch. The rotating interpolation function appears to be feasible. Sketch

How does it work? First, we define the rotation matrix fwd that rotates the data CCW. The second line takes each point in the data, shifts is by the amount {x0,0} and the applies the rotation matrix. The third line uses the built-in Interpolation function. The fourth line requires a little sketch to explain.

The sketch attempts to show 4 data points, the abscissa $x_0$ and the desired interpolant $y_0$. When the data points are rotated the new abscissa will be $x_1$ and the new interpolant will be $y_1$. Since we are rotating 45 degrees, we will have $x_1=y_1$. We use FindRoot to find the abscissa that satisfies this condition. Thus, the angle is hardwired into the argument of FindRoot.

Having found $x_1$, we can find $y_1$ either by interpolation or by geometry. Finally, we rotate the point $(x_0,y_0)$ back (using the same fwd rotation matrix) and shift it back. Thus, inside the rifn[] function, pt0 is the interpolated point $(x_0,y_0)$ that we were looking for. The rifn function returns only $y_0$.

This approach appears to be quite feasible for the given data.

$\endgroup$
  • $\begingroup$ Interesting approach, very similar to MichaelE2's approach (squaring the function, interpolating that and square rooting the result), which I find a little bit easier. $\endgroup$ – user16320 Apr 11 '17 at 18:45

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.