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) 
  

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

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