javascriptpolynomial-math

How to calculate coefficients of polynomial using Lagrange interpolation


I need to calculate coefficients of polynomial using Lagrange interpolation polynomial, as my homework, I decide to do this in Javascript.

here is definition of Lagrange polynomial (L(x))

enter image description here

Lagrange basis polynomials are defined as follows

enter image description here

Calculate y value for specific X (W(x) function) is simple but I need to calculate coefficients of polynomial (array of [a0, a1, ..., an]) I need to do this to n<=10 but it will be nice to have arbitrary n, then I can put that function into horner function and draw that polynomial.

enter image description here

I have function that calculate denominator in first equation

function denominator(i, points) {
   var result = 1;
   var x_i = points[i].x;
   for (var j=points.length; j--;) {
      if (i != j) {
        result *= x_i - points[j].x;
      }
   }
   return result;
}

and function that return y using horner method (I also have drawing function using canvas)

function horner(array, x_scale, y_scale) {
   function recur(x, i, array) {
      if (i == 0) {
         return x*array[0];
      } else {
         return array[i] + x*recur(x, --i, array);
      }
   }
   return function(x) {
      return recur(x*x_scale, array.length-1, array)*y_scale;
   };
}

anybody know algorithm to do this, or idea how to calculate those coefficients


Solution

  • Well, you can do it the naive way. Represent a polynomial by the array of its coefficients, the array

    [a_0,a_1,...,a_n]
    

    corresponding to a_0 + a_1*X + ... + a_n*X^n. I'm no good with JavaScript, so pseudocode will have to do:

    interpolation_polynomial(i,points)
        coefficients = [1/denominator(i,points)]
        for k = 0 to points.length-1
            if k == i
                next k
            new_coefficients = [0,0,...,0] // length k+2 if k < i, k+1 if k > i
            if k < i
                m = k
            else
                m = k-1
            for j = m downto 0
                new_coefficients[j+1] += coefficients[j]
                new_coefficients[j] -= points[k]*coefficients[j]
            coefficients = new_coefficients
        return coefficients
    

    Start with the constant polynomial 1/((x_1-x_0)* ... *(x_i-x_{i-1})*(x_i-x_{i+1})*...*(x_i-x_n)) and multiply with X - x_k for all k != i. So that gives the coefficients for Li, then you just multiply them with yi (you could do that by initialising coefficients to y_i/denominator(i,points) if you pass the y-values as parameters) and add all the coefficients together finally.

    polynomial = [0,0,...,0] // points.length entries
    for i = 0 to points.length-1
        coefficients = interpolation_polynomial(i,points)
        for k = 0 to points.length-1
            polynomial[k] += y[i]*coefficients[k]
    

    Calculating each Li is O(n²), so the total calculation is O(n³).

    Update: In your jsFiddle, you had an error in the polynomial multiplication loop in addition to (the now corrected) mistake with the start index I made, it should be

    for (var j= (k < i) ? (k+1) : k; j--;) {
         new_coefficients[j+1] += coefficients[j];
         new_coefficients[j] -= points[k].x*coefficients[j];
    }
    

    Since you decrement j when testing, it needs to start one higher.

    That doesn't produce a correct interpolation yet, but it's at least more sensible than before.

    Also, in your horner function,

    function horner(array, x_scale, y_scale) {
       function recur(x, i, array) {
          if (i == 0) {
             return x*array[0];
          } else {
             return array[i] + x*recur(x, --i, array);
          }
       }
       return function(x) {
          return recur(x*x_scale, array.length-1, array)*y_scale;
       };
    }
    

    you multiply the highest coefficient twice with x, it should be

    if (i == 0) {
        return array[0];
    }
    

    instead. Still no good result, though.

    Update2: Final typo fixes, the following works:

    function horner(array, x_scale, y_scale) {
       function recur(x, i, array) {
          if (i == 0) {
             return array[0];
          } else {
             return array[i] + x*recur(x, --i, array);
          }
       }
       return function(x) {
          return recur(x*x_scale, array.length-1, array)*y_scale;
       };
    }
    
    // initialize array
    function zeros(n) {
       var array = new Array(n);
       for (var i=n; i--;) {
         array[i] = 0;
       }
       return array;
    }
    
    function denominator(i, points) {
       var result = 1;
       var x_i = points[i].x;
       for (var j=points.length; j--;) {
          if (i != j) {
            result *= x_i - points[j].x;
          }
       }
        console.log(result);
       return result;
    }
    
    // calculate coefficients for Li polynomial
    function interpolation_polynomial(i, points) {
       var coefficients = zeros(points.length);
        // alert("Denominator " + i + ": " + denominator(i,points));
       coefficients[0] = 1/denominator(i,points);
        console.log(coefficients[0]);
        //new Array(points.length);
       /*for (var s=points.length; s--;) {
          coefficients[s] = 1/denominator(i,points);
       }*/
       var new_coefficients;
    
       for (var k = 0; k<points.length; k++) {
          if (k == i) {
            continue;
          }
          new_coefficients = zeros(points.length);
           for (var j= (k < i) ? k+1 : k; j--;) {
             new_coefficients[j+1] += coefficients[j];
             new_coefficients[j] -= points[k].x*coefficients[j];
          }   
          coefficients = new_coefficients;
       }
       console.log(coefficients);
       return coefficients;
    }
    
    // calculate coefficients of polynomial
    function Lagrange(points) {
       var polynomial = zeros(points.length);
       var coefficients;
       for (var i=0; i<points.length; ++i) {
         coefficients = interpolation_polynomial(i, points);
         //console.log(coefficients);
         for (var k=0; k<points.length; ++k) {
           // console.log(points[k].y*coefficients[k]);
            polynomial[k] += points[i].y*coefficients[k];
         }
       }
       return polynomial;
    }