Chapter 3.   Piecewise Polynomial Interpolation

Piecewise polynomial interpolation:

Given a partition on an interval [a, b]:

a=x₁<x₂<…<x_n=b

The piecewise interpolation means interpolating a function on each subinterval [x_i, x_i+1] with a polynomial.

3.1   Piecewise Linear Interpolation

Assume that x(1:n) and y(1:n) are given. If you connect the points (x₁, y₁) … (x_n, y_n)with straight lines, then it is the piecewise linear interpolation.

1. How TO SET UP the piecewise linear polynomials? 

On each subinterval, we use (ai, bi) to represent the linear function

Li(x)=ai+bi(z-xi).

Then the vectors a=[a(1),a(2),…,a(n-1)] and [b(1),b(2),…,b(n-1)] represent the interpolating piecewise polynomial.

ai=y_i ,   bi=(y_i+1-y_i)/(x_i+1-x_i)

The function pwL performs piecewise linear interpolation.

function [a,b] = pwL(x,y)

%

% Pre:

%   x   column n-vector with x(1) < x(2) < ... < x(n)

%   y   column n-vector.

%

% Post:

%  a,b  column (n-1)-vectors with the property that 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);

Example.

Interpolate sin(2*pi*x) on a uniform, nine-point partition of [0,1].

x=linspace(0,1,9);

[a,b]=pwl(x,sin(2*pi*x))

a =

  Columns 1 through 7

         0    0.7071    1.0000    0.7071    0.0000   -0.7071   -1.0000

  Column 8

   -0.7071

b =

  Columns 1 through 7

    5.6569    2.3431   -2.3431   -5.6569   -5.6569   -2.3431    2.3431

  Column 8

    5.6569

2. EVALUATION

Evaluating a piecewise linear polynomial contains two steps.  First, determine the location of x. Second, use the corresponding linear polynomial to obtain the value.

· Find the location of x.

Method 1. Search the subinterval which houses x successively.

if z==x(n)

     i=n-1;

else

i=sum(x<=z);

end

The flops time is ~i.

Method 2. Binary search.

if z==x(n)

i=n-1;

else

Left=1;

Right=n;

while Right>Left+1

mid=floor((Left+Right)/2);

if  z<x(mid)

Right=mid;

else

Left=mid;

end

end

i=left;

end

The flops time ~log₂(n).

The search can be added a guess number.

function i = Locate(x,z,g)

%

% Pre:

%   x   column n-vector with x(1) < x(2) <...<x(n).

%   z   scalar with x(1) <= z <= x(n).

%   g   optional 3rd argument that satisfies 1 <= g <= n-1.

%

% Post:

%   i   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

· Evaluate.

When x is in the i^th subinterval, use a(i)+b(i)*(z-x(i))to evaluate the function value. For vector z, we use vectorization.

function LVals = pwLEval(a,b,x,zvals)

%

% Pre:

%      x   column n-vector with x(1) < ... < x(n).

%  zvals   column m-vector with each component in [x(1),x(n)].

%    a,b   column (n-1)-vectors

%

% Post:

%  LVals   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

Example: Interpolating the function Humps.

The function humps  is defined by the following.

function [out1,out2] = humps(x)

%HUMPS  A function used by QUADDEMO, ZERODEMO and FPLOTDEMO.

%   Y = HUMPS(X) is a function with strong maxima near x = .3

%   and x = .9. 

%

%   [X,Y] = HUMPS(X) also returns X.  With no input arguments,

%   HUMPS uses X = 0:.05:1.

%

%   Example:

%      plot(humps)

%

%   See QUADDEMO, ZERODEMO and FPLOTDEMO.

%   Copyright (c) 1984-97 by The MathWorks, Inc.

%   $Revision: 5.3 $  $Date: 1997/04/08 05:34:37 $

if nargin==0, x = 0:.05:1; end

y = 1 ./ ((x-.3).^2 + .01) + 1 ./ ((x-.9).^2 + .04) - 6;

if nargout==2,

  out1 = x; out2 = y;

else

  out1 = y;

end

x=linspace(0,3,200);

plot(x,humps(x))

Its piecewise linear interpolation is shown as follows.

% Script File: ShowPWL1

%

% Convergence of the piecewise linear interpolant to

% humps(x) on [0,2*pi].

z = linspace(0,3,200)';

  fvals = humps(z);

  for n = [5 10 25 50]

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

     pause

  end

3. A Priori Determination of the Braekpoits.

