Lecture Handouts and Examples

| CS 322 Home | 

Lecture 10 Wednesday, February 25

 

% Script Lec10
% Looks at second derivative continuity of pwcubic hermite interpolant
bigscreen
close all
for s2 = [-10  -5  -1  -5/6  1   5  10 100]
    figure
    ShowSmallSpline(s2)
    pause(2)
end

  function ShowSmallSpline(s2)
% Displays the piecewise cubic hermite interpolant of the data
%
%                      (0,1,2), (1,3,s2), (2,-1,1)
%
% On [0,1] we have the cubic q1(z) = 1 + 2z +(s2-2)z^2(z-1) which satisfies
%
%         q1(0) = 1, q1'(0) = 2, q1(1) = 3, q1'(3) = s2
%
% It has second derivative
%
%         q1"(z) = 2(s2-2) + 6(s2-2)z


% The cubic q2(z) = 3 + s2(z-1) + (-4 - s2)(z-1)^2  + (9+s2)(z-1)^2(z-2)
%
%         q2(1) = 3, q1'(1) = s2,  q2(2)=-1, q2'(2) = 1
%
% It has second derivative
%
%         q2"(z) = 2(-4-s2) + (9+s2)(4(z-1) + 2(z-2))
%
% The equation q1"(1) = q2"(1) says
%
%               2(s2-2) +6(s2-2)  =  2(-4-s2) + (9+s2)(-2)
%
% i.e., s2 = 5/6. When this is the case we have a cubic spline interpolant.

close all

% Evaluate the two local cubics...
zL = linspace(0,1);
q1Vals = 1 + 2*zL + (s2-2)*(zL.^2).*(zL-1);
zR = linspace(1,2);
q2Vals = 3 + s2*(zR-1) + (-4-s2)*(zR-1).^2  + (9+s2)*((zR-1).^2 ).*(zR-2);

% Display
subplot(2,1,1)
plot(zL,q1Vals,zR,q2Vals,[0 1 2],[1 3 -1],'o')
legend('q1(z)','q2(z)')
title(sprintf('An n=2 piecewise cubic hermite interpolant with s2 = %5.2f',s2),'fontsize',18)

% Display their second derivatives
subplot(2,1,2)
plot(zL,2*(s2-2) + 6*(s2-2)*zL, zR,2*(-4-s2) + (9+s2)*(4*(zR-1)+2*(zR-2)))
legend('q1"(z)','q2"(z)')
title('It''s second derivative...','fontsize',18)

Lecture 4. Wednesday, February 4

Read about InterpV and HornerV in Chapter 2 of SCMV.

% Lec4A
% Plots random cubic interpolants of sin(x) on [0,2pi].
% Uses the Vandermonde method.

close all
x0 = linspace(0,2*pi,100)';
y0 = sin(x0);
for eg=1:30
   x = 2*pi*sort(rand(4,1));
   y = sin(x);
   a    = InterpV(x,y);
   pVal = HornerV(a,x0);
   plot(x0,y0,x0,pVal,'--',x,y,'r*')
   axis([0 2*pi -2 2])
   title('A Random Cubic Interpolant of sin(x) on [0,2\pi]' )
   legend('sin(x)','Cubic Interpolant')
   grid on
   shg
   pause
end
clc
grid off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Lec4B
% Equally spaced interpolants of sin(x) on [0,2pi].
% Uses the Vandermonde method.

close all
x0 = linspace(0,2*pi,100)';
y0 = sin(x0);
for deg=1:20
   x = linspace(0,2*pi,deg+1)';
   y = sin(x);
   a    = InterpV(x,y);
   pVal = HornerV(a,x0);
   subplot(2,1,1)
   plot(x0,y0,x0,pVal,'--',x,y,'r*')
   axis([0 2*pi -2 2])
   s = sprintf('deg %1d Interpolant of sin(x)',deg);
   title(s)
   legend('sin(x)','Polynomial Interpolant')
   subplot(2,1,2)
   absErr = abs(pVal-y0);
   plot(x0,absErr)
   axis([0 2*pi -inf inf])
   s = sprintf('Approx maximum absolute error = %5.1e',max(absErr));
   title(s,'fontsize',24)
   shg
   pause
