Chapter 3 M-Files

M-File Home

Script File Index

 ShowPWL1 Illustrates pwL and pwLEval. ShowPWL2 Compares pwLStatic and pwLAdapt. ShowHermite Illustrates the Hermite interpolation idea. ShowpwCH Illustrates pwC and pwCEval. ShowSpline Illustrates CubicSpline. ShowSplineErr Explores the not-a-knot spline interpolant error. ShowSplineTools Illustrates \Matlab spline tools.

Function File Index

 pwL Sets up a piecewise linear interpolant. Locate Determines the subinterval in a mesh that houses a given x-value. pwLEval Evaluates a piecewise linear function. pwLStatic A priori determination of a mesh for a pwL approximation. pwLAdapt Dynamic determination of a mesh for a pwL approximation. HCubic Constructs the cubic Hermite interpolant. pwC Sets up a piecewise cubic Hermite interpolant. pwCEval Evaluates a piecewise cubic function. CubicSpline Constructs complete, natural, or not-a-knot spline.

% 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)))
subplot(1,2,2)
plot(x,y,'.');
title(sprintf('delta = %8.4f  Adapt n= %2.0f',delta,length(x)))
end


% 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


% 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


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


   function [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.
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);


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