Chapter 2: Problem Solutions

Problems Home

 P2.1.1 P2.2.1 P2.3.1 P2.4.1 P2.1.2 P2.2.2 P2.3.2 P2.4.2 P2.1.3 P2.2.3 P2.3.3 P2.4.3 P2.1.4 P2.2.4 P2.4.4 P2.1.5 P2.2.5 P2.4.5 P2.4.6 P2.4.7

P2.1.1

```% Problem P2_1_1
%
% Modified Vandermonde interpolation of census data
% Census Data from 1820 through 1980.

x = (1820:10:1980)';
y = [  9638453; ...
12866020; ...
17068953; ...
23191876; ...
31443321; ...
38558371; ...
50189209; ...

62979766; ...
75994575; ...
91972266; ...
105710620; ...
122775046; ...
131669275; ...
150697361; ...
179323175; ...
203235298; ...
226547082];
clc

disp('Predict  1980 population based on census data')
disp('for years k:10:1970 where k is the starting year.')
disp(' ')
yr = input('Enter starting year: ');
z = linspace(yr,1980)';
k = 16 - (1970-yr)/10;
close all
format long
figure
a = MyInterpV(x(k:16),y(k:16));
pval = HornerV(a,z);
plot(z,pval,x(k:17),y(k:17),'o')
title('Prediction of 1980 population using basis x^j')

figure
u = yr;
a = MyInterpV(x(k:16),y(k:16),u);
pval = HornerV(a,z-u);
plot(z,pval,x(k:17),y(k:17),'o')
title(sprintf('Prediction of 1980 population using basis (x-%4.0f)^j',u))

figure
u = (yr+1970)/2;
a = MyInterpV(x(k:16),y(k:16),u);
pval = HornerV(a,z-u);
plot(z,pval,x(k:17),y(k:17),'o')
title(sprintf('Prediction of 1980 population using basis (x-%4.0f)^j',u))
figure
u = (yr+1970)/2;
v = (1970-yr)/2;
a = MyInterpV(x(k:16),y(k:16),u,v);
pval = HornerV(a,(z-u)/v);
plot(z,pval,x(k:17),y(k:17),'o')
title(sprintf('Prediction of 1980 population using basis ((x-%4.0f)/%2.0f)^j',u,v))

function a = MyInterpV(x,y,u,v)
%
%   x is a column n-vector with distinct components.
%   y is a column n-vector.
%   u is a scalar (optional).
%   v is a scalar (optional).
%
%   a is a column n-vector with the property that
%   if p(x) = a(1) + a(2)xbar + ... a(n)xbar^(n-1) with
%       xbar = (x-u)/v then p(x(i)) = y(i), i=1:n.
%
%   A call of the form MyInterpV(x,y,u) assumes v=1.
%   A reference of the form MyInterpV(x,y) assumes u=0 and v=1.

n = length(x);
V1 = ones(n,n);
if nargin <= 3
v = 1;
end

if nargin == 2
u = 0;
end

xbar = (x-u)/v;
for j=2:n
% Set up column j.
V1(:,j) = xbar.*V1(:,j-1);
end
a = V1\y;
```

P2.1.2

```% Problem P2_1_2
%
% Polynomial evaluation for odd and even polynomials.

close all

z  = linspace(-pi,pi)';
c = [1 1 2 6 24 120 720 5040 40320 322880  3228800];
a = [1 -1/6  1/120  -1/5040 1/322880];
y = HornerVgen(a,z,'odd');
subplot(2,1,1)
plot(z,y)
title('First 5 terms of Taylor series for sin(x)')

a = [1 -1/2  1/24  -1/720 1/40320 ]';
y = HornerVgen(a,z,'even');
subplot(2,1,2)
plot(z,y)
title('First 5 terms of Taylor series for cos(x)')

function pval = HornerVgen(a,z,s)
%
% a is a column n-vector.
% z is a column m-vector.
% s is a string equal to 'odd' or 'even'.
%
% pval is a column m-vector with the property that if
% pval(i) = p(z(i)) for i=1:m where
%
%     p(x) = a(1)x + a(2)x^3 + a(3)x^5 +...+ a(n)x^(2n-1) if s = 'odd'
%     p(x) = a(1)  + a(2)x^2 + a(3)x^4 +...+ a(n)x^(2n-2) if s = 'even'

n = length(a);
m = length(z);

z2 = z.*z;

pval = a(n)*ones(m,1);
for k=n-1:-1:1
pval = z2.*pval + a(k);
end
if (length(s) == 3)
if s == 'odd'
pval = z.*pval;
end
end
```

