Chapter 2 M-Files

M-File Home

Script File Index

Function File Index

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