Lagrange Interpolator Function

In astronomy, interpolation within tables is often a very useful procedure.  One very effective method of tabular interpolation is Lagrangian interpolation, covered here.  The Lagrangian Interpolator function was originally designed to interpolate within astronomical tabulations, but it can be used for interpolating within just about any continuous numerical data table.

We are given two columns of matching XY data pairs.  Given any general value of X within the range of a given table, the idea is to compute the corresponding value of Y.  In astronomy this could be a table of times vs some astronomical coordinate within which we want to interpolate for some specific non-tabulated value for a given moment.

For example, we may have a table of astronomical positions of the sun at intervals of 1 hour, but we may need to find the position of the sun at the time 16:25, which is not in the table of hourly values.  To do this, we interpolate, and this Lagrangian method is one of several methods we may apply, depending on the situation.  It is a very simple process to implement in just about any programming language, especially PHP.

Mathematically, the data is treated like a two-column table of X-data vs Y-data values.  There must be a corresponding Y-data point for every given X-data point.  The intervals between the XY data points do not have to be linear or uniform.  Lagrangian interpolation determines the polynomial that relates or fits all the given paired XY data points exactly.  Once the function determines this polynomial, we may then interpolate across the table so that, given any general X-argument, we may directly compute the corresponding value of Y from it.

Such an XY data table might look something like:
 X    Y
 19   23.9294430
 20   23.9902584
 21   24.0510412
 22   24.1117964
Given a value of X within the table, we interpolate the corresponding Y value.
 If X = 20.72, then Y = 24.03402501312


Given n data points consisting of n correlated XY data pairs at arbitrary XY intervals, Lagrange interpolation will define an (n-1) degree polynomial that passes through all n points.  We can use this polynomial to compute any general XY values in between the tabulated values.

Mathematically, it can be formally expressed as:


In the case of (j == i), we equate the inner product to 1

Internet references consulted while designing the following Lagrangian interpolator function:
http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html
http://en.wikipedia.org/wiki/Lagrange_polynomial



The following listing is a PHP function to perform instant Lagrangian interpolation.  Given an XY data table and an X-argument, it returns the corresponding interpolated Y value.  It is the foundation of the Lagrangian Interpolator Tool program.  Internally it uses arbitrary precision arithmetic.
/*
   Lagrangian interpolator function for any number of XY data pairs.

   Version  : 2.1
   Author   : Jay Tanner 2014
   Language : PHP v5.4.7
   License  : This source code is released under the provisions of the
              GNU Affero General Public License (AGPL), version 3.
              http://www.gnu.org/licenses/agpl.html
*/

 function Lagrange_Interpolate ($XYDataStr, $xArgStr, $nDecimals=16)