P2.1.3

```% Problem P2_1_3
%
% Odd and even polynomials.

clc
z = linspace(0,1,10)';
a = rand(9,1);

% Part A.
%
% Write p(x) = (pe(x) + po(x))/2 where pe(x) = (p(x)+p(-x)) and
% po(x) = (p(x)-p(-x)). Note that pe and po are even and odd.
% See P2.1.2. It follows that
%            p(x)/p(-x) = (pe(x)+po(x))/(pe(x)-po(x))
% The cost of evaluating pe AND po is roughly equal to the cost
% of evaluating p. Thus, we can evaluate the given rational with
% little more work than required to evaluate p.

z2 = z.^2;
n = length(a);
pe =    HornerV(a(1:2:n),z2);
po  = z.*HornerV(a(2:2:n),z2);
partA = (pe + po)./(pe - po);
CheckA = (HornerV(a,z)./HornerV(a,-z)) - partA

% Part B.
%
% Compute the even polynomial p(x) + p(-x). See above.

partB = 2*pe;
CheckB = (HornerV(a,z) + HornerV(a,-z)) - partB

% Part C.
%
% Check the derivative calculation via divided differences:
partC = HornerV(a(2:n).*(1:(n-1))',z);
CheckC = ((HornerV(a,z+10^(-8))-HornerV(a,z))/10^(-8)) - partC

partD = HornerV([0 ; a./(1:n)'],1)
partE = HornerV([0 ; a./(1:n)'],z) - HornerV([0 ; a./(1:n)'],-z)
```

P2.1.4

```% Problem P2_2_1
%
% Conversion from the Newton to the 1,x,x^2,.. repesentation

clc

z = linspace(0,2*pi)';
n = 7;
x = linspace(0,2*pi,n)';
y = cos(x);
c = InterpN(x,y);
a1 = N2V(c,x(1:n-1));
a2 = InterpV(x,y);
home
disp(sprintf('Coefficients for interpolant of cos(z) at z = linspace(0,2pi,%1.0f)',n))
disp(' ')
disp('       via  N2V                via InterpV')
disp('------------------------------------------------------')
for i=1:n
disp(sprintf(' %20.16f    %20.16f',a1(i),a2(i)))
end

function a = N2V(c,x)
%
% c is a column n-vector
% x is a column (n-1)-vector
%
% a is a column n-vector so that if
%       p(z) = c(1) + c(2)(z-x(1)) + ... + c(n)(z-x(1))...(z-x(n-1))
%

z = linspace(min(x),max(x),length(c))';

% Evaluate p(z) at the n points z(1),...,z(n) and construct the
% Vandermonde interpolant of (z(i),p(z(i)), i=1:n>

a = InterpV(z,HornerN(c,x,z));
```

P2.1.5

```z
```

P2.2.1

