Back to tutorial index

Lab Solutions


Topics in this lab

Introduction

clear all;
close all;

Back to the top

Problem 1

% get coefficients for the polynomial
clf;
c = [5, -2, 4, -1];
xs = linspace(-5,5,100);
ys = polyval(c,xs);
p1 = plot(xs,ys,'k');
hold on;

c_deriv = polyder(c);
ys_deriv = polyval(c_deriv,xs);

p2 = plot(xs,ys_deriv,'r');
legend([p1 p2],{'p(x)','p''(x)'});
snapnow;

Back to the top

Problem 2

clear all;
clf;

c = [5, -2, 4, -1];
xs = linspace(-5,5,100);
ys = polyval(c,xs-1);
p1 = plot(xs,ys,'k');
hold on;

c_deriv = polyder(c);
ys_deriv = polyval(c_deriv,xs-1);

p2 = plot(xs,ys_deriv,'r');
legend([p1 p2],{'p(x)','p''(x)'});
snapnow;

Back to the top

Problem 3

clf;
clear all;
v = [-1, 2, 5];
format long;
disp(norm(v,2))
disp(sqrt(dot(v,v)));
disp(sqrt(v(1)^2 + v(2)^2 + v(3)^2));
   5.477225575051660

   5.477225575051661

   5.477225575051661

Back to the top

Problem 4

clf;
clear all;

xdata = [-2,-1,0,1,2];
ydata = xdata.*(xdata.^2 - 1);
pp = spline(xdata,ydata);
xs = linspace(-2,2,100);
ys = ppval(pp,xs);

p1 = plot(xs,ys,'b','linewidth',2);
hold on;
plot(xs,xs.*(xs.^2 - 1),'r','linewidth',1);

d2_coefs = zeros(pp.pieces,3);

for i = 1:pp.pieces
    d2_coefs(i,:) = polyder(pp.coefs(i,:));
end

dpp = mkpp(pp.breaks,d2_coefs);

ys_deriv = ppval(dpp,xs);

plot(xs,ys_deriv,'k','linewidth',2);
plot(xs,3*xs.^2 - 1,'r','linewidth',1);
snapnow;

Back to the top

Problem 5

clear all;
clf;

N = 8;
xdata = linspace(-5,5,N);
ydata = -2 + 4*rand(1,N);

pp = spline(xdata,ydata);
xs = linspace(-5,5,100);
ys = ppval(pp,xs);
plot(xs,ys);
hold on;
plot(xdata,ydata,'.','markersize',30);

Back to the top

Problem 6

clear all;
clf;

N = 10;
xdata = linspace(-5,5,N);
ydata = -2 + 4*rand(1,N);
pp = spline(xdata,ydata);
xs = linspace(-5,5,200);
ys = ppval(pp,xs);
plot(xs,ys,'r','linewidth',2);
hold on;

d_coefs = zeros(pp.pieces,3);
d2_coefs = zeros(pp.pieces,2);

for j = 1:pp.pieces,
    c = pp.coefs(j,:);
    d_coefs(j,:) = polyder(c);
    d2_coefs(j,:) = polyder(d_coefs(j,:));
end

% First and second derivatives of the spline.
dpp = mkpp(pp.breaks,d_coefs);
d2pp = mkpp(pp.breaks,d2_coefs);

% Compute the tangents and normals at the data points.  We will make
% r'(x) and r''(x) 3-component vectors so that the norm and cross product
% functions in matlab can be used. The third component will just be zero
% since the curve is restricted to the x-y plane.
fprime = ppval(dpp,xdata);
rprime = [ones(1,N); fprime; zeros(1,N)];
f2prime = ppval(d2pp,xdata);
r2prime = [zeros(1,N); f2prime; zeros(1,N)]; % Cross product requires

