Prelim 1 Review Problem Solutions

% Problem 1
%
% The spacing between base-2 floating point numbers of the form m*2^e
% is 2^(-t)*2^e where t is the length of the mantissa in bits.
% If M and M+2 are floating point numbers and M+1 is not, then the
% spacing of floating point numbers between M and M+2 is >= 2.
% Suppose M = m*2^e. It follows that
%
%    M*2^{-t} = m*(2^e)*(2^{-t}) >= 2m
%
% and so
%
%            2^{t} <= M/m <= 2M.
%
% Thus, t <= log_{2}(M) - 1

M = 1;
k = 0;
while ((M < M+1) & ((M+1)<(M+2)))
   M = 2*M;
end
t = log2(M)-1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Problem 2 

n = 20;
A = zeros(n,n);
for i=1:n
   for j=1:n
      if abs(i-j)<=2
         A(i,j) = cos(abs(i-j)*pi/n);
      end
   end
end
x = randn(n,1);
error = max(abs(CosProd(x)-A*x))


function t = CosProd(x)
% x is a column n-vector and n > 3.
% y = A*x where A is the n-by-n matrix defined by
%
%           0                 if |i-j|>2
% A(i,j) =
%           cos(|i-j)|*pi/n)  if |i-j|<=2
   
n = length(x);
u = cos(pi/n);
v = cos(2*pi/n);
   
% The pattern...
   
%     t(1)   =                      x(1)   + ux(2)   + vx(3)
%     t(2)   =            ux(1)   + x(2)   + ux(3)   + vx(4)
%     t(3)   =    vx(1) + ux(2)   + x(3)   + ux(4)   + vx(5)
%
%     t(n-2) =  vx(n-4) + ux(n-3) + x(n-2) + ux(n-1) + vx(n)     
%     t(n-1) =  vx(n-3) + ux(n-2) + x(n-1) + ux(n) 
%     t(n)   =  vx(n-2) + ux(n-1) + x(n) 
   
% The solution
   
t = x;
t(1)     = t(1) + u*x(2) + v*x(3);
t(2)     = t(2) + u*x(1) + u*x(3) + v*x(4);
t(3:n-2) = t(3:n-2) + u*(x(2:n-3)+x(4:n-1)) + v*(x(1:n-4)+x(5:n));
t(n-1)   = t(n-1) + v*x(n-3) + u*x(n-2) + u*x(n);
t(n)     = t(n) + v*x(n-2) + u*x(n-1);
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
% Problem 3

n = 10;
x = linspace(0,pi,n);
y = exp(-x).*sin(x);
m = MaxJump(x,y)

  function m = MaxJump(x,y)
% x and y are column n-vectors and x(1) < x(2) <...< x(n).
% Let S be the cubic spline produced by spline(x,y).
% m is the maximum value of f(z) on the interval [x(2) , x(n-1)] where
% f(z) is the limit of |S'''(z+delta) - S'''(z-delta)| as delta goes to zero.

S = spline(x,y);

[breaks,coeff,L] = unmkpp(S)

% Compute the 3rd derivative of each local cubic.

Der3  = 6*coeff(:,1)             
m = max(abs(Der3(1:L-1)-Der3(2:L)));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Problem 4

a = [2;5;-1;2]
alfa = 3;
beta = 1;
c = Convert(a,alfa,beta);

z = linspace(0,5)';

zalfa = z-alfa;
pvals = ((a(4)*zalfa + a(3)).*zalfa + a(2)).*zalfa + a(1);

zbeta = z-beta;
qvals = ((c(4)*zbeta + c(3)).*zbeta + c(2)).*zbeta + c(1);

error = max(abs(qvals - pvals))


  function c = Convert(a,alfa,beta)
% a is a column 4-vector and alfa and beta are scalars.
% c is a column 4-vector so that the cubic polynomials
%
%       p(x) = a(1) + a(2)(x-alfa) + a(3)(x-alfa)^2 + a(4)(x-alfa)^3
%
%       q(x) = c(1) + c(2)(x-beta) + c(3)(x-beta)^2 + c(4)(x-beta)^3

% Note that if q interpolates p at 4 distinct points then q = p.

x = [1;2;3;4]
xalfa = x - alfa;
pvals = ((a(4)*xalfa + a(3)).*xalfa + a(2)).*xalfa + a(1); 

% Let's determine q so that it interpolates (x(i),pvals(i)), i=1:4

% q(x(i)) = p(x(i)) is a linear equation in c(1),c(2),c(3), and c(4):

%      c(1) + c(2)(x(i) - beta) + c(3)(x(i)-beta)^2 + c(4)(x(i) - beta)^3 = p(x(i))

% Set up the matrix of coefficients and solve for c(1:4):

xbeta = x-beta;
V = [ones(4,1) xbeta  xbeta.^2  xbeta.^3];
c = V\pvals;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Problem 5

n = 7;
p = 2;
A = triu(tril(randn(n,n),p),-1);
[L,U] = SpecLU(A,p)
Error = A - L*U

function [L,U] = SpecLU(A,p)
% A is an n-by-n matrix with lower bandwidth 1 and upper bandwidth p.
% Assume that A has an LU factorization.
% Computes the factorization A = LU where L is an n-by-n unit lower unit 
% triangular matrix and U is an n-by-n upper triangular. 

   [n,n] = size(A); v = zeros(n,1);
   for k=1:n-1
      v(k+1) = A(k+1,k)/A(k,k);
      r = min(n,k+p);
      A(k+1,k:r) = A(k+1,k:r) - v(k+1)*A(k,k:r);
   end
   L = diag(ones(n,1)) + diag(v(2:n),-1); U = triu(A);


%%%%%%%%%%%%%%%%%%%%%%%%

% Problem 6

a = 10;
b = 30;
T = pi;
tol = .0001;
fname = 'P6F';  % sin(2x)

% Let I(L,R) be the integral from L to R of P6F(x)

k = floor((b-a)/T);         % Number of whole periods encountered from a to b
c = b-k*T;      

% I(a,b) = k*I(a,a+T) + I(a,c)
%        = k*(I(a,c) + I(c,a+T)) + I(a,c)
%        = (k+1)*I(a,c) + k*I(c,a+T)

% If  e1 bound the absolute error in our QUAD estimate of I(a,c)
% and e2 bound the absolute error in our QUAD estimate of I(c,a+T)
% then total error <= (k+1)*e1 + k*e2

% We want the total absolute error to be <= tol. 

% One choice: e1 = e2 = .5*tol/(k+1);

I_clever = (k+1)*QUAD('P6F',a,c,[0 .5*tol/(k+1)]) + k*QUAD('P6F',c,a+T,[0 .5*tol/(k+1)])

% Compare...

I_simple = QUAD('P6F',a,b,[0 tol])
I_exact  = .5*(-cos(2*b) + cos(2*a))