% 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 % %

%

%
%%
%
% -
%
- Use the functions polyval and polyder to plot the following polynomial and its % derivative. % $$p(x) = 5x^3 - 2x^2 + 4x - 10$$ %
- 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 % $$ %
- 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). %
- 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. % - Construct and plot cubic spline through data points given by
%
% N = 8; % xdata = linspace(-5,5,N); % ydata = -2 + 4*rand(1,N); %

% - 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 : %% %

% - 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$$. %
- 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. % - Show that % $${\bf T} \cdot {\bf B} = {\bf N} \cdot {\bf B} = 0$$ %
- Show that % $$ % \parallel {\bf B} \parallel = % \parallel {\bf N} \parallel = % \parallel {\bf T} \parallel = 1 % $$ %

Compare your answers with the solutions.

%