Given the tolerance d, try to find the piecewise linear interpolation such that the error is no more that d.

Error Estimate.   |f(z)-L(z)|<=M₂h²/8.

When the partition is uniform, h=(b-a)/(n-1).  Hence,

n>=1+(b-a)*sqrt(m2/(8*d)).

The following function finds the breakpoints for a given function with the designed torelance.

function [x,y] = pwLStatic(fname,M2,alpha,beta,delta)

%

% Pre:

%  fname        string that names an available function f(x).

%               Assume that f can take vector arguments.

%  M2           an upper bound for|f"(x)| on [alpha,beta].

%  alpha,beta   scalars with alpha<beta

%  delta        positive scalar

%

% Post:

%  x,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(x) - f(x)| <= delta, where

%               x = linspace(alpha,beta,n).

%

  n = max(2,ceil(1+(beta-alpha)*sqrt(M2/(8*delta))));

  x = linspace(alpha,beta,n)';

  y = feval(fname,x);

4. Adaptive Piecewise Linear Interpolation

Purpose: We try to use less breakpoints while get the good interpolation.

Fact: Breakpoints that need in an even area of the function are less that those in the rapidly changed area.

· How to determine a sub-partition is necessary or not for a subinterval.

(1) Tolerance: If |f((xL+xR)/2-(f(xL)+f(xR))/2|<=d ?

(2) Minimum length of the interval: If xR-xL<=h_min?

function [x,y] = pwLAdapt(fname,xL,fL,xR,fR,delta,hmin)

%

%   Pre:

%  fname       string that names an available function of the form

%              y = f(u).

%     xL       real

%     fL       real, fL = f(xL)

%     xR       real

%     fR       real, fR = f(xR)

%  delta       positive real

%   hmin       positive real

%

%  Post:

%      x       column n-vector with the property that

%              xL = x(1) < ... < x(n) = xR. Each subinterval

%              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 (xR,fR).

%                

%      y       column n-vector with the property  that

%                       y(i) = f(x(i)).

%          

%

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

Example.

% Script File: ShowPWL2

%

% Compares pwLStatic and pwLAdapt on [0,1] 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,1,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),1,humps(1),delta,.001);

     subplot(1,2,2)

   plot(x,y,'.');

     title(sprintf('delta = %8.4f  Adapt n= %2.0f',delta,length(x)))

  end

3.2 Piecewise Cubic Hermite Interpolation

Assume that the function y=f(x) is given on an interval [x₁,x_n] and f(x) is differentialable. Let

   x₁<x₂ <^…<x_n

is a partition of [x₁, x_n]. The piecewise cubic Hermite interpolation for the function f(x) on the partition is the piecewise cubic polynomial H(x), which has the following properties.

1. On each subinterval [x_i,x_i+1], H(x)  is a cubic polynomial.

2. Over the interval [x₁,x_n], H'(x) exists and is continuous.

3. H(x_i )=f(x_i) and H'(x_i)=f'(x_i), for i=1,2,^…,n.

3.2.1   Construct Cubic Hermite Polynomial Interpolation On A Subinterval

1. Represent a cubic polynomial on [x_i, x_i+1] in the following form

q(z)=a+b(x-x_i)+c(x-x_i)²+d(x-x_i)²(x-x_i+1) (=H(z) on [x_i, x_i+1] )

2. Use the conditions q(x_i)=f(x_i), q(x_i+1)=f(x_i+1), q'(x_i)=f'(x_i), and q'(x_i+1)=f(x_i+1) to set a linear system for [a, b, c, d].

a=f(x_i)     a+b(x_i+1-x_i)+c(x_i+1-x_i)²=f(x_i+1)

b=f'(x_i)      b+2c(x-x_i)+d(x-x_i)²=f'(x_i+1)

3. Solve the linear system.

a=f(x_i)

b=f'(x_i)

y'_i=(f(x_i+1)-f(x_i))/(x_i+1-x_i)

