%% Splines and the geometry of curves %% Introduction % % Cubic splines are piecewise polynomials which are continuous and have continuous % derivatives. In this lab, we will investigate how to use the Matlab % functions pchip and spline. As we will see, pchip is good for interpolating physical data which % should be constrained in some way (should always be positive, or between 0 % and 1, for example) while the spline function is ideally % suited for geometric data. % clear all; close all; format short; %% Single polynomial interpolant (polyfit) % % Before starting with splines, we illustrate how to fit a single % polynomial through a set of data points. Recall that we can uniquely % interpolate N+1 data points with a polynomial of degree N. % As an example, suppose we wish to put a smooth curve through a set of % random data points. % N = 13; xdata = linspace(-5,5,N+1); ydata = rand(1,N+1); xs = linspace(-5,5,500); plot(xdata,ydata,'k.','markersize',30); hold on; plot(xdata,0*xdata,'k--'); plot(xdata,0*xdata + 1,'k--'); ylim([-1 2]); %% % % We first construct the polynomial using the polyfit function. This function returns the % coefficients of the interpolating polynomial. We then evaluate the % polynomial using polyval. % p = polyfit(xdata,ydata,N); % Coefficients of polynomial of degree (N-1). ys = polyval(p,xs); hold on; plot(xs,ys,'b','linewidth',2); plot(xdata,ydata,'.k','markersize',30); plot(xdata,0*xdata,'k--'); plot(xdata,0*xdata + 1,'k--'); title('Interpolating function (polyfit)','fontsize',18); ylim([-1 2]); %% Piecewise Polynomial (spline) % % We can improve the behavior of the single polynomial interpolant by % adjusting the location of the points in xdata. % But often, we do not have a choice as to where to place the data points. % A second option is to construct a piecewise polynomial. A common % approach is to put a cubic polynomial over each interval defined % by the N data points. The function spline is an % example of piecewise polynomial interpolant. We can simply use the % spline command with a single call. For example % ys = spline(xdata,ydata,xs); clf; plot(xs,ys,'r','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); plot(xdata,0*xdata,'k--'); plot(xdata,0*xdata + 1,'k--'); title('Interpolating function (spline)','fontsize',18); ylim([-1 2]); %% % % Notice that the curve appears to have fewer oscillations than the % single polynomial case. % %% % % As with the polyfit, we can also call the spline using a two step procedure. First, we compute % the coefficients of the polynomial using the spline function, then we evaluate the polynomial using % ppval. % pp_spline = spline(xdata,ydata) %% % % The pp_spline is a Matlab structure that % stores several pieces of information. In this case, the structure stores % the break points (defined by xdata), the % coefficients of each polynomial in the segments defined by the % breakpoints, the dimension of the polynomial, the number of pieces, and % the order. %

% To evaluate the piecewise polynomial, we use ppval % ys = ppval(pp_spline,xs); %% % % This would produce the same y values that we computed above with the % single call to spline. % %% First derivative (polyfit) % % To compute derivatives of the polynomial interpolant, we need to be able % to compute derivatives of a polynomial described by a set of coefficients. % For example, to compute the derivative of the polynomial obtained using % polyfit, we use polyder as % clf; p = polyfit(xdata,ydata,N-1); p_deriv = polyder(p); ys_deriv = polyval(p_deriv,xs); ydata_deriv = polyval(p_deriv,xdata); plot(xs,ys_deriv,'b','linewidth',2); hold on; plot(xdata,ydata_deriv,'k.','markersize',30); title('First derivative (polyfit)','fontsize',18); ylim([-10 10]); %% % % Notice the scale on the y-axis. The derivative is quite large. % %% First derivative (spline) % % To compute the derivative of a piecewise polynomial, we first need % access to the coefficients of each piecewise segment. These % coefficients can be accessed from the structure pp_spline. % pp_spline.coefs %% % % The rows of this array are the coefficients of the polynomial segments. % The form that each polynomial takes is given by % $$% p_k(x) = a_3(x-x_k)^3 + a_2(x-x_k)^2 + a_1(x-x_k) + a_0 %$$ % where the $x_k$ are the breakpoints % typically defined by the x data points. % We then evaluate the interpolant using ppval. %

% Using the structure pp_spline in combination with % polyder, we can compute the coefficients of the % of each the derivative of each % polynomial segment. Using these coefficients, we can then create a new % piecewise polynomial structure which we can easily evaluate. % % loop over each interval to get the derivatives in that interval. coeffs_deriv = zeros(N,3); % store the coefficients of the derivative for k = 1:N, c = pp_spline.coefs(k,:); % pp is a structure coeffs_deriv(k,:) = polyder(c); end pp_deriv = mkpp(xdata,coeffs_deriv); ys_deriv = ppval(pp_deriv,xs); ydata_deriv = ppval(pp_deriv,xdata); clf; plot(xs,ys_deriv,'r','linewidth',2); hold on; plot(xdata,ydata_deriv,'k.','markersize',30); plot(xdata,0*xdata,'k--'); title('First derivative (spline)','fontsize',18); ylim([-2 2]); %% Second derivative (polyfit) % % Continuing in this fashion, we can plot the second derivatives of each of % these interpolants. % clf; p = polyfit(xdata,ydata,N-1); p_deriv2 = polyder(polyder(p)); ys_deriv2 = polyval(p_deriv2,xs); ydata_deriv2 = polyval(p_deriv2,xdata); plot(xs,ys_deriv2,'b','linewidth',2); hold on; plot(xdata,ydata_deriv2,'k.','markersize',30); title('Second derivative (polyfit)','fontsize',18); ylim([-10 10]); %% Second derivatives (spline) % % The cubic spline has a continuous second derivative. % % loop over each interval to get the derivatives in that interval. coeffs_deriv2 = zeros(N,2); for k = 1:N, c = pp_deriv.coefs(k,:); % pp is a structure coeffs_deriv2(k,:) = polyder(c); end pp_deriv2 = mkpp(xdata,coeffs_deriv2); ys_deriv2 = ppval(pp_deriv2,xs); ydata_deriv2 = ppval(pp_deriv2,xdata); clf; plot(xs,ys_deriv2,'r.','markersize',10); hold on; plot(xdata,ydata_deriv2,'k.','markersize',30); title('Second derivative (spline)','fontsize',18); ylim([-10 10]); %% The Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) % % While both the single polynomial and the spline function produce smooth % curves that interpolate our data, we may not be so happy with the fact % that in some cases, the curve may introduce new extrema that were not % present in the original data. This is especially problematic if the data % is only physically meaningful in a limited range, say, [0,1]. In this % situation, Matlab provides a useful command that prevents the % introduction of new extrema into the data. This is called the pchip and can be used in a manner very similar to the % spline comand. Below, we illustrate the pchip % interpolant and the first two derivatives. You'll notice that the second % derivative is no longer continuous. % ys = pchip(xdata,ydata,xs); clf; plot(xs,ys,'g','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); plot(xlim,[0 0],'k--'); plot(xlim,[1 1],'k--'); title('Interplating function (pchip)','fontsize',18); ylim([-1 2]); %% % % Notice that all of the values are between 0 and 1. The values % returned from evaluating a pchip curve will % never exceed range of the original data. % %% First derivative (PCHIP) % % The pchip spline, unlike the spline derivative, % uses derivatives that are computed strictly from the data. Derivative % values at extrema are set to 0, whereas derivatives at intermediate data % points are taken to be a weighted average of the two one-sided slopes % based on y-data. %

