User:Rcampbell1/sandbox
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 j≠i, 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)