Lecture Handouts and Examples

| CS 322 Home

 Lec 3 Feb 2 Lec 4 Feb 4 Lec 10 Feb 25

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