tangents = zeros(3,N);
normals = zeros(3,N);
% Loop over the values in xdata to compute the tangents and normals.
for j = 1:N
    g = norm(rprime(:,j),2);
    tangents(:,j) = rprime(:,j)/g;
    gprime = dot(rprime(:,j),r2prime(:,j))/g;
    kappa = norm(cross(r2prime(:,j),rprime(:,j)),2)/g^3;
    normals(:,j) = (g*r2prime(:,j)-gprime*rprime(:,j))/(g^3*kappa);
end

% Plot the original data
plot(xdata,ydata,'k.','markersize',30);

% Plot the tangents, use only the first two components since this is a
% 2D curve
quiver(xdata,ydata,tangents(1,:),tangents(2,:),0.25,'k-','LineWidth',1);

% Plot the normals, use only the first two components since this is a
% 2D curve
quiver(xdata,ydata,normals(1,:),normals(2,:),0.25,'b-','LineWidth',1);

title('Tangents and Normals (spline)','fontsize',18);
xlim([-5.5 5.5]);
ylim([-3 3]);
daspect([1 1 1]);
snapnow;

Notice how the normal always points towards the radius of curvature.

Back to the top

Problem 7

for i = 1:N-1,
    fprintf('T(i) dot N(i) = %12.4e\n',dot(tangents(:,i),normals(:,i)));
end
T(i) dot N(i) =  -2.7756e-17
T(i) dot N(i) =   4.9960e-16
T(i) dot N(i) =   1.6653e-16
T(i) dot N(i) =  -3.3307e-16
T(i) dot N(i) =   0.0000e+00
T(i) dot N(i) =   1.1102e-16
T(i) dot N(i) =   0.0000e+00
T(i) dot N(i) =   5.5511e-17
T(i) dot N(i) =  -1.1102e-16

Back to the top

Problem 8

binormals = zeros(3,N);
for i = 1:N-1,
    b = cross(tangents(:,i),normals(:,i));
    binormals(:,i) = b(:)';
    fprintf('%24.16f %24.16f %24.16f\n',b(1),b(2),b(3));

end
      0.0000000000000000       0.0000000000000000      -1.0000000000000000
      0.0000000000000000      -0.0000000000000000      -1.0000000000000000
      0.0000000000000000      -0.0000000000000000       1.0000000000000000
      0.0000000000000000       0.0000000000000000      -1.0000000000000002
     -0.0000000000000000       0.0000000000000000       1.0000000000000000
     -0.0000000000000000       0.0000000000000000       1.0000000000000000
      0.0000000000000000      -0.0000000000000000       1.0000000000000000
      0.0000000000000000      -0.0000000000000000      -1.0000000000000002
     -0.0000000000000000       0.0000000000000000       1.0000000000000000

Back to the top

Problem 9

for i = 1:N-1,
    fprintf('B(i) dot N(i) = %12.4e\n',dot(binormals(:,i),normals(:,i)));
end
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
B(i) dot N(i) =   0.0000e+00
fprintf('\n');
for i = 1:N-1,
    fprintf('B(i) dot T(i) = %12.4e\n',dot(binormals(:,i),tangents(:,i)));
end
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00
B(i) dot T(i) =   0.0000e+00

Back to the top

Problem 10

for i = 1:N-1,
    d1 = norm(binormals(:,i),2);
    d2 = norm(tangents(:,i),2);
    d3 = norm(normals(:,i),2);
    fprintf('%24.16f %24.16f %24.16f\n',d1,d2,d3);
end
      1.0000000000000000       1.0000000000000000       1.0000000000000000
      1.0000000000000000       0.9999999999999999       1.0000000000000002
      1.0000000000000000       1.0000000000000000       1.0000000000000002
      1.0000000000000002       1.0000000000000000       1.0000000000000004
      1.0000000000000000       1.0000000000000000       1.0000000000000000
      1.0000000000000000       1.0000000000000000       0.9999999999999999
      1.0000000000000000       1.0000000000000000       1.0000000000000000
      1.0000000000000002       1.0000000000000000       1.0000000000000000
      1.0000000000000000       1.0000000000000000       1.0000000000000000

Back to lab exercises.

Back to the top

Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_15_solns.m)

Powered by MathJax

Published with MATLAB® 8.4