P4: Given Scripts/Functions

  function G = CholSax(A)
% G = CholSax(A)
% Cholesky factorization of a symmetric and positive definite matrix A.
% G is lower triangular so A = G*G'.

[n,n] = size(A);
G = zeros(n,n);
s = zeros(n,1);
for j=1:n
   s(j:n) = A(j:n,j);
   for k=1:j-1
      s(j:n) = s(j:n) - G(j:n,k)*G(j,k);
   end
   G(j:n,j) = s(j:n)/sqrt(s(j));
end

  function x = LTriSolBAnd(L,b,p)
% x = LTriSolBand(L,b,p)
%
% Solves the nonsingular lower triangular system  Lx = b 
% where L is n-by-n with lower bandwidth p, b is n-by-1, 
% and x is n-by-1.

n = length(b);
x = zeros(n,1);
for j=1:n-1
   x(j) = b(j)/L(j,j);
   r = min(j+p,n);
   b(j+1:r) = b(j+1:r) - L(j+1:r,j)*x(j);
end
x(n) = b(n)/L(n,n);

  function x = UTriSolBand(U,b,p)
% x = UTriSol(U,b)
%
% Solves the nonsingular upper triangular system  Ux = b.
% where U is n-by-n with upper bandwidth p, b is n-by-1, 
% and x is n-by-1.

n = length(b);
x = zeros(n,1);
for j=n:-1:2
   x(j) = b(j)/U(j,j);
   r = max(1,j-p);
   b(r:j-1) = b(r:j-1) - x(j)*U(r:j-1,j);
end
x(1) = b(1)/U(1,1);

% P4A

% Check CholSaxBand with an n = 7, p = 3 example...

clc
n = 7; p = 3;
randn('state',0); B = randn(n,n);
A = B'*B + n*eye(n,n);
A = tril(triu(A,-p),p);
G1 = CholSax(A)
G2 = CholSaxBand(A,p)

% Script P4B

close all

% Check that BEval is working.

a = -3; b = 3;
n = 7;
plot([-3 3],[0 0])
axis([-3 3 -.5 1.25])
hold on
M = 100;
for k=-3:2
   z = linspace(k,k+1,M+1)';
   F = zeros(M,4);
   for i=1:M
      [F(i,:),j] = BEval(z(i),a,b,n);
   end
   plot(z(1:M),F,'k')
end
hold off
title('Basis Splines B_{-3},...,B_{3} for the case x(1) = 0, h = 1')

% Now work with the Sunspot data

y = GetSunSpotData;
m = length(y);
z = linspace(1874,2002,m);

for n=[20 22 50 ]
   % Fit the data with a linear combination of n+2 basis splines
   alfa = LSSpline(z,y,n);
   for i=1:m
      % Evaluate those basis splines that are nonzero at z(i).
      [f,j] = BEval(z(i),z(1),z(m),n);
      g(i) = f*alfa(j:j+3);
   end
   % Display the LS spline on top of the sunspot data...
   figure
   plot(z,y,'y',z,g)
   axis([1875 2002 0 4000])
   title(sprintf('n = %2d',n))
end

function z = GetSunSpotData
% for y = 1875:2001 and j=1:12
%   z(12*(y-1875)+j) = amount of sunspot activity during month j
%   of year y.
fid = fopen('SunSpotArea.txt');
n = 12*(2001-1875+1);
z = zeros(n,1);
for k=1:n
   x = fscanf(fid,'%g',1);
   y = fscanf(fid,'%g',1);
   z(k) = fscanf(fid,'%g',1);
end

    function [f,j] = BEval(zVal,a,b,n)
  % zVal is a scalar that satisfies a <= zVal <= b.
  % n is an integer with 2 <= n.
  % Let h = (b-a)/(n-1) and define x(k) = a + (k-1)h for any integer k.
  % j is an integer with 1 <= j <= n-1 such that
  %              x(j) <= zVal <= x(j+1).
  % f is a row 4-vector that houses the values of the four nonzero 
  % basis splines at zVal, i.e.,
  %      f = [B_{j-1} (zVal) B_{j}(zVal)  B_{j+1}(zVal)  B_{j+2}(zVal) ].

    function alfa = LSSpline(z,y,n)
  % z and y are column m-vectors and 2 <= n <= m-2.
  % Assume z(1) < z(2) < ... < z(m).
  % alfa is a column (n+2)-vector with the property that if
  %
  %                 alfa' = [alpha(0),...,alpha(n+1)]
  %
  % then 
  %
  %           alpha(0)B_{0}(x) + alpha_{1}B_{1}(x) + ... + alpha_{n+1}B_{n+1}(x)
  %
  % is the least squares fit of (z(1),y(1)),...,(z(m),y(m)). Here, B_{0},...,B_{n+1}
  % are the basis splines associated with x(j) = z(1) + (j-1)h, j=0:n+1, 
  % h = (z(m)-z(1))/(n-1).