```% Problem P2_2_1
%
% Conversion from the Newton to the 1,x,x^2,.. repesentation

clc

z = linspace(0,2*pi)';
n = 7;
x = linspace(0,2*pi,n)';
y = cos(x);
c = InterpN(x,y);
a1 = N2V(c,x(1:n-1));
a2 = InterpV(x,y);
home
disp(sprintf('Coefficients for interpolant of cos(z) at z = linspace(0,2pi,%1.0f)',n))
disp(' ')
disp('       via  N2V                via InterpV')
disp('------------------------------------------------------')
for i=1:n
disp(sprintf(' %20.16f    %20.16f',a1(i),a2(i)))
end

function a = N2V(c,x)
%
% c is a column n-vector
% x is a column (n-1)-vector
%
% a is a column n-vector so that if
%       p(z) = c(1) + c(2)(z-x(1)) + ... + c(n)(z-x(1))...(z-x(n-1))
%

z = linspace(min(x),max(x),length(c))';

% Evaluate p(z) at the n points z(1),...,z(n) and construct the
% Vandermonde interpolant of (z(i),p(z(i)), i=1:n>

a = InterpV(z,HornerN(c,x,z));
```

P2.2.2

```% Problem P2_2_2
%
% Illustrate recursive evaluation of interpolating polynomial.

close all
z = linspace(0,2*pi)';
x = linspace(0,2*pi,11)';
y = sin(x);
pval = RecurEval(x,y,z);
plot(z,pval,x,y,'o')
title('Interpolant of sin(x) at x = linspace(0,2pi,11)')

function pval = RecurEval(x,y,z)
% x is a column n-vector with distinct entries
% y is a column n-vector
% z is a column m-vector
%
% pval is a column m-vector with  the property that pval(i) = p(x(i))
%    where p(x) is the degree n-1 interpolant of (x(i),y(i)), i=1:n.

n = length(x);
if n==1
pval(1) = y(1);
else
pvalL = RecurEval(x(1:n-1),y(1:n-1),z);
pvalR = RecurEval(x(2:n),y(2:n),z);
pval = ((z-x(n)).*pvalL - (z-x(1)).*pvalR)/(x(1)-x(n));
end
```

P2.2.3

```% Problem P2_2_3
%
% Interactive plot of a function and a chosen interpolant.

clc

s = input('Enter the name of a function in quotes: ');
n = input('Enter the number of interpolation points: ');
L = input('Enter the left endpoint of the interval: ');
R = input('Enter the right endpoint of the interval: ');
x = linspace(L,R,n)';
y = feval(s,x);
z = linspace(L,R)';
c = InterpN(x,y);
pval = HornerN(c,x,z);
fval = feval(s,z);
plot(z,fval,z,pval,x,y,'o')
t = ['Plot of ' s ' and a ' sprintf('%1.0f',n) '-point interpolant'];
title(t)
```

P2.2.4

```z
```

P2.2.5

```z
```

P2.3.1

```% Problem P2_3_1
%
% Flop comparison of HornerN and HornerV

F = zeros(9,4);
a = linspace(0,1,10)';
c = linspace(0,1,10)';
x = linspace(0,1,10)';
z = linspace(0,1,200)';
i=1;
for n=2:10
j=1;
for m=50:50:200
flops(0)
pval = HornerV(a(1:n),z(1:m));
fV = flops;
flops(0)
pval = HornerN(c(1:n),x(1:n),z(1:m));
fN = flops;
F(i,j) = fN/fV;
j=j+1;
end
i=i+1;
end

clc
disp('(Flops HornerN)/(Flops HornerV)')
disp(' ')
disp('  n      m=50    m=100    m=150    m=200')
disp('--------------------------------------------')
for i=1:9
s1 =sprintf('  %5.3f  ',F(i,:));
disp([sprintf('%3.0f   ',i+1) s1])
end
```

P2.3.2

```% Problem 2_3_2
%
% Explore the function e(s)

close all
for n=2:10
z = linspace(0,n-1)';
e = z;
for k=1:n-1
e = e.*(z-k)/(k+1);
end
plot(z,e)
title(sprintf('n=%2.0f  max|e(s)| = %4.2f',n,max(abs(e)) ))
pause(1)
end
```

P2.3.3