{

 $XDataStr = $YDataStr = "";

 $d = $nDecimals;

 $z5 = "0." . str_repeat("0", $d) . "5";

 $XY = preg_split("[ ]", preg_replace("/\s+/", " ", trim($XYDataStr)));

 $TotalDataCount = count($XY);

 $n = $TotalDataCount / 2;

 if ($n != floor($n + 0.5))
    {return "ERROR: XY Data Count Mismatch. Odd data element.";}

 if ($TotalDataCount < 4)
    {return "ERROR: There must be at least two XY data pairs.";}

 for($i=0;   $i < $TotalDataCount;   $i+=2)
    {
     $XDataStr .= $XY[$i]   . " ";
     $YDataStr .= $XY[$i+1] . " ";
    }

 $X = preg_split("[ ]", trim($XDataStr));
 $Y = preg_split("[ ]", trim($YDataStr));

 $x = trim($xArgStr);  if ($x == "") {$x = "0";}

 $y = 0;

 for ($i=0;   $i < $n;   $i++)
     {
      $Li = "1";

      for ($j=0;   $j < $n;   $j++)
          {
           if ($j != $i) // Skip this cycle if j == i
              {
               $u  = bcmul($Li, bcsub($x, $X[$j], 32), 32);
               $w  = bcsub($X[$i], $X[$j], 32);
               $Li = bcdiv($u, $w, 32);
              }
          }

      $y = bcadd($y, bcmul($Y[$i], $Li, 32), 32);
     }

 $s = ($y < 0)? -1:1;
 $y = bcmul($s, bcadd(bcmul($s, $y, 32), "$z5", $d), $d);

 if (strpos($y, ".") !== FALSE) {$y = rtrim(rtrim($y, "0"), ".");}

 return $y;

}
The following is a fully commented version of the above function.
/*
   Lagrangian interpolator function for any number of XY data pairs.

   Author   : Jay Tanner - 2014
   Language : PHP v5.4.7
   License  : This source code is released under the provisions of the
              GNU Affero General Public License (AGPL), version 3.
              http://www.gnu.org/licenses/agpl.html
*/

   function Lagrange_Interpolate ($XYDataStr, $xArgStr, $nDecimals=16)
{

// -----------
// Initialize.

   $XDataStr = $YDataStr = "";

// ------------------------------------------------------
// Create format string for rounding off the output based
// on the number (1 to 16) of decimals specified.

   $d = $nDecimals;

   $z5 = "0." . str_repeat("0", $d) . "5";

// --------------------------------------------------------------
// Read and split XY data pairs into a work array.  In the array,
// even number indices = X-data, odd number indices = Y-data.

   $XY = preg_split("[ ]", preg_replace("/\s+/", " ", trim($XYDataStr)));

// ---------------------------------------------
// Count the total number of data elements. This
// value should always be an even number.

   $TotalDataCount = count($XY);

// Number of data pairs.  This value should be an integer value
// exactly equal to 1/2 the total number of data points. If not,
// then there is an odd mismatched data point.

   $n = $TotalDataCount / 2;

// -----------------------------------------------------------
// Return error message if data vector element count mismatch.
// For every X value there must be a corresponding Y value or
// an XY data count mismatch error occurs.  An error will also
// occur if insufficient data.  There must be at least two XY
// data pairs given.

   if ($n != floor($n + 0.5))
      {return "ERROR: XY Data Count Mismatch. Odd data element.";}

   if ($TotalDataCount < 4)
      {return "ERROR: There must be at least two XY data pairs.";}

// -------------------------------------------------------
// Construct separate XY data strings from the array data.
// The XY data strings should each contain the same number
// of data elements.

   for($i=0;   $i < $TotalDataCount;   $i+=2)
      {
       $XDataStr .= $XY[$i]   . " ";
       $YDataStr .= $XY[$i+1] . " ";
      }

// --------------------------------------------------------------
// Split the created XY data vector strings into matching indexed
// arrays.  For every X value there must be a matching Y value
// and no two X values can be identical.

   $X = preg_split("[ ]", trim($XDataStr));
   $Y = preg_split("[ ]", trim($YDataStr));

// ----------------------------------------
// Read X argument for which to interpolate
// the Y value from the given XY data.

   $x = trim($xArgStr);  if ($x == "") {$x = "0";}

// -----------------------------------
// Initialize Y summation accumulator.
   $y = 0;

// -----------------------------------------------------
// Compute Lagrangian product (Li) for given X argument.

   for ($i=0;   $i < $n;   $i++)
       {
        $Li = "1";

        for ($j=0;   $j < $n;   $j++)
            {
             if ($j != $i) // Skip this cycle if j == i
                {
                 $u  = bcmul($Li, bcsub($x, $X[$j], 32), 32);
                 $w  = bcsub($X[$i], $X[$j], 32);
                 $Li = bcdiv($u, $w, 32);
                }
            } // Next j

//      --------------------------------------
//      Accumulate sum of Yi polynomial terms.
//      Y += Yi*Li

        $y = bcadd($y, bcmul($Y[$i], $Li, 32), 32);

       } // Next i

// ---------------------------------------------------------------
// Round off interpolated Y value to specified number of decimals.

   $s = ($y < 0)? -1:1;
   $y = bcmul($s, bcadd(bcmul($s, $y, 32), "$z5", $d), $d);

// --------------------------------------------------------------
// Cut off any redundant zeros and/or decimal point from Y value.
// This extra step is simply to tidy up the returned Y value.

   if (strpos($y, ".") !== FALSE) {$y = rtrim(rtrim($y, "0"), ".");}

// Done.
   return $y;

} // End of  Lagrange_Interpolate(...)
©Jay Tanner - 2017