% First, we have to create a spline for the first and second derivatives. % pp_pchip = pchip(xdata,ydata); coeffs_deriv = zeros(N,3); coeffs_deriv2 = zeros(N,2); for k = 1:N, c = pp_pchip.coefs(k,:); coeffs_deriv(k,:) = polyder(c); coeffs_deriv2(k,:) = polyder(polyder(c)); end pp_deriv = mkpp(xdata,coeffs_deriv); pp_deriv2 = mkpp(xdata,coeffs_deriv2); %% % % Then we can evaluate and plot these splines representing the derivatives. % clf; ys_deriv = ppval(pp_deriv,xs); ydata_deriv = ppval(pp_deriv,xdata); plot(xs,ys_deriv,'g','linewidth',2); hold on; plot(xdata,ydata_deriv,'k.','markersize',30); plot(xlim,[0 0],'k--'); title('First derivative (pchip)','fontsize',18); %% Second derivative (PCHIP) % % The second derivative of the PCHIP spline is not continuous, as you can % see from this graph % clf; ys_deriv2 = ppval(pp_deriv2,xs); ydata_deriv2 = ppval(pp_deriv2,xdata); plot(xs,ys_deriv2,'g.','markersize',10); hold on; plot(xdata,ydata_deriv2,'k.','markersize',30); title('Second derivative (pchip)','fontsize',18); %% Computing tangents to the curve % % We can compute tangents to the curve using a formula for the tangent. % For a function, we assume the data points have the form % ${\bf r}(x) = (x,f(x))$. Then we can % apply the formula for the tangent, given by % $$% {\bf T}(x) = \frac{{\bf r}'(x)}{g(x)} %$$ % where % $$% g(x) = \parallel \!\! {\bf r}'(x) \!\! \parallel = \sqrt{1 + f'(x)^2} %$$ % We will test this out on the interpolating spline. % f_spline = spline(xdata,ydata); coeffs_deriv = zeros(N,3); coeffs_deriv2 = zeros(N,2); for k = 1:N, c = f_spline.coefs(k,:); coeffs_deriv(k,:) = polyder(c); coeffs_deriv2(k,:) = polyder(polyder(c)); end f_deriv = mkpp(xdata,coeffs_deriv); f_deriv2 = mkpp(xdata,coeffs_deriv2); ys = ppval(f_spline,xs); clf; plot(xs,ys,'r','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); %% % % Now we will use the formulas to compute the tangents to the curve % at our data points. % % Compute tangent vector r = [xdata; ydata]; rp = [ones(1,N+1); ppval(f_deriv,xdata)]; g = sqrt(dot(rp,rp)); T = bsxfun(@times,rp,1./g); %% % % We can plot the tangents as arrows using the quiver command. % quiver(xdata,ydata,T(1,:),T(2,:),0.25,'k','linewidth',2); hold on; plot(xdata,ydata,'k.','markersize',30); title('Tangents (spline)','fontsize',18); xlim([-5 5]); ylim([-.5 1.5]); %% Lab exercises % %
%
%
1. Use the functions polyval and polyder to plot the following polynomial and its % derivative. % $$p(x) = 5x^3 - 2x^2 + 4x - 10$$ %
2. Use the functions polyval and polyder to plot the following polynomial and its % derivative. % $$p(x) = 5(x-1)^3 - 2(x-1)^2 + 4(x-1) - 10 %$$ %
3. Using the norm function to compute % the 2-norm $\parallel \cdot \parallel$, % show that % $$\parallel {\bf v} \parallel = % \sqrt{v_1^2 + v_2^2 + v_3^2}$$ % where ${\bf v} = (-1,2,5)$. You % can compute the 2-norm of a vector (i.e. the length of the vector) with % the command norm(v,2). %
4. Construct a cubic spline through the polynomial % $$p(x) = x(x^2-1)$$ % at the points $[-2,-1,0,1,2]$. Use % the function ppval to interpolate the spline % at a set of points on the interval $[-2,2]$. Use polyder % to compute the derivative of the polynomial. % You can start with this code fragment. %
% xdata = [-2,-1,0,1,2];
% ydata = xdata.*(xdata.^2 - 1);
% pp = spline(xdata,ydata);
% xs = linspace(-2,2,100);
% ys = ppval(pp,xs);
% 
%
% Be sure you understand the structure pp. % Plot the original polynomial and its derivative on top of your spline results % to see that your spline results are correct. %
5. Construct and plot cubic spline through data points given by %
% N = 8;
% xdata = linspace(-5,5,N);
% ydata = -2 + 4*rand(1,N);
% 
%
6. Using the the data points from the previous example, compute and % plot the tangents and normals to the curves. Use the formulas % $$% {\bf T}(x) = \frac{{\bf r}'(x)}{g(x)} %$$ % and % $$% {\bf N}(x) = \frac{g(x) {\bf r}''(x) - g'(x){\bf r}'(x)}{g(x)^3 \kappa(x)} %$$ % where % $$% g(x) = \parallel {\bf r}'(x) \parallel \quad \mbox{and} \quad % \kappa(x) = \frac{\parallel {\bf r}''(x) \times {\bf r}'(x) \parallel}{g(x)^3} %$$ % To compute the 2-norm $\parallel {\bf v} \parallel$, % use norm(v,2). To compute the cross product ${\bf u}\times {\bf v}$, use the command cross(u,v). %
% Plot the tangents and the normals. Your final plot should look like : %
% %
%
7. Use the dot function to show that the normals and the % tangent vectors you created are perpendicular. That is, show that % $${\bf T} \cdot {\bf N} = 0$$. %
8. Show that the binormal vector, computed as % $${\bf B} = {\bf T} \times {\bf N}$$ % has zeros in the x- and y- components, is only non-zero in the % z-component. %
9. Show that % $${\bf T} \cdot {\bf B} = {\bf N} \cdot {\bf B} = 0$$ %
10. Show that % $$% \parallel {\bf B} \parallel = % \parallel {\bf N} \parallel = % \parallel {\bf T} \parallel = 1 %$$ %
% %
% %% % %