/*
Largrange Interpolation Algorithm
(1) give n, give x[0..n], give y[0..n]
(2) Ln,k(x) = (x-x0)*(x-x1)*...*(x-xn) / (xk-x0)*(xk-x1)-...*(xk-xn)
(3) P(x) = Ln1(x)f(x1)+Ln2(x)f(x2)+...+Lnn(x)f(xn)
assume function = -4x^4 + 2x^3-3x^2+4x-1
-by EdisonX
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define N 4
double func(double x)
{
double sum=0.0;
double x2=x*x, x3=x2*x,x4=x3*x;
sum = -4*x4+2*x3-3*x2+4*x-1;
return sum;
}
double LargrangeInterpolation(double x,double *xn, double *yn, unsigned n)
{
double numerator, denomenator; // 分子,分母
double P; // result
unsigned i, j;
for(i=0, P=0.0; i<n; i++){
if(yn[i]==0.0) continue;
numerator=1.0, denomenator=1.0;
for(j=0; j<n&& numerator!=0.0; j++){
if(i!=j) {
numerator*=(x-xn[j]);
denomenator*=(xn[i]-xn[j]);
}
}
P += numerator*yn[i]/denomenator ;
}
return P;
}
int main()
{
unsigned i;
double x=1.5;
double xn[N+1]={1.234, 3.782, 4.563, 5.987, -1.157};
double yn[N+1];
for(i=0; i<N+1; i++) yn[i]=func(xn[i]); // reference point
printf("correct: %lf\n", func(x));
printf("Lagrange: %lf\n",LargrangeInterpolation(x, xn, yn, N+1));
return 0;
}
output :
correct: -15.250000
Lagrange: -15.250000
留言列表