Back to November Calendar            ←Previous Entry                             Next Entry→

November 15, 2005

Spline Representations VII

 

The advantage of natural splines is that the curve meeting these conditions minimizes the total curvature as measured by

                            

 The text, Applied Numerical Methods for Engineers is accompanied by a disc containing numerical libraries NLIB.   These work in tandem with MatLab.  After installing the files as a subdirectory of MatLab, you launch MatLab and enter the command, browse, to implement the NLIB/MatLab interface.  This runs the file browse.m, which presents a user with the menu options

 

Run : examples, selected problems, exit

View : examples, selected problems

User : create script, run, view, add, remove

Help : NLIB toolbox, software updates

 

The run option allows the user to execute the computational examples of the text.  The view option allows the user to also see source code and data files. The user option provides a simple way to include user scripts. You can create, run, and view user scripts as well as add them to or remove them from the browser menus. That makes browse an “integrated development environment” (IDE) for high-level program development.  The help option has documentation on all of the functions in the NLIB toolbox. 

 

As such, NLIB is, at first blush, a promising competitor for the more established  Numerical Recipes texts published by Cambridge University Press.   Certainly the link with MatLab sets it apart.  The menu window looks like this:

 

If you choose run/example/example 4.4.2: Cubic Spline Fit, the window disappears and the following code appears in the MatLab command window:

 

Example 4.4.2: Cubic Spline Fit

 

     -------------

    |          -4 |

    |          -3 |

    |          -2 |

    |          -1 |

x = |           0 |

    |           1 |

    |           2 |

    |           3 |

    |           4 |

     -------------

1 => Save, 2 => Append, Enter => continue:

 

Note the include file in the c++ program below, util.h, contains the links to accommodate the user types “vector” and “matrix.”  I’ll take a look.  Here is util.h:

 

/*------------------------------------------------------------------*/

/* File:        util.h */

 

/* Description: Include this file to insert the NLIB header file and

                source files. By using this include file, there is

                no need to use a compiled library.  This is more

                portable than nlib.h, but it generates larger .exe

                files (about 50 %) and it takes somewhat longer to

                compile. */

/*-------------------------------------------------------------------*/

 

/* NLIB header file */

 

#include "c:\nlib\nlib.h"

 

/* Core module, main-program support */

 

#include "c:\nlib\nlib.c"       /* NLIB core functions */

#include "c:\nlib\show.c"       /* tabular display */

#include "c:\nlib\draw.c"       /* graphical display */                   

 

/* Application modules (Chapters 2-10) */

 

#include "c:\nlib\linear.c"     /* linear algebraic systems, 2 */         

#include "c:\nlib\eigen.c"      /* eigenvalues and eigenvectors */        

#include "c:\nlib\curves.c"     /* curve fitting */                       

#include "c:\nlib\roots.c"      /* root finding */                        

#include "c:\nlib\optim.c"      /* optimizaition */                       

#include "c:\nlib\integ.c"      /* differentiation and integrition */     

#include "c:\nlib\ode.c"        /* ordinary differential equations */     

#include "c:\nlib\pde.c"        /* partial differential equations */      

#include "c:\nlib\dsp.c"        /* digital signal processing */           

 

I find the fact that my files are not located at "c:\nlib\nlib.h" just a little disturbing…

 

I can confirm that nlib.c contains all the low level core functions, including functions for creating and initializing vectors and matrices, while show.c, to quote the program itself,

 

provides console support functions to interactively display scalars, vectors, and matrices on the console screen.  It uses platform-specific functions (not ANSI C) that are contained in the console configuration file "console.c".  An appropriate version of console.c must be copied from the console subdirectory to the NLIB directory.

 

In fact, there is no console.c in my NLIB directory.  Maybe the install was not complete?

 

The draw.c function, describes itself this way:

 

This file contains NLIB functions that write data into Matlab executable M-files.  These M-files can be run from Matlab to produce plots and perform additional Matlab processing.

 

In the main program, shown immediately below, the call to spline (whose code follows) comes after the two calls to showvec:

 

/* -------------------------------------------------------------- */   
/* Example 4.4.2 Cubic Spline Fit                                 */

/* -------------------------------------------------------------- */

#include "c:\nlib\util.h"

int  main (void)

