3
$\begingroup$

I have a set of data

table = Reverse[{
   {0, -947}, {-212, -947}, {-424, -950}, {-635, -963}, {-845, -995}, {-1051, -1044},
   {-1248, -1119}, {-1432, -1224}, {-1591,-1365}, {-1715, -1537}, {-1796, -1732}, 
   {-1828, -1942}, {-1810,-2153}, {-1755, -2357}, {-1668, -2551}, {-1556, -2730}, 
   {-1423,-2895}, {-1278, -3050}, {-1126, -3197}, {-973, -3343}, {-803,-3506}
}];

which I need to interpolate with a parametric curve. Endpoints and tangent vectors at them are critical. Trying independent polynomial fits for X and Y leads to a wavy curve which is unacceptable. BSplineFunction of degree around 3 leads to a very good-looking curve:

interp = BSplineFunction[table, SplineDegree -> 3];    
Show[ListPlot[table], ParametricPlot[interp[t], {t, 0, 1}], AspectRatio -> 1]

enter image description here

but has a discontinuous derivative (btw: why? help promises that it should be a BSplineFunction of one lesser degree but it's a scary jump-step function if you plot it):

x1[t_?NumericQ] := Module[{val}, val = interp[t]; First@val]
y1[t_?NumericQ] := Module[{val}, val = interp[t]; Last@val]    
Plot[Norm[{x1'[t], y1'[t]}], {t, 0, 1}, PlotRange -> All]

enter image description here

I then intend to use these functions in NDSolve, NIntegrate etc., therefore behaviour like above is too bad. What parametric interpolation is better in this case?

$\endgroup$
5
  • $\begingroup$ What are x1 and y1? $\endgroup$
    – VLC
    Commented Feb 7, 2013 at 15:21
  • $\begingroup$ @VLC check the x and y in mathematica.stackexchange.com/questions/19229/… $\endgroup$ Commented Feb 7, 2013 at 15:41
  • $\begingroup$ @PlatoManiac Ah, ok. $\endgroup$
    – VLC
    Commented Feb 7, 2013 at 15:55
  • $\begingroup$ By extracting the values from the BSplineFunction and then differentiating, you are taking a finite difference derivative sampled at whatever points Plot feels like using. Try e.g. Plot[Norm[interp'[t]], {t, 0, 1}], which I think you will find much better. $\endgroup$ Commented Feb 8, 2013 at 1:41
  • $\begingroup$ You may want to look at this method mathematica.stackexchange.com/a/33262/1364 It works for me. $\endgroup$ Commented May 29, 2014 at 20:08

2 Answers 2

4
$\begingroup$

Borrowing the parametrizeCurve[] routine from this answer (it's amazing what searching can do...), you can do this:

tvals = parametrizeCurve[table];
{ft, gt} = Interpolation[Transpose[{tvals, #}], Method -> "Spline"] & /@ Transpose[table];

ParametricPlot[{ft[u], gt[u]}, {u, 0, 1}, Axes -> None, Frame -> True,
               Epilog -> {AbsolutePointSize[4], Magenta, Point /@ table}]

curve and points

Verify the $C^2$ continuity:

{Plot[{ft'[u], gt'[u]}, {u, 0, 1}], Plot[{ft''[u], gt''[u]}, {u, 0, 1}]} // GraphicsRow

derivative plots

$\endgroup$
2
  • $\begingroup$ There turns out to be some substantial wavering at the lower end, so I had to come up with some crazy Fit with basis {1,t,t^2,...,t^6,Exp[0.05 t],Cos[t]} which magically led to incredibly smooth curve with constant-signed curvature. $\endgroup$ Commented Feb 20, 2013 at 14:14
  • $\begingroup$ All this is so important because of its architectural nature. $\endgroup$ Commented Feb 20, 2013 at 14:15
2
$\begingroup$

(1) I get a reasonable result interpolating separately.

(2) I don't think your derivatives are doing what you expect.

x1 = Interpolation[table[[All, 1]]];
y1 = Interpolation[table[[All, 2]]];

Show[ListPlot[table], ParametricPlot[{x1[t], y1[t]}, {t, 1, 21}], 
 AspectRatio -> 1]

enter image description here

Here are those derivatives.

Plot[{x1'[t], y1'[t]}, {t, 0, 21}, PlotRange -> All]

enter image description here

If you raise the InterpolationOrder then those derivative curves get smoother.

$\endgroup$
1
  • $\begingroup$ Problem with Interpolation is that as soon as no points are close to each other, everything seems fine. But if data is not equidistant like this, the curve breaks between close points - something that doesn't happen with BezierFunction for example (though it's not an interpolation at all). $\endgroup$ Commented Feb 8, 2013 at 6:30

Your Answer

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

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