Problem Set 2 - Given Files

CSInterp.m:

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.m:

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

CubicSpline.m:

function [a,b,c,d] = CubicSpline(x,y,derivative,muL,muR)
% [a,b,c,d] = CubicSpline(x,y,derivative,muL,muR)
% Cubic spline interpolation with prescribed end conditions.
%
% x,y are column n-vectors. It is assumed that n >= 4 and x(1) < ... x(n).
% derivative is an integer (1 or 2) that specifies the order of the
% endpoint derivatives.
% muL and muR are the endpoint values of this derivative.
%
% a,b,c, and d are column (n-1)-vectors that define the spline S(z). On
% [x(i),x(i+1)],
%
%          S(z) =  a(i) + b(i)(z-x(i)) + c(i)(z-x(i))^2 +
% d(i)(z-x(i))^2(z-x(i+1).
%
% Usage:
%   [a,b,c,d] = CubicSpline(x,y,1,muL,muR)   S'(x(1))  = muL, S'(x(n))  = muR
%   [a,b,c,d] = CubicSpline(x,y,2,muL,muR)   S''(x(1)) = muL, S''(x(n)) = muR
%   [a,b,c,d] = CubicSpline(x,y)             S'''(z) continuous at x(2) and x(n-1)

% First, set up all but the first and last equations that
% define the vector of interior knot slopes s(2:n-1).

n = length(x);
Dx = diff(x);
yp = diff(y) ./ Dx;
T = zeros(n-2,n-2);
r = zeros(n-2,1);
for i=2:n-3
   T(i,i)   = 2*(Dx(i) + Dx(i+1));
   T(i,i-1) = Dx(i+1);
   T(i,i+1) = Dx(i);
   r(i)     = 3*(Dx(i+1)*yp(i) + Dx(i)*yp(i+1));
end

% For each of the 3 cases, finish setting up the linear system,
% solve the system, and set s(1:n) to be the vector of slopes.

if nargin==5
   %Derivative information available.
   if derivative==1
      % End values for S'(z) specified.
      T(1,1) = 2*(Dx(1) + Dx(2));
      T(1,2) = Dx(1);
      r(1) = 3*(Dx(2)*yp(1)+Dx(1)*yp(2)) - Dx(2)*muL;
      T(n-2,n-2) = 2*(Dx(n-2)+Dx(n-1));
      T(n-2,n-3) = Dx(n-1);
      r(n-2) = 3*(Dx(n-1)*yp(n-2) + Dx(n-2)*yp(n-1)) -Dx(n-2)*muR;
      s = [muL; T\r; muR];
   end
  
   if derivative==2
      % End values for S''(z) specified.
      T(1,1) = 2*Dx(1) + 1.5*Dx(2);
      T(1,2) = Dx(1);
      r(1) = 1.5*Dx(2)*yp(1) + 3*Dx(1)*yp(2) + Dx(1)*Dx(2)*muL/4;
      T(n-2,n-2) = 1.5*Dx(n-2)+2*Dx(n-1);
      T(n-2,n-3) = Dx(n-1);
      r(n-2) = 3*Dx(n-1)*yp(n-2) + 1.5*Dx(n-2)*yp(n-1)-Dx(n-1)*Dx(n-2)*muR/4;
      stilde = T\r;
      s1 = (3*yp(1) - stilde(1) - muL*Dx(1)/2)/2;
      sn = (3*yp(n-1) - stilde(n-2) + muR*Dx(n-1)/2)/2;
      s = [s1;stilde;sn];
   end;
else
   % No derivative information. Compute the not-a-knot spline.
   q = Dx(1)*Dx(1)/Dx(2);
   T(1,1) = 2*Dx(1) +Dx(2) + q;
   T(1,2) = Dx(1) + q;
   r(1) = Dx(2)*yp(1) + Dx(1)*yp(2)+2*yp(2)*(q+Dx(1));
   q = Dx(n-1)*Dx(n-1)/Dx(n-2);
   T(n-2,n-2) = 2*Dx(n-1) + Dx(n-2)+q;
   T(n-2,n-3) = Dx(n-1)+q;
   r(n-2) = Dx(n-1)*yp(n-2) + Dx(n-2)*yp(n-1) +2*yp(n-2)*(Dx(n-1)+q);
   stilde = T\r;
   s1 = -stilde(1)+2*yp(1);
   s1 = s1 + ((Dx(1)/Dx(2))^2)*(stilde(1)+stilde(2)-2*yp(2));
   sn = -stilde(n-2) +2*yp(n-1);
   sn = sn+((Dx(n-1)/Dx(n-2))^2)*(stilde(n-3)+stilde(n-2)-2*yp(n-2));
   s = [s1;stilde;sn];
end   

% Compute the a,b,c,d vectors.
      
a = y(1:n-1);
b = s(1:n-1);
c = (yp - s(1:n-1)) ./ Dx; 
d = (s(2:n) + s(1:n-1) - 2*yp) ./ (Dx.* Dx);

HornerN.m:

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

HornerV.m:

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

InterpN.m:

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;

InterpV.m:

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;

Locate.m:

function i = Locate(x,z,g)
% i = Locate(x,z,g)
% Locates z in a partition x.
% 
% x is column n-vector with x(1) < x(2) <...<x(n) and
% z is a scalar with x(1) <= z <= x(n).
% g (1<=g<=n-1) is an optional input parameter
% 
% i is an integer such that x(i) <= z <= x(i+1). Before the general
% search for i begins, the value i=g is tried.

if nargin==3
   % Try the initial guess.
   if (x(g)<=z) & (z<=x(g+1))
      i = g;
      return
   end
end   
n = length(x);
if z==x(n)
   i = n-1;
else
   % Binary Search
   Left = 1;
   Right = n;
   while Right > Left+1
      % x(Left) <= z <= x(Right)
      mid = floor((Left+Right)/2);
      if z < x(mid)
         Right = mid; 
      else
         Left = mid;
      end
   end   
   i = Left;
end

pwC.m:

function [a,b,c,d] = pwC(x,y,s)
% [a,b,c,d] = pwC(x,y,s)
% Piecewise cubic Hermite interpolation.
%
% x,y,s  column n-vectors with x(1) < ... < x(n)
%
% a,b,c,d  column (n-1)-vectors that define a continuous, piecewise
% cubic polynomial q(z) with the property that for i = 1:n,
%
%             q(x(i)) = y(i) and q'(x(i)) = s(i).
%
% On the interval [x(i),x(i+1)],
%
%      q(z) = a(i) + b(i)(z-x(i)) + c(i)(z-x(i))^2  + d(i)(z-x(i))^2(z-x(i+1)).

n  = length(x);
a  = y(1:n-1);
b  = s(1:n-1);
Dx = diff(x);
Dy = diff(y);
yp = Dy ./ Dx;
c  = (yp - s(1:n-1)) ./ Dx;
d  = (s(2:n) + s(1:n-1) - 2*yp) ./ (Dx.* Dx);

pwCEval.m:

function Cvals = pwCEval(a,b,c,d,x,zVals)
% Cvals = pwCEval(a,b,c,d,x,zVals)
%
% Evaluates the piecewise cubic polynomial defined by the column
% (n-1)-vectors a,b,c, and d and the column n-vector x. It is assumed that
% x(1) < ... < x(n).
% zVals  is a column m-vector with each component in [x(1),x(n)].
%
% CVals is a column m-vector with the property that CVals(j) = C(zVals(j)) 
% for j=1:m where on the interval [x(i),x(i+1)]
%
%   C(z)= a(i) + b(i)(z-x(i)) + c(i)(z-x(i))^2 + d(i)(z-x(i))^2(z-x(i+1))
 
m = length(zVals); 
Cvals = zeros(m,1); 
g=1;
for j=1:m
   i = Locate(x,zVals(j),g);
   Cvals(j) = d(i)*(zVals(j)-x(i+1)) + c(i);
   Cvals(j) = Cvals(j)*(zVals(j)-x(i)) + b(i);
   Cvals(j) = Cvals(j)*(zVals(j)-x(i)) + a(i);
   g = i;
end