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
```

```% 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
```

```% 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))"
```

```% 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')
```

```% 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
```

```  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));
```

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

```  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
```