Chapter 3 M-Files

M-File Home

Script File Index

Function File Index

ShowPWL1

% Script File: ShowPWL1
% Convergence of the piecewise linear interpolant to
% humps(x) on [0,3]

close all
z = linspace(0,3,200)';
fvals = humps(z);
for n = [5 10 25 50]
   figure
   x = linspace(0,3,n)';
   y = humps(x);
   [a,b] = pwL(x,y);
   Lvals = pwLEval(a,b,x,z);
   plot(z,Lvals,z,fvals,'--',x,y,'o');
   title(sprintf('Interpolation of humps(x) with pwL, n = %2.0f',n))
end

ShowPWL2

% Script File: ShowPWL2
% Compares pwLStatic and pwLAdapt on [0,3] using the function
%
%   humps(x) = 1/((x-.3)^2 + .01)  +  1/((x-.9)^2+.04)
%

close all

% Second derivative estimate based on divided differences
z = linspace(0,1,101);
humpvals = humps(z);
M2 = max(abs(diff(humpvals,2)/(.01)^2));
for delta = [1 .5 .1 .05 .01]
   figure
   [x,y] = pwLStatic('humps',M2,0,3,delta);
   subplot(1,2,1)
   plot(x,y,'.');
   title(sprintf('delta = %8.4f Static  n= %2.0f',delta,length(x)))
   [x,y] = pwLAdapt('humps',0,humps(0),3,humps(3),delta,.001);
   subplot(1,2,2)
   plot(x,y,'.');
   title(sprintf('delta = %8.4f  Adapt n= %2.0f',delta,length(x))) 
end 

ShowHermite

% Script File: ShowHermite
% Plots a succession of cubic interpolants to cos(x).
% x(2) converges to x(1) = 0 and x(3) converges to x(4) = 3pi/2.

close all
z = linspace(-pi/2,2*pi,100);
CosValues = cos(z);
for d = [1 .5  .3  .1  .05 .001]
   figure
   xvals = [0;d;(3*pi/2)-d;3*pi/2];
   yvals = cos(xvals);
   c = InterpN(xvals,yvals);
   CubicValues = HornerN(c,xvals,z); 
   plot(z,CosValues,z,CubicValues,'--',xvals,yvals,'*')
   axis([-.5 5 -1.5 1.5])
   title(sprintf('Interpolation of cos(x). Separation = %5.3f',d))
end

ShowPWCH

% Script File: ShowpwCH
% Convergence of the piecewise cubic hermite  interpolant to
% exp(-2x)sin(10*pi*x) on [0,1].)

close all
z = linspace(0,1,200)';
fvals = exp(-2*z).*sin(10*pi*z);
for n = [4 8 16 24]
   x = linspace(0,1,n)';
   y = exp(-2*x).*sin(10*pi*x);
   s = 10*pi*exp(-2*x).*cos(10*pi*x)-2*y;
   [a,b,c,d] = pwC(x,y,s);
   Cvals = pwCEval(a,b,c,d,x,z);
   figure
   plot(z,fvals,z,Cvals,'--',x,y,'*');
   title(sprintf('Interpolation of exp(-2x)sin(10pi*x) with pwCH, n = %2.0f',n))
end
legend('e^{-2z}sin(10\pi z)','The pwC interpolant')

ShowSpline

 
% Script File: ShowSpline
% Illustrates CubicSpline with various end conditions, some not so good.