end
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Lec4c
% Equally spaced interpolants of 1/(1+25x^2) on [-1 1].
% Uses the Vandermonde method.

close all
f = inline('1./(1+25*x.^2)');
x0 = linspace(-1,1,100)';
y0 = feval(f,x0);
for deg=1:20
   x = linspace(-1,1,deg+1)';
   y = feval(f,x);
   a    = InterpV(x,y);
   pVal = HornerV(a,x0);
   subplot(2,1,1)
   plot(x0,y0,x0,pVal,'--',x,y,'r*')
   s = sprintf('deg %1d Interpolant of 1/(1+25x^{2})',deg);
   title(s,'fontsize',24)
   legend('1/(1+25x^{2})','Polynomial Interpolant',0)
   subplot(2,1,2)
   absErr = abs(pVal-y0);
   plot(x0,absErr)
   s = sprintf('Approx maximum absolute error = %5.1e',max(absErr));
   title(s,'fontsize',24)
   shg
   pause
end
clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Script Lec4d
  close all
  
  ShowInterp('sin',0,2*pi,5)
  
  figure
  g = inline('sin(2*x)');
  ShowInterp(g,0,2*pi,5)
  
  figure
  ShowInterp(@sin,0,2*pi,5)
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function ShowInterp(fName,a,b,deg)
% fName is a string that names a function f defined on
% [a,b]. Assume a < b
% deg is a positive integer
% Displays the equal-spacing, degree d polynomial interpolant 
% of f on [a,b].

% Display the function over  the interval...
  x = linspace(a,b)';
  y = feval(fName,x);
  plot(x,y);
  axis([a b -inf inf])

% Display the points of interpolation
  hold on
  x0 = linspace(a,b,deg+1)';
  y0 = feval(fName,x0);
  plot(x0,y0,'or')
  
% Generate the interpolant
  a = InterpV(x0,y0);
  
% Evaluate the interpolant across the interval and display.
  pVals = HornerV(a,x);
  plot(x,pVals,'r')
  hold off
  shg

Lecture 3. Monday, February 2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script Lec3A

x = linspace(-7,7);
y = MyExp(x); 
shg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function y = MyExp(x)
% x is a vector
% y has the same length and orientation with 
% y(i) = exp(x(i)), i = 1:length(x)

y = ones(size(x));    % y = vector of all 1's, same size & shape as x
termIndex = 1;
termNext  = x;
yNext = y + termNext;
while(any(y~=yNext))

    % Iterate as long as the addition of the next term makes a difference
    y = yNext;
    termIndex = termIndex + 1;
    termNext = (termNext.*x)/termIndex;
    yNext = yNext + termNext;
    
    % Display the relative error
    relErr = abs(yNext - exp(x))./exp(x);
    % Add eps to relErr so don't evaluate log at zero.
    semilogy(x,relErr+eps,[min(x) max(x)] ,[10^-15 10^-15],'r')
    title(sprintf('Number of terms = %1d',termIndex))
    ylabel('Relative Error')
    xlabel('x')
    shg
    pause
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Script Lec3B
% Explores running sums in Taylor expansion for exp(x)

close all
x = -7;
nTerms = 30;

sum  = zeros(1,nTerms);
term = zeros(1,nTerms);

sum(1)  = 1 + x;
term(1) = x;

for k=2:nTerms
    term(k) = x*term(k-1)/k;
    sum(k) = sum(k-1) + term(k);
end
semilogy(1:nTerms,abs(sum),1:nTerms,abs(term))
title(sprintf('exp(%5.1f) via taylor series',x))
xlabel('k')
legend('1 + term(1) + ... + term(k)',' | term(k) |')
shg