{

   int n = 9, m = 10*n, i,j;

   float d = 4, a;                                                        

   vector       x = vec (n,""),

                y = vec (n,""),

                c = vec (n,"");

   matrix       X = mat (2,m,""), Y = mat (2,m,"");

   printf ("\nExample 4.4.2: Cubic Spline Fit\n");

   /* Generate the control points coordinates */

   for (i = 1; i <= n; i++)    {

      a = -d + 2*d*(i-1)/(n-1);

      x[i] = a;

      y[i] = 1/(1 + a*a);

   }

   showvec ("x",x,n,0);

   showvec ("y",y,n,0);

   /* ic = 1 means compute new values for spline coefficients c */

   /* The reason you only need to do this once is that there are */
   /* only n-1 = 8 of these: one for each cubic. */

   spline (0,x,y,c,n,1);

   for (i = 1; i <= m; i++)

   {

      a = x[1] + (i-1)*(x[n] - x[1])/(m-1);

      X[1][i] = a;

      Y[1][i] = spline (a,x,y,c,n,0);

      j = MIN (i,n);

      X[2][i] = x[j];

      Y[2][i] = y[j];

   }

   /* This is the call to plot the graph  subst. OpenGL? */

   graphmm (X,Y,2,m,"Cubic Spline Fit","x","y","fig442",1);

   wait();

   return 0;

}

/* -------------------------------------------------------------- */

 

The function spline appears in the curves.c  file.  Here’s the code for spline:

 

/* -------------------------------------------------------------------- */

float spline (float a, vector x, vector y, vector c, int n, int ic) {

/* -------------------------------------------------------------------- */

 

/* Description: Evaluate a cubic spline function at x.

 

   On entry:    a  = evaluation point

                x  = n by 1 vector of strictly increasing independent

                     variables

                y  = n by 1 vector of dependent variables

 

                c  = n by 1 vector containing the spline coefficients,

                n  = number of sample points (n >= 3)

                ic = initial condition code defined as follows:

 

                    0  =>  use previously computed c

                    1  =>  compute vector c

 

   On exit:    spline = natural cubic spline function evaluated at a

 

   Notes:       1. spline is initialized using ic <> 0 and then

                   subsequent calls use ic = 0.

                2. spline uses the function tridec in linear.c to

                   compute the cubic spline coefficients. */   

           

   int  k = 1, i;

   float  s, b0, b1, p;

   vector  h, w, v;

   matrix  T;

 

   if (n < 3)

      nerror (100,"spline",5,(float) n,"");

  

   if (ic) {   /* Compute spline coefficients */

      h = vec (n-1,"");                                                                                                           

      w = vec (n-1,"");

      v = vec (n,"");                                                                                                             

      T = mat (3,n,"");                                                                                                           

 

      /* Initialize right-hand side v */

      /* First compute the delta x’s*/

      for (i = 1; i < n; i++)

      {

         h[i] = x[i+1] - x[i];

         if (h[i] > 0)  /* then the slopes, w */

            w[i] = (y[i+1] - y[i])/h[i];

         else

         {

            printf ("\nThe independent variables in spline must be");

            printf ("\nstrictly increasing.\n");

            return 0;

         }

      }

      v[1] = 0;

      v[n] = 0;

      for (i = 2; i < n; i++) /* 6 times the change in the slopes */

         v[i] = 6*(w[i] - w[i-1]);

 

      /* Solve equations */

      for (i = 1; i < n; i++)

      {

          if (i > 1)

          {

              T[1][i] = h[i];

              T[2][i] = 2*(h[i-1] + h[i]);

         }

         if (i < n-1)

            T[3][i] = h[i];

      }

      T[2][1] = T[2][n] = 1;

      /* tridec is a routine specially optimized for solving tridiagonal systems */

      tridec (T,v,n,c);

      freemat (T,3);

      freevec (v);

      freevec (w);

      freevec (h);

   }  // end if(ic) only done once

 

   /* Evaluate spline function at a */

 

   while ((a > x[k+1]) && (k < n-1)) k++;  // find segment index k for given a

   p = x[k+1] - x[k];      // p is like h: the current delta x                                                                    

   b0 = (6*y[k+1] - p*p*c[k+1])/(6*p);  // the coefficients of the linear

   b1 = (6*y[k]   - p*p*c[k])/(6*p);   // terms in the Hultquist form for cubic

   s = b0*(a - x[k]) + b1*(x[k+1] - a);

   s += (float) (pow(a-x[k],3)*c[k+1] + pow(x[k+1]-a,3)*c[k])/(6*p);

   return s;

}

 