close all
xvals = linspace(0,4*pi,100);
yvals = sin(xvals);
for n = 4:3:16
   x = linspace(0,4*pi,n)';
   y = sin(x);
   [a,b,c,d] = CubicSpline(x,y,1,1,1);
   svals = pwCEval(a,b,c,d,x,xvals);
   figure
   plot(xvals,yvals,xvals,svals,'--')
   title(sprintf('''Good'' Complete spline interpolant of sin(x) with %2.0f subintervals',n-1))
end
for n = 4:3:16
   x = linspace(0,4*pi,n)';
   y = sin(x);
   [a,b,c,d] = CubicSpline(x,y,1,-1,-1);
   svals = pwCEval(a,b,c,d,x,xvals);
   figure
   plot(xvals,yvals,xvals,svals,'--')
   title(sprintf('''Bad'' Complete spline interpolant of sin(x) with %2.0f subintervals',n-1))
end
for n = 4:3:16
   x = linspace(0,4*pi,n)';
   y = sin(x);
   [a,b,c,d] = CubicSpline(x,y,2,0,0);
   svals = pwCEval(a,b,c,d,x,xvals);
   figure
   plot(xvals,yvals,xvals,svals,'--')
   title(sprintf('Natural spline interpolant of sin(x) with %2.0f subintervals',n-1))
end
for n = 4:3:16
   x = linspace(0,4*pi,n)';
   y = sin(x);
   [a,b,c,d] = CubicSpline(x,y);
   svals = pwCEval(a,b,c,d,x,xvals);
   figure
   plot(xvals,yvals,xvals,svals,'--')
   title(sprintf('Not-a-Knot spline interpolant of sin(x) with %2.0f subintervals',n-1))
end

ShowSplineErr

% Script File: ShowSplineErr
% Examines error in the not-a-knot spline interpolant of 
% sin(2pi*x) across the interval [0,1]. Two equally-spaced knot
% cases are plotted: h =.05  and h = .005.  

close all
z = linspace(0,1,500)';
SinVals = sin(2*pi*z);
k=0;
for n=[21 201 ]

   x = linspace(0,1,n)';
   y = sin(2*pi*x);
   [a,b,c,d] = CubicSpline(x,y);
   Svals = pwCEval(a,b,c,d,x,z);
   k=k+1;
   subplot(1,2,k)
   semilogy(z,abs(SinVals-Svals)+eps);
   axis([0 1 10^(-12) 10^(-3)])
   xlabel(sprintf('Knot Spacing = %5.3f',1/(n-1)))
end

ShowSplineTools

% Script File: ShowSplineTools
% Illustrates the Matlab functions spline, ppval, mkpp, unmkpp

close all

% Set Up Data:
n = 9;
x = linspace(-5,5,n);
y = atan(x);

%   Compute the spline interpolant and its derivative:
S  = spline(x,y);	
[x,rho,L,k] = unmkpp(S);
drho = [3*rho(:,1) 2*rho(:,2) rho(:,3)];
dS = mkpp(x,drho);
 
%   Evaluate S and dS:
z = linspace(-5,5);
Svals = ppval(S,z);
dSvals = ppval(dS,z);

%   Plot: 
atanvals = atan(z);	
figure
plot(z,atanvals,z,Svals,x,y,'*');
title(sprintf('n = %2.0f Spline Interpolant of atan(x)',n))
	
datanvals = ones(size(z))./(1 + z.*z);
figure
plot(z,datanvals,z,dSvals)
title(sprintf('Derivative of n = %2.0f Spline Interpolant of atan(x)',n)) 

pwL

function [a,b] = pwL(x,y)
% [a,b] = pwL(x,y)
% 
% Generates the piecewise linear interpolant of the data specified by the
% column n-vectors x and y. It is assumed that x(1) < x(2) < ... < x(n).
%
% a and b are column (n-1)-vectors with the property that for i=1:n-1, the
% line L(z) = a(i) + b(i)z passes though the points (x(i),y(i)) and 
(x(i+1),y(i+1)).

n = length(x); 
a = y(1:n-1);
b = diff(y) ./ diff(x);

Locate

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 

pwLEval

  function LVals = pwLEval(a,b,x,zVals)
% LVals = pwLEval(a,b,x,zvals)
% Evaluates the piecewise linear polynomial defined by the column (n-1)-vectors
% a and b 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)].
%
% LVals is a column m-vector with the property that LVals(j) = L(zVals(j)) 
% for j=1:m where L(z)= a(i) + b(i)(z-x(i)) for x(i)<=z<=x(i+1).