```% Problem P2_3_3
%
% Using the interpolation polynomial error bound.
% f(x) = exp(ax). For the n-point interpolant, the n-th
% derivative is important and it equals a^n exp(ax).
% It follows that it is no bigger than
%    derBound(n) =  (abs(a)^n)*max([exp(a*L) exp(a*R)]);
%
% From (2.2) it follows that the error in the interpolant is bounded by
%
%     e(n) = (derBound(n)/4n)( (R-L)/(n-1) )^n

clc
disp('   R    a : delta=10^-3   delta=10^-4   delta=10^-5   delta=10^-6')
disp('--------------------------------------------------------------------')
for a=[-5 -1 1 5]
L = 0;
for R = [1 2 4 8];
n1 = nbest(L,R,a,.001);
n2 = nbest(L,R,a,.0001);
n3 = nbest(L,R,a,.00001);
n4 = nbest(L,R,a,.00001);
disp(sprintf('%4.0f %4.0f :    %4.0f         %4.0f          %4.0f          %4.0f',a,R,n1,n2,n3,n4))
end
end

function n = nbest(L,R,a,delta)
% Smallest n so that the n-point equal spacing interpolant of
% exp(ax) on [L,R} has error less than delta.

n = 2;
while ((abs(a)^n)*max([exp(a*L) exp(a*R)])/(4*n))*((R-L)/(n-1))^n  > delta
n=n+1;
end
```

P2.4.1

```% Problem P2_4_1
%
% Compare InterpN and InterpN2

n=21;
format long
x = linspace(-1,1,n)';
y = 1./(1+25*x.^2);
c1 = InterpN(x,y);
c2 = InterpN2(x,y);
clc
home
disp('abs(InterpN(x,y) - InterpN2(x,y)) = ')
disp(' ')
disp(abs(c1 - c2))
```

P2.4.2

```% Problem P2_4_2
%
% Compare InterpN and InterpNEqual

clc

n = 10;
L = 0;
R = 2*pi;
x = linspace(L,R,n)';
y = sin(x);

c = InterpN(x,y);
[cE xE yE] = InterpNEqual('sin',L,R,n);

disp('Compare InterpNEqual and InterpN.')
disp(' ')
disp('abs(c-cE) = ')
disp(' ')
disp(abs(c-cE))

function [c,x,y] = InterpNEqual(fname,L,R,n)
% fname is a string that names a function of the form f(x)
% L < R and n is a positive integer.
%
% x = linspace(L,R,n)'
% y = f(x)
% c = the coefficients of the Newton interpolant of f at x.

x = linspace(L,R,n)';
h = (R-L)/(n-1);
y = feval(fname,x);
c = y;
for k=2:n
c(k:n) = diff(c(k-1:n))/(h*(k-1));
end
```

P2.4.3

```% Problem P2_4_3
%
% Inverse Quadratic Interpolation. Find zero of sin function
% in the vicinity of PI.

clc
x = [3 ; 3.5 ; 4];
y = sin(x);
root = InverseQ(x,y);
disp(sprintf('Root = %8.6f', root))

function root = InverseQ(x,y)
% x and y are column 3-vectors with x(1) < x(2) < x(3) and y(1)y(3)<0
% root = q(0) where q is the quadratic interpolant
%       of (y(1),x(1)), (y(2),x(2)), (y(3),x(3))

a = InterpV(y,x);
root = HornerV(a,0);
```

P2.4.4

