Hermite Cubic and
B-spline Interpolations (Lec 11)

yukita@k.hosei.ac.jp

•Hermitian Cubic Interpolation

We wish to find a piecewise polynomial of degree 3 that passes through the data points and is also continuously differentialble at these points.  We also prescribe  the values of the derivatives at the given data points.

Definition:

Her[t_]=Which[ 0<=t<1, 2t^3-3t^2+1,
               -1<t<0, -2t^3-3t^2+1,
               True, 0 ];

DHer[t_]=Which[ 0<=t<1, t^3-2t^2+t,
                -1<t<0, t^3+2t^2+t,
                True, 0 ];

Plot[ {Her[t], DHer[t]}, {t,-1,1} ];

[Graphics:HTMLFiles/index_1.gif]

valuesAt[t_]={ {Her[t], Derivative[1][Her][t] },
              {DHer[t], Derivative[1][DHer][t]} };

valuesAt[-1]

{{0, 0}, {0, 0}}

valuesAt[0]

{{1, 0}, {0, 1}}

valuesAt[1]

{{0, 0}, {0, 0}}

With the above observations in mind, we introduce a set of functions
{ Her[[t-j], DHer[t-j] ;   j=0,  1, 2, 3, ....  }

The functions {Her[t-j]} give a partition of unitiy on a finite interval.

Plot[ Sum[ Her[t-j], {j,0,10}], {t,-2, 12} ];

[Graphics:HTMLFiles/index_2.gif]

Let's plot and see them with a concrete example.

•Example

We will interpolate  the following data.

points = {{-1,-1},{0, -1},{2,1},{4,4},{5,2},{7,0}};

velocities = {{1,1},{1,0},{2,1},{2,0},{1,-1},{1,0}};

P[t_]=Sum[ Her[t-j] points[[j+1]] +
           DHer[t-j] velocities[[j+1]], {j,0,5}];

g1=ParametricPlot[ Evaluate[P[t]], {t,0,5},
                   DisplayFunction->Identity];

g00=Graphics[ Table[ Circle[points[[j+1]], 0.1 ],{j,0,5}] ];

Needs["Graphics`Arrow`"];
g01=Graphics[ Table[ Arrow[ points[[j+1]],
                            points[[j+1]]+velocities[[j+1]] ],
                     {j,0,5} ] ];

Show[{g1,g00,g01}, DisplayFunction->$DisplayFunction];

[Graphics:HTMLFiles/index_3.gif]

The following properties give us a nice localized control over the interplotation.

Her[t-j]        /. t->j                   ==> 1
D[Her[t-j],t]   /. t->j                   ==> 0
Her[t-j]        /. t->i (!=j)             ==> 0
D[Her[t-j],t]   /. t->i (!=j)             ==> 0

DHer[t-j]       /. t->j                   ==> 0
D[DHer[t-j],t]  /. t->j                   ==> 1
DHer[t-j]       /. t->i (!=j)             ==> 0
D[DHer[t-j],t]  /. t->i (!=j)             ==> 0
(* psuedo codes, don't evaluate them *)

We consider the two linear mappings.

interp : R^(2n+2) ->  PCF
      {y0,y1, ..., yn, y0',y1',...,yn'}
      |->  P[t]=Sum[ yj Her[t-j] , {j,0,n}]
               +Sum[ yj' DHer[t-j],{j,0,n}]
(* pseudo code, don't evaluate it *)

where PCF  is a real vector space of piecewise cubic functions.

proj : PCF -> R^(2n+2)
         P[t] |-> { P[0], P[1], ..... , P[n],
                    P'[0],P'[1],..... , P'[n] }
(* pseudo code, don't evaluate it *)

It is easy to see that

proj@@interp : R^(2n+2) -> R^(2n+2)
(* pseudo code, don't evaluate it *)

is the identity map.

•B-spline Interpolation

Recall the definition of basis functions for B-splines.

B[1][t_]=Which[ -1<t<0, t+1, 0 <= t <1, -t+1, True, 0];

B[q_][t_]:=1/q ( (t+(q+1)/2) B[q-1][t+1/2] +
                ((q+1)/2-t) B[q-1][t-1/2] );

We will interpolate the data below.

data={{0,0},{1,3},{3,3},{5,-2},{6,-2},{7,+1},{9,4},{10,4}};

cp[j_]={cpx[j],cpy[j]};

P[t_]=Sum[ B[2][t-j] cp[j], {j,0,7}]+
           B[2][t+1] cp[0] +
           B[2][t-8] cp[7];

eq = Table[ P[j]==data[[j+1]], {j,0,7} ]

  7 cpx[0]   cpx[1]  7 cpy[0]   cpy[1]
{{-------- + ------, -------- + ------} == {0, 0}, 
     8         8        8         8
 
   cpx[0]   3 cpx[1]   cpx[2]
  {------ + -------- + ------, 
     8         4         8
 
    cpy[0]   3 cpy[1]   cpy[2]
    ------ + -------- + ------} == {1, 3}, 
      8         4         8
 
   cpx[1]   3 cpx[2]   cpx[3]
  {------ + -------- + ------, 
     8         4         8
 
    cpy[1]   3 cpy[2]   cpy[3]
    ------ + -------- + ------} == {3, 3}, 
      8         4         8
 
   cpx[2]   3 cpx[3]   cpx[4]
  {------ + -------- + ------, 
     8         4         8
 
    cpy[2]   3 cpy[3]   cpy[4]
    ------ + -------- + ------} == {5, -2}, 
      8         4         8
 
   cpx[3]   3 cpx[4]   cpx[5]
  {------ + -------- + ------, 
     8         4         8
 
    cpy[3]   3 cpy[4]   cpy[5]
    ------ + -------- + ------} == {6, -2}, 
      8         4         8
 
   cpx[4]   3 cpx[5]   cpx[6]
  {------ + -------- + ------, 
     8         4         8
 
    cpy[4]   3 cpy[5]   cpy[6]
    ------ + -------- + ------} == {7, 1}, 
      8         4         8
 
   cpx[5]   3 cpx[6]   cpx[7]
  {------ + -------- + ------, 
     8         4         8
 
    cpy[5]   3 cpy[6]   cpy[7]
    ------ + -------- + ------} == {9, 4}, 
      8         4         8
 
   cpx[6]   7 cpx[7]  cpy[6]   7 cpy[7]
  {------ + --------, ------ + --------} == {10, 4}}
     8         8        8         8

sol=Solve[ eq, Join[Table[cpx[j],{j,0,7}],
                    Table[cpy[j],{j,0,7}] ] ]

              28729              2381735
{{cpx[0] -> -(------), cpx[7] -> -------, 
              235416             235416
 
               115315             924269            201103
   cpy[0] -> -(------), cpy[7] -> ------, cpx[1] -> ------, 
               235416             235416            235416
 
             705439            1216247            1413719
   cpx[2] -> ------, cpx[3] -> -------, cpx[4] -> -------, 
             235416            235416             235416
 
             1601407            2161135            807205
   cpx[5] -> -------, cpx[6] -> -------, cpy[1] -> ------, 
             235416             235416             235416
 
             922069              689635
   cpy[2] -> ------, cpy[3] -> -(------), 
             235416              235416
 
               550915             228469            1063429
   cpy[4] -> -(------), cpy[5] -> ------, cpy[6] -> -------}
               235416             235416            235416
 
   }

g1=ParametricPlot[ Evaluate[ P[t]/. sol[[1]] ],
                   {t,0,7},
                   DisplayFunction->Identity,
                   PlotLabel->"B-spline" ];

g0=Graphics[ Table[ Circle[data[[j+1]], 0.1 ],{j,0,7}] ];

gB=Show[{g1,g0} ];

Show[gB,DisplayFunction->$DisplayFunction];

[Graphics:HTMLFiles/index_4.gif]

It looks nicer than the Lagrange approximation.

Compare the B-spline and Lagrange interpolations.

Show[ GraphicsArray[{gB, gLag}],
      DisplayFunction->$DisplayFunction ];

[Graphics:HTMLFiles/index_5.gif]

•Exercises

•Problem 1

Prove that the base function used in the Hermite interpolation are uniquely determined under a reasonable set of conditions; find two cubic functions f[t] and g[t] that satisfy the following conditions.

f[0]==1, f'[0]=0, f[1]== f'[1]==0,
g[0]==0, g'[0]==1, and g[1]==g'[1]==0.

•Solution

•Problem 2

Using the result of Problem 1, that is, f[t] and g[t], give the definitions of Hermite interpolation basis.

•Solution

Converted by Mathematica  (July 8, 2003)