m = length(zVals); 
LVals = zeros(m,1); 
g = 1;
for j=1:m
   i = Locate(x,zVals(j),g);
   LVals(j) = a(i) + b(i)*(zVals(j)-x(i));
   g = i;
end 

pwLStatic

  function [x,y] = pwLStatic(fname,M2,alpha,beta,delta)
% [x,y] = pwLStatic(fname,M2,alpha,beta,delta)
% Generates interpolation points for a piecewise linear approximation of
% prescribed accuracy. 
%
% fname is string that names an available function f(x).
% Assume that f can take vector arguments.
% M2 is an upper bound for|f"(x)| on [alpha,beta].
% alpha and beta are scalars with alpha<beta.
% delta is a positive scalar.
%
% x and y  column n-vectors with the property that y(i) = f(x(i)), i=1:n. 
% The piecewise linear interpolant L(x) of this data satisfies 
% |L(z) - f(z)| <= delta for x(1) <= z <= x(n).

n = max(2,ceil(1+(beta-alpha)*sqrt(M2/(8*delta))));
x = linspace(alpha,beta,n)';
y = feval(fname,x); 

pwLAdapt

   function [x,y] = pwLAdapt(fname,xL,fL,xR,fR,delta,hmin)
% [x,y] = pwLAdapt(fname,xL,fL,xR,fR,delta,hmin)
% Adaptively determines interpolation points for a piecewise linear
% approximation of a specified function.
%   
% fname is a string that specifies a function of the form y = f(u).
% xL and xR are real scalars and fL = f(xL) and fR = f(xR).
% delta and hmin are positive real scalars that determine accuracy.
%
% x and y are column n-vector with the property that
%              xL = x(1) < ... < x(n) = xR
% and y(i) = f(x(i)), i=1:n. Each subinterval [x(i),x(i+1)] is
% either <= hmin in length or has the property that at its midpoint m,
% |f(m) - L(m)| <= delta where L(x) is the line that connects (x(i),y(i)) 
% and (x(i+1),y(i+1)).
 
if (xR-xL) <= hmin
   % Subinterval is acceptable
   x = [xL;xR]; 
   y = [fL;fR];
else
   mid  = (xL+xR)/2; 
   fmid = feval(fname,mid);
   if (abs(((fL+fR)/2) - fmid) <= delta )
      % Subinterval accepted. 
      x = [ xL;xR];
      y = [fL;fR];
   else
      % Produce left and right partitions, then synthesize.
      [xLeft,yLeft]   = pwLAdapt(fname,xL,fL,mid,fmid,delta,hmin);
      [xRight,yRight] = pwLAdapt(fname,mid,fmid,xR,fR,delta,hmin);
      x = [ xLeft;xRight(2:length(xRight))];
      y = [ yLeft;yRight(2:length(yRight))];
   end
end

HCubic

  function [a,b,c,d] = HCubic(xL,yL,sL,xR,yR,sR)
% [a,b,c,d] = HCubic(xL,yL,sL,xR,yR,sR)
% Cubic Hermite interpolation
% 
% (xL,yL,sL) and (xR,yR,sR) are x-y-slope triplets with xL and xR distinct.
% a,b,c,d are real numbers with the property that if
%            p(z) = a + b(z-xL) + c(z-xL)^2 + d(z-xL)^2(z-xR)
% then p(xL)=yL, p'(xL)=sL, p(xR)=yR, and p'(xR)=sR.

a = yL; 
b = sL;
delx = xR - xL;  
yp = (yR - yL)/delx;
c  = (yp - sL)/delx;
d  = (sL - 2*yp + sR)/(delx*delx); 

pwC

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

  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 

CubicSpline

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);