When you first run example program e442.m you have the options of appending the x-vector.  Then the x and y vector are initialized with this code:

 

for (i = 1; i <= n; i++)

   {

      a = -d + 2*d*(i-1)/(n-1);

      x[i] = a;

      y[i] = 1/(1 + a*a);

   }

 

which hard-wires the sampling the function

                                                                

at n points between d and d, according to the indexing formula:

                                                              

In this code, the distance from the center to the boundary of the input interval d=4 and the number of control points n=9 are the default values, so after pressing “Enter” to continue, the command window is filled with both the 9 integer values of x from -4 to 4 and corresponding values for y:

 

Example 4.4.2: Cubic Spline Fit

 

     -------------

    |          -4 |

    |          -3 |

    |          -2 |

    |          -1 |

x = |           0 |

    |           1 |

    |           2 |

    |           3 |

    |           4 |

     -------------

1 => Save, 2 => Append, Enter => continue:

 

     -------------

    |     0.05882 |

    |         0.1 |

    |         0.2 |

    |         0.5 |

y = |           1 |

    |         0.5 |

    |         0.2 |

    |         0.1 |

    |     0.05882 |

     -------------

1 => Save, 2 => Append, Enter => continue:

 

The showvec function has this save/append/continue option at it’s end since these displays result from the commands:

 

showvec ("x",x,n,0);

showvec ("y",y,n,0);

 

Having the x and y vectors defined, we are ready for the call

 

spline (0,x,y,c,n,1);

 

which primes the pump for computing the splines by evaluating the spline coefficients (for the cubic terms in the Hurlquist cube)  ci  for  0<i<9  This only needs doing once.

 

Now that the linear coefficients c have been computed, m = 10*n = 90 is the number of times the function will be sampled to draw the curve. Let’s see:

 

for (i = 1; i <= m; i++)

   {

      a = x[1] + (i-1)*(x[n] - x[1])/(m-1);

      X[1][i] = a;

      Y[1][i] = spline (a,x,y,c,n,0);

      j = MIN (i,n);

      X[2][i] = x[j];

      Y[2][i] = y[j];

   }

 

Let’s try to thread it through and hope that’ll show what’s going on.

 

Mindful of the boundary conditions x[1] = -4 and x[n] = 4, when i = 1, we have

      X[1][1] = a = -4  

and so

               Y[1][1] = spline (-4,x,y,c,9,0)

will return 1/(1+16) ~ 0.05882, right?  As for j, well, on the first pass through, i = 1 and n = 9, so the minimum is…(drum roll, please) j = MIN(1,9) = 1!, but i quickly ascends to 9 and continues onward to 90.  That means the MIN j value goes from 1 to 9 and stays there while i marches on to 90.  More about this in a minute…

 

Once the 2x90 X and Y matrices are set up, we’re ready for the final big function call

  

graphmm (X,Y,2,m,"Cubic Spline Fit","x","y","fig442",1);

 

Which is a NLIB function found in draw.c whose comment section explaining the meaning of the input parameters is show below.  Note that it seems like m has switched from the number of rows to the number of columns, and this is supported by the final comment “X and Y are transposed.”

 

/* -------------------------------------------------------------------- */

void graphmm (matrix X, matrix Y, int m, int n, char *title,

              char *xaxis, char *yaxis, char *fname, int dm)

/* -------------------------------------------------------------------- */

{

/* Description: Create a Matlab M-file which graphs the rows of the

                m by n matrix Y vs. the rows of the m by n matrix X

                and labels the graph.

On entry:    X     = m by n matrix of independent variables

             Y     = m by n matrix of dependent variables

             m     = number of rows of X and Y

             n     = number of columns of X and Y

             title = graph title

             xaxis = x axis label

             yaxis = y axis label

             fname = Matlab M-file name

             dm    = display mode. Normally, the points are

                     connected with solid straight line segments.

                     If dm != 0, then the last row of Y is

                     graphed using isolated plot symbols drawn

                     at each point.

 

   Notes:    1) X and Y are transposed in the M-file */

 

Aha!   The display mode parameter is set to dm = 1 (not 0) and that’s how the last (2nd) rows are using different (isolated) plot symbols!  So the broken line plot is X[1] v. Y[1] and the isolated points plot is X[2] v. Y[2] .   Part of what was hard to understand is how, after i grows larger than n on its way to 10n, the second rows of X and Y are all the same…no trouble of course, since it’s costs so little, why not?