Approximating unit circle with a cubic spline

Hey all. We're covering cubic splines in class and I've been messing around in Matlab a little bit and I've run into a question I can't seem to answer. I'm trying to approximate the unit circle with cubic splines. I've written the following function to do this. It consumes two row vectors, x and y. Using the built-in spline function, I've parameterized x and y with respect to a variable t that I let range from 0 to 1. x and y have two extra elements compared to t to allow for the use of endslopes. My problem is that I'm not getting a symmetrical curve and I'm wondering why? Does the error stem from the improper use of the spline function or is there a mathematical reason for this that I'm currently missing? Here's the relevant code:

function parametric_spline(x, y)
% Draws a parametric curve (x(t), y(t)

t = linspace(0,1,length(x)-2);
tplot = linspace(0,1,100);

csx = spline(t, x, tplot);
csy = spline(t, y, tplot);

figure(1)
plot(csx,csy)
hold on
plot(cos(2*pi*tplot),sin(2*pi*tplot), 'color', 'red')
xlim([-2 2]);
ylim([-2 2]);
axis equal

figure(2)
plot(tplot, csx);
hold on
plot(tplot, cos(tplot*2*pi), 'color', 'red');
xlim([0 1]);
ylim([-2 2]);

figure(3)
plot(tplot, csy);
hold on
plot(tplot, sin(tplot*2*pi), 'color', 'red');
xlim([0 1]);
ylim([-2 2]);

end

As input, I've used:

x = [0, 1, 0, -1, 0, 1, 0]; y = [1, 0, 1, 0, -1, 0, 1];

I've included the figures produced. Is there a reason the fit is so much better to cos(t) than to sin(t)?

EDIT:

As a followup, I tried setting the first and last values in y as 5 and 5 - it seems to work much better now. I can't see why, however. The slope of y(t) at the beginning and end should be 1, should it not?

Answer 1

[John D'Errico]

You were close. First, I've be careful about the parameter t. Since the points are equally spaced around the circle, there is no problem. But if the points were not equally spaced in angle around the circle, then t must reflect that, else it would create artifacts in the curve shape. (My recommendation is usually to use a cumulative chordal arc length, so the cumulative piecewise linear distance along the connect-the-dots curve between the points.) I'll use your parameter t here though.

x = [0, 1, 0, -1, 0, 1, 0]; y = [1, 0, 1, 0, -1, 0, 1];
t = linspace(0,1,length(x)-2);

The problem is what have you provided for information about the end slopes?

For the x curve, x is essentially cos(2*pi*t). So it has zero end slopes. You provided a vector x that has first and last elements as zero, so that suggests the cosine curve should fit very nicely.

Now, lets look at the y curve. What is the derivative of sin(2*pi*t) with respect to t? It would be:

2*pi*cos(2*pi*t)

So at 0 and 1 for t, we would expect a derivative of 2*pi, NOT 1. So this should work better:

y = [2*pi, 0, 1, 0, -1, 0, 2*pi];

I'll let you verify my claim, but in my own tests, it worked very nicely.

[Mark Soric]

This is definitely a facepalm moment for me - I completely forgot about the scaling there. Thanks John!

Followup question about the spacing of t. I was thinking of trying to calculate the curvature of the cubic between each point and space t values according to the magnitude of the curvature.

Is there a built-in Matlab function that would create a vector t from 0 to 1 where the spacing isn't linear?

[John D'Errico]

The clue about what you did was to look at the plots. The sine wave one was wrong, but not the cosine one. And the sine curve was wrong mainly at the ends. So I immediately knew where to look. And, well, after using splines since the dark ages, I've already made all possible mistakes.

Of course, you can also use tools like my interparc to interpolate general curves like this. You seem to be making good progress on your own though, and writing code like this is a great way to learn. The reason I point out interarc is that tool tries to locate points that are equally spaced in distance along a curve. I'm not sure if that relates to what you are trying to do though, and I really do like to see people wanting to learn on their own. Getting your hands dirty is the very best way to learn.

If the points are not perfectly uniformly spaced as you have, the nice way to create a parametric t (assuming row vectors x and y):

t = cumsum([0,sqrt(diff(x).^2 + diff(y).^2)]);
t = t/t(end);

You can see what it does: forming the distance between each pair of points, then using cumsum.

A problem with using known end conditions as you did, where you force the slope to a known value, is you needed to know the slope in advance. Get it wrong and as you saw, the curve looks a bit wonky.

I suppose you might have used a spline where the sine wave used zero end conditions on the second derivative at those ends, known as a "natural" cubic spline. The problem there is it would have failed to work well on the cosine curve, as there the second derivative is definitely non-zero at the ends.

The Carl de Boor trick, of using not-a-knot end conditions is often a nice one. (Enforce third derivative continuity at the second and next to last knots.)