Chapter 2: Problem Solutions
| 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 |
% 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;
% 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
% 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)
% 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));
z
% 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));
% 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
% 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)
z
z
% 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
% 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
% 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
% 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))
% 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
% 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);
% 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)]);
% Adjust for edge effects
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
% 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)
z
z