Chapter 2 M-Files

M-File Home

Script File Index

ShowV Illustrates InterpV and HornerV.
ShowN Illustrates InterpN and HornerN.
InterpEff Compares Vandermonde and Newton approaches.
RungeEg Examines accuracy of interpolating polynomial.
TableLookUp2D Illustrates SetUp and LinInterp2D.
ShowPolyTools Illustrates Matlab's polynomial tools.
Pallas Fits periodic data with CSInterp and CSEval.

Function File Index

InterpV Construction of Vandermonde interpolating polynomial.
HornerV Evaluates the Vandermonde interpolating polynomial.
InterpNRecur Recursive construction of the Newton interpolating polynomial.
InterpN Nonrecursive construction of the Newton interpolating polynomial.
InterpN2 Another nonrecursive construction of the Newton interpolating polynomial.
HornerN Evaluates the Newton interpolating polynomial.
SetUp Sets up matrix of f(x,y) evaluation.
LinInterp2D 2-Dimensional Linear Interpolation.
Humps2D A sample function of two variables.
CSInterp Fits periodic data with sines and cosines.
CSEval Evaluates sums of sines and cosines.

 


ShowV

% Script File: ShowV 
% Plots 4 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:4
   x = 2*pi*sort(rand(4,1));
   y = sin(x);
   a = InterpV(x,y);
   pVal = HornerV(a,x0);
   subplot(2,2,eg)
   plot(x0,y0,x0,pVal,'--',x,y,'*')
   axis([0 2*pi -2 2])
end

ShowN

% Script File: ShowN
% Random cubic interpolants of sin(x) on [0,2pi].
% Uses the Newton representation.

close all
x0 = linspace(0,2*pi,100)';
y0 = sin(x0);
for eg=1:20
   x = 2*pi*sort(rand(4,1));
   y = sin(x);
   c = InterpN(x,y);
   pvals = HornerN(c,x,x0);
   plot(x0,y0,x0,pvals,x,y,'*')
   axis([0 2*pi -2 2])
   pause(1)
end

InterpEff

% Script File: InterpEff
% Compares the Vandermonde and Newton Approaches

clc
disp('Flop Counts:')
disp(' ')
disp('  n   InterpV   InterpN     InterpNRecur')
disp('--------------------------------------')
for n = [4 8 16]
   x = linspace(0,1,n)'; y = sin(2*pi*x);
   flops(0); a = InterpV(x,y);      f1 = flops;
   flops(0); c = InterpN(x,y);      f2 = flops;
   flops(0); c = InterpNRecur(x,y); f3 = flops;
   disp(sprintf('%3.0f %7.0f    %7.0f     %7.0f',n,f1,f2,f3));
end

RungeEg

% Script File: RungeEg
% For n = 7:2:15, the equal-spacing interpolants of f(x) = 1/(1+25x^2) on [-1,1]'
% are of plotted.

close all
x = linspace(-1,1,100)'; 
y = ones(100,1)./(1 + 25*x.^2);
for n=7:2:15
   figure
   xEqual = linspace(-1,1,n)'; 
   yEqual = ones(n,1)./(1+25*xEqual.^2);
   cEqual=InterpN(xEqual,yEqual);
   pValsEqual = HornerN(cEqual,xEqual,x);
   plot(x,y,x,pValsEqual,xEqual,yEqual,'*')
   title(sprintf('Equal Spacing (n = %2.0f)',n))
end

TableLookUp2D

% Script File: TableLookUp2D
% Illustrates SetUp and LinearInterp2D.

close all
clc
a = 0; b = 5; n = 11;
c = 0; d = 3; m = 7;
fA = SetUp('Humps2D',a,b,n,c,d,m);
contour(fA)
x = input(sprintf('Enter x (%3.1f < = x < ="%3.1f" ):',a,b)); 
y = input(sprintf('Enter" y (%3.1f < ="y" < ="%3.1f" ):',c,d)); 
z= LinInterp2D(x,y,a,b,c,d,fA);
disp(sprintf('f(x,y)="%20.16f\n',z))" 