```% Script P2_4_4
%
% Illustrates SetUp and CubicInterp2D.

a = 0; b = 5; n = 16;
c = 0; d = 3; m = 10;
fA = SetUp('Humps2D',a,b,n,c,d,m);
E = zeros(n-1,m-1);
hx = (b-a)/(n-1);
hy = (d-c)/(m-1);
clc

disp('E(i,j) is the error in 2D cubic interpolant of humps2D(xc,yc)')
disp('at the center of the rectangle {(u,v)|x(i)<=u<=x(i+1), y(j)<=v<=y(j+1)}')
disp('where x = linspace(0,5,16) and y = linspace(0,3,10).')
disp(' ')

disp('E = ')
disp(' ')
for i=1:n-1
for j=1:m-1
xc = a+(i-.5)*hx;
yc = c+(j-.5)*hy;
exact = Humps2D(xc,yc);
approx = CubicInterp2D(xc,yc,a,b,n,c,d,m,fA);
E(i,j) = abs(approx-exact);
end

disp(sprintf('%6.4f  ',E(i,:)))
end

function z = CubicInterp2D(x,y,a,b,n,c,d,m,fA)
%
%  n,m:      integers >=2
%  x,a,b:    scalars that satisfy a<=x<\=b
%  y,c,d:    scalars that satisfy c<=y<=d
%  fA:       an n-by-m matrix with the property that
%              fA(i,j) = f(a+(i-1)(b-a)/(n-1),c+(j-1)(d-c)/(m-1))
%              for some function f(.,.).
% Post:
%    z:        cubically interpolated value for f(x,y).
%
% Determine where in the grid the point is located.

hx = (b-a)/(n-1);
i  = max([1 ceil((x-a)/hx)]);
hy = (d-c)/(m-1);
j  = max([1 ceil((y-c)/hy)]);
if (1 < i)&(i < n-1)
% not on an edge
ivals = (i-1):(i+2);
elseif 1==i
% on left edge
ivals = 1:4;
else
% on right edge
ivals = n-3:n;
end
if (1 < j)&(j < m-1)
% not on an edge
jvals = (j-1):(j+2);
elseif 1==j
% on bottom edge
jvals = 1:4;
else
% on top edge
jvals = m-3:m;
end
% Determine coordinates of the grid points.
xvals = c + (ivals-1)*hx;
yvals = c + (jvals-1)*hy;
% call Small2D with the appropriate 4x4 subgrid:

z = Small2D(x,y,xvals',yvals',fA(ivals,jvals));

function z = Small2D(x,y,xvals,yvals,fvals)
%
% x and y are scalars representing the point at which to interpolate f
% xvals and yvals are column 4-vectors with
%     xvals(1) < xvals(2) < xvals(3) < xvals(4)
%     yvals(1) < yvals(2) < yvals(3) &ly yvals(4)
% fvals is a 4x4 matrix of evaluations of f at the xvals,yvals coords.
%
% Let p1 be the cubic interpolant of (xvals(q),fvals(q,1)), q=1:4.
% Let p2 be the cubic interpolant of (xvals(q),fvals(q,2)), q=1:4.
% Let p3 be the cubic interpolant of (xvals(q),fvals(q,3)), q=1:4.
% Let p4 be the cubic interpolant of (xvals(q),fvals(q,4)), q=1:4.
%
% Let p be the cubic interpolant of the four points
%                    (yvals(1),p1(x))
%                    (yvals(2),p2(x))
%                    (yvals(3),p3(x))
%                    (yvals(4),p4(x))
%
%      z = p(y).

fy = zeros(4,1);
for k=1:4
c = InterpN(xvals,fvals(:,k));  % interpolate horizontally
fy(k) = HornerN(c,xvals,x);     % evaluate interpolants at x
end

c = InterpN(yvals,fy);   % interpolate vertically
z = HornerN(c,yvals,y);  % evaluate interpolant at y
```

P2.4.5

```% Problem P2_4_5
%
% Plotting the derivative of the polynomial interpolant.

n = 6;
x = linspace(0,2*pi,n);
y = sin(x);

z = linspace(0,2*pi);
yp = cos(z);
plot(z,yp,'.')

hold on
PlotDerPoly(x,y)
hold off

function PlotDerPoly(x,y)
% x and y are column n-vectors with x(1) < ... < x(n)
% Plots p'(z) where p is the polynomial interpolant of
%       (x(i),y(i)), i=1:n

n = length(x);
a = polyfit(x,y,n-1);

der_a = a(1:n-1).*((n-1):-1:1);
z = linspace(x(1),x(n));
der_pvals = polyval(der_a,z);
plot(z,der_pvals)
```

P2.4.6

```z
```

P2.4.7

```z
```
`  `