User:Rcampbell1/sandbox

From Wikibooks, open books for an open world
Jump to navigation Jump to search

Lagrange interpolation is an algorithm which returns the polynomial of minimum degree which passes through a given set of points (xi, yi).

Algorithm[edit | edit source]

Given the n points (x0, y0), ..., (xn-1, yn-1), compute the Lagrange polynomial . Note that the ith term in the sum, is constructed so that when xj is substituted for x to have a value of zero whenever ji, and a value of yj whenever j = i. The resulting Lagrange polynomial is the sum of these terms, so has a value of p(xj) = 0 + 0 + ... + yj + ... + 0 = yj for each of the specified points (xj, yj).

In both the pseudocode and each implementation below, the polynomial p(x) = a0 + a1x + a2x2 + ... + an-1xn-1 is represented as an array of it's coefficients, (a0, a1, a2, ..., an-1).

Pseudocode[edit | edit source]

algorithm lagrange-interpolate is
    input: points (x0, y0), ..., (xn-1, yn-1) 
    output: Polynomial p such that p(x) passes through the input points and is of minimal degree

    for each point (xi, yi) do
        compute tmp := 
        compute term := tmp*          

    return p, the sum of the values of term

In sample implementations below, the polynomial p(x) = a0 + a1x + a2x2 + ... + an-1xn-1 is represented as an array of it's coefficients, (a0, a1, a2, ..., an-1).

While the code is written to expect points taken from the real numbers (aka floating point), returning a polynomial with coefficients in the reals, this basic algorithm can be adapted to work with inputs and polynomial coefficients from any field, such as the complex numbers, integers mod a prime or finite fields.

C[edit | edit source]

#include <stdio.h>
#include <stdlib.h>

void interpolate(int numpts, float* xval, float* yval, float* thepoly){
   float theterm[10];
   float prod;
   int i,j,k;
   for (i=0; i<numpts; i++){thepoly[i] = 0.0;};
   for (i=0; i<numpts; i++){
      prod = 1.0;
      for (j=0; j<numpts; j++){theterm[j] = 0.0;};
      // Compute Prod_{j != i} (x_i - x_j)
      for (j=0; j<numpts; j++){
         if (i==j) continue;
	 prod *= (xval[i] - xval[j]);
      };
      // Compute y_i/Prod_{j != i} (x_i - x_j)
      prod = yval[i]/prod;
      theterm[0] = prod;
      // Compute theterm := prod*Prod_{j != i} (x - x_j)
      for (j=0; j<numpts; j++){
	 if (i==j) continue;
	 for (k=numpts-1; k>0; k--){
            theterm[k] += theterm[k-1];
            theterm[k-1] *= (-xval[j]);
	 };
      };
      // thepoly += theterm (as coeff vectors)
      for (j=0; j<numpts; j++){thepoly[j] += theterm[j];};
   };
}
         
int main(){
   float xval[10], yval[10];
   float thepoly[10], tmpx, tmpy;
   int numpts, i;
   printf("number of points: "); scanf("%d",&numpts);
   printf("Enter %d points (x,y):\n",numpts);
   for (i=0; i<numpts; i++){
      printf("x y: "); scanf("%f %f",&tmpx,&tmpy);
      xval[i] = tmpx; yval[i] = tmpy;
   };
   interpolate(numpts,xval,yval,thepoly);
   printf("p(x) = (%6.2f)",thepoly[0]);
   if (thepoly[1] != 0.0) {printf(" + (%6.2f)x",thepoly[1]);};
   for (i=2; i<numpts; i++){
      if (thepoly[i] != 0.0){
         printf(" + (%6.2f)x^%d",thepoly[i],i);
      };
   };
   printf("\n");
   exit(0);
}

Python[edit | edit source]

def interpolate(inpts):
    n = len(inpts)
    thepoly = n*[0]
    for i in range(len(inpts)):
        prod = 1
        # Compute Prod_{j != i} (x_i - x_j)
        for j in (j for j in range(n) if (j != i)):
            prod *= (inpts[i][0] - inpts[j][0])
        # Compute y_i/Prod_{j != i} (x_i - x_j)
        prod = inpts[i][1]/prod
        theterm = [prod] + (n-1)*[0]
        # Compute theterm := prod*Prod_{j != i} (x - x_j)
        for j in (j for j in range(n) if (j != i)):
            for k in range(n-1,0,-1):
                theterm[k] += theterm[k-1]
                theterm[k-1] *= (-inpts[j][0])
        # thepoly += theterm
        for j in range(n):
            thepoly[j] += theterm[j]
    return(thepoly)