c=(y'_i-b)/(x_i+1-x_i) 

d=(f'(x_i+1)+f'(x_i)-y'_i)/(x_i+1-x_i)

The following is the corresponding program.

function [a,b,c,d] = HCubic(xL,yL,sL,xR,yR,sR)

%

% Pre:

%    xL,xR   distinct reals

%

% Post:

%  a,b,c,d   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, 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);

· The interpolation error:

Theorem: Suppose f(z) and its first four derivatives are continuous on [x_L, x_R] and  |f⁽⁴⁾(z)|<= M₄ for all z in [x_L, x_R]. Then

|f(z)-q(z)|<=M₄(x_R-x_L)⁴/384.

3.3.2 Put All Pieces Together 

Using the vectorization, we have the following.

  function [a,b,c,d] = pwC(x,y,s)

%

% Pre:

%      x,y,s  column n-vectors with x(1) < ... < x(n)

%

% Post:

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

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

3.2.3 Evaluation

The evaluation contains two steps, as we illustrated in the last section. First, locate the positions of x-coordinates of the input. Second, using the right piece of the polynomial to obtain the value.

function Cvals = pwCEval(a,b,c,d,x,z)

%

% Pre:

%    a,b,c,d   column (n-1)-vectors that represent a piecewise

%              cubic polynomial.

%    x         column n-vector with x(1) < ... < x(n).

%    z         column m-vector with each z(i) in [x(1),x(n)].

%

% Post:

%    Cvals     column m-vector with the property that

%              Cvals(j) = C(z(j)), j=1:m. If z(j) is in

%              the interval [x(i),x(i+1)], then

%              C(z(j)) = a(i) +

%                        b(i)(z(j)-x(i)) +

%                        c(i)(z(j)-x(i))^2 +

%                        d(i)(z(j)-x(i))^2(z(j)-x(i+1))

%

  m = length(z);

  Cvals = zeros(m,1);

  g=1;

  for j=1:m

     i = Locate(x,z(j),g);

     Cvals(j) = d(i)*(z(j)-x(i+1)) + c(i);

     Cvals(j) = Cvals(j)*(z(j)-x(i)) + b(i);

     Cvals(j) = Cvals(j)*(z(j)-x(i)) + a(i);

     g = i;

end

3.3 Cubic Splines

3.3.1      Cubic Spline Interpolation

Assume that the function y=f(x) is given on an interval [x₁,x_n] and f(x) is differentiable. Let

   x₁<x₂ <^…<x_n

is a partition of [x₁,x_n]. The cubic spline interpolation for the function f(x) on the partition is the piecewise cubic polynomial S(x) which has the following properties.

1. On each subinterval [x_i, x_i+1], S(x)  is a cubic polynomial.

2. Over the interval [x₁, x_n], S''(x) exists and is continuous.

3. S(x_i )=f(x_i) for i=1,2,^…,n.

4. At the end points x₁ and x_n, S(x) satisfies the end conditions which are explained in the next subsections.

3.3.2 The End Conditions

1. Complete Spline: S'(x₁)=muL, S'(x_n)=muR.

2. Second –type End Conditions and The Natural Splines: S''(x₁)=muL, S''(x_n)=muR.

When we set S''(x₁)=0, S''(x_n)=0, we get the natural spline.

3. The Not-a-Knot Spline: Let S'''(x) be continuous at x₂ and x_n-1. ( This   removes the breakpoints at x₂ and x_n-1).

3.3.3 Usage Of The Function Cubicspline (built-in)

Function cubicspline(x, y, derivative, muL, muR) is used to construct the cubic spline interpolant with any of the three types of end conditions. The explanation of the inputs can refer the following.

[a, b, c, d]=cubicspline(x, y, derivative, muL, muR)

% x,y  column n-vectors. n >= 4 and x(1) < ... x(n)

% derivative is the order of the spline's derivative that      % are used in the end conditions (= 1 or 2).

% muL,muR are the value of the specified derivative at the

% left and right endpoint.

%  a,b,c,d all are column (n-1)-vectors that define the % % % spline. On [x(i),x(i+1)], the spline S(z) is specified by % the cubic

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

3.3.4 MATLAB Spline Tools

Function spline(x,y,z) can be used to create not-a-knot spline interpolants, where x and y are the data to be inperplated, and z be the vector on which the spline takes the values.

Example.

z=linspace(-5,5);

x=linspace(-5,5,9);

y=atan(x);

Svals=spline(x,y,z);

plot(z,Svals)

Function ppval(pp,x)evaluates piecewise polynomial. It returns the value of the piecewise polynomial pp at the points x.  The piecewise polynomial form (pp-form) is

    returned by function spline.

Example.

The following will perform the same as in the last example.

x=linspace(-5,5,9);

y=atan(x);

S=spline(x,y);

z=linspace(-5,5);

Svals=ppval(S,z);

plot(z,Svals)

Please use help to find the usage of other two useful functions mkpp and unmkpp.