ShowPolyTools

% Script File: ShowPolyTools

close all

x = linspace(0,4*pi);
y = exp(-.2*x).*sin(x);

x0 = [.5 4  8 12];
y0 = exp(-.2*x0).*sin(x0);

p = polyfit(x0,y0,length(x0)-1);
pvals = polyval(p,x);

plot(x,y,x,pvals,x0,y0,'o')
axis([0 4*pi -1 1])
set(gca,'XTick',[])
set(gca,'YTick',[])
title('Interpolating {\ite}^{-.2{\itx}}sin({\itx}) with a Cubic Polynomial')

Pallas

% Script File: Pallas
% Plots the trigonometric interpolant of the Gauss Pallas data.

A = linspace(0,360,13)';
D = [ 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804 408]';

Avals = linspace(0,360,200)'; 
F = CSInterp(D(1:12));
Fvals = CSeval(F,360,Avals);

plot(Avals,Fvals,A,D,'o')
axis([-10 370 -200 1700])
set(gca,'xTick',linspace(0,360,13))
xlabel('Ascension (Degrees)')
ylabel('Declination (minutes)')

InterpV

  function a = InterpV(x,y)
% a = InverpV(x,y) 
% This computes the Vandermonde polynomial interpolant where
% x is a column n-vector with distinct components and y is a 
% column n-vector.
%
% a is a column n-vector with the property that if
%
%         p(x) = a(1) + a(2)x + ... a(n)x^(n-1) 
% then
%         p(x(i)) = y(i), i=1:n

n = length(x);
V = ones(n,n);
for j=2:n
   % Set up column j.
   V(:,j) = x.*V(:,j-1);
end
a = V\y;

HornerV

  function pVal = HornerV(a,z)
% pVal = HornerV(a,z)
% evaluates the Vandermonde interpolant on z where
% a is an n-vector and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
%          p(x) = a(1) + .. +a(n)x^(n-1)
% then
%          pVal(i) = p(z(i))  ,  i=1:m.

n = length(a); 
m = length(z);
pVal = a(n)*ones(size(z));
for k=n-1:-1:1
   pVal = z.*pVal + a(k);
end

InterpNRecur

  function c = InterpNRecur(x,y)
% c = InterpNRecur(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector. c is  a column n-vector with the property that if
%
%      p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
%      p(x(i)) = y(i), i=1:n.

n = length(x);
c = zeros(n,1);
c(1) = y(1);
if n > 1
   c(2:n) = InterpNRecur(x(2:n),(y(2:n)-y(1))./(x(2:n)-x(1)));
end

InterpN

  function c = InterpN(x,y)
% c = InterpN(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector. c is  a column n-vector with the property that if
%
%      p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
%      p(x(i)) = y(i), i=1:n.

n = length(x);
for k = 1:n-1
   y(k+1:n) = (y(k+1:n)-y(k)) ./ (x(k+1:n) - x(k));
end
c = y;

InterpN2

  function c = InterpN2(x,y)
% c = InterpN2(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector. c is  a column n-vector with the property that if
%
%      p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
%      p(x(i)) = y(i), i=1:n.

n = length(x);
for k = 1:n-1
   y(k+1:n) = (y(k+1:n)-y(k:n-1)) ./ (x(k+1:n) - x(1:n-k));  
end
c = y;

HornerN

  function pVal = HornerN(c,x,z)
% pVal = HornerN(c,x,z)
% Evaluates the Newton interpolant on z where
% c and x are n-vectors and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
%         p(x) = c(1) +  c(2)(x-x(1))+ ... + c(n)(x-x(1))...(x-x(n-1))
% then 
%         pval(i) = p(z(i)) , i=1:m.

n = length(c); 
pVal = c(n)*ones(size(z));
for k=n-1:-1:1
   pVal = (z-x(k)).*pVal + c(k);
end

SetUp

  function fA = SetUp(fname,a,b,n,c,d,m)
% fA = SetUp(fname,a,b,n,c,d,m)
% Sets up a matrix of f(x,y) evaluations.
% fname is a string that names a function of the form f(x,y).
% a, b, c, and d  are scalars that satisfy a<=b and c<="d." % n and m are integers>=2.
% fA is an n-by-m matrix with the property that
%
%     A(i,j) = f(a+(i-1)(b-a)/(n-1),c+(j-1)(d-c)/(m-1)) , i=1:n, j=1:m

x = linspace(a,b,n);
y = linspace(c,d,m);
fA = zeros(n,m);
for i=1:n
   for j=1:m
      fA(i,j) = feval(fname,x(i),y(j));
   end
end

LinInterp2D

  function z = LinInterp2D(xc,yc,a,b,c,d,fA)
% z = LinInterp2D(xc,yc,a,b,n,c,d,m,fA) 
% Linear interpolation on a grid of f(x,y) evaluations.
% xc, yc, a, b, c, and d  are scalars that satisfy a<=xc<=b and c<=yc<=d. 
% fA is an n-by-m matrix with the property that 
% 
% A(i,j)=f(a+(i-1)(b-a)/(n-1),c+(j-1)(d-c)/(m-1)) , i=1:n, j=1:m 
% 
% z is a linearly interpolated value of f(xc,yc). 
[n,m]=size(fA); 
% xc=a+(i-1+dx)*hx 0<=dx<=1 
hx=(b-a)/(n-1); 
i=max([1" ceil((xc-a)/hx)]); 
dx=(xc (a+(i-1)*hx))/hx;
 
% yc=c+(j-1+dy)*hy 0<=dy<=1 
hy=(d-c)/(m-1); 
j=max([1 ceil((yc-c)/hy)]); 
dy=(yc (c+(j-1)*hy))/hy; 
z=(1-dy)*((1-dx)*fA(i,j)+dx*fA(i+1,j)) + dy*((1-dx)*fA(i,j+1)+dx*fA(i+1,j+1)); 

Humps2D

  function z = Humps2D(x,y);
z = 1/( (.2*(x-3)^2 +.1) + (.3*(y-1)^2 + .1) );

CSInterp

  function F = CSInterp(f)
% F = CSInterp(f)
% f is a column n vector and n = 2m.
% F.a is a column m+1 vector and F.b is a column m-1 vector so that if 
% tau = (0:n-1)'*pi/m, then
%
%         f = F.a(1)*cos(0*tau) +...+ F.a(m+1)*cos(m*tau) + 
%             F.b(1)*sin(tau)   +...+ F.b(m-1)*sin((m-1)*tau)

n = length(f); 
m = n/2;
tau = (pi/m)*(0:n-1)';
P = [];
for j=0:m,   P = [P cos(j*tau)]; end
for j=1:m-1, P = [P sin(j*tau)]; end
y = P\f;
F = struct('a',y(1:m+1),'b',y(m+2:n));

CSEval

  function Fvals = CSEval(F,T,tvals)
% Fvals = CSEval(F,T,tvals)
% F.a is a length m+1 column vector, F.b is a length m-1 column vector, 
% T is a positive scalar, and tvals is a column vector.
% If 
%   F(t)  = F.a(1) + F.a(2)*cos((2*pi/T)*t) +...+ F.a(m+1)*cos((2*m*pi/T)*t) + 
%                    F.b(1)*sin((2*pi/T)*t) +...+ F.b(m-1)*sin((2*m*pi/T)*t)
%
% then Fvals = F(tvals).    

Fvals = zeros(length(tvals),1);
tau = (2*pi/T)*tvals;
for j=0:length(F.a)-1, Fvals = Fvals + F.a(j+1)*cos(j*tau); end
for j=1:length(F.b),   Fvals = Fvals + F.b(j)*sin(j*tau); end