Polinomio de interpolación de Lagrange: C++ Linux

Si tenemos n pares ordenados correspondientes a Y como una función de X, siempre es posible unirlos mediante una función polinómica a la cual denominaremos polinomio de interpolación. El procedimiento numérico de Lagrange, que da lugar a un polinomio de grado n-1, es el que nos ocupa para el siguiente código:

#include <iostream>
#include <fstream>

using namespace std;

#include <iomanip>

void monomio2(int n, double *x, double *D);

int main () {
	
	cout << "Polinomio de interpolacion de Lagrange\n\n";

	ifstream label1 ("datos//datos.in");
	
	int n, i, j, k;

	label1 >> n;

	double *x, *xp, *f, *L, *D, producto, sum, *P;

	x = new double [n], xp = new double [n], f = new double [n];
	L = new double [n], D = new double [n], P = new double [n];

	cout << "Numero de pares ordenados = " << n << "\n\n";

	cout << "Valores de x y f(x)\nx   f(x)\n";

	for (i= 0; i < n; i++){

		label1 >> x[i] >> f[i];
		cout << x[i] << setw(5) << setiosflags(ios::right) << f[i] << "\n";

	}

	cout << endl;

	for (k = 0; k < n; k++){

	producto = 1;

		for (i = 0; i < n; i++){

			if (i != k) producto = producto * (x[k] - x[i]);

		}
	
	L[k] = f[k]/producto;

	}

	cout << "Coeficientes de interpolacion de Lagrange\n\n";

	cout.setf(ios::fixed);

	cout.precision(1);


	for (i= 0; i < n; i++){

		cout << "L(" << n-i-1 << ") = " << setw(4) << setiosflags(ios::right) << L[i] << endl;

	}

	cout << endl;

	for (k = 0; k < n; k++){
	
		j = 0;	
		
		for (i = 0; i < n; i++){

			if (i != k) {xp[j] = x[i]; j += 1;}

		}
	
		monomio2 (n-1, xp, D);

		for (i= 0; i < n; i++){

			P[i] = P[i] + L[k]*D[i];

		}
	 
	}

	cout << "Polinomio de interpolacion de Lagrange (grado " << n-1 << ")\n\n";

	for (i= 0; i < n; i++){

		cout << "P(" << n-i-1 << ") = " << setw(4) << setiosflags(ios::right) << P[i] << endl;
			
	}

	return 0;

}

void monomio2(int n, double *x, double *D){

	double *E;

	E = new double [n];

	D[0] = 1;

	D[1] = -x[0];

	for (int i = 1; i < n; i++) {

		for (int k =1; k < i+1; k++) {

			E[k] = D[k] + D[k-1]*(-x[i]);

		}

		D[i+1] = D[i]*(-x[i]);

		for (int j = 1; j < i+1; j++) {

			D[j] = E[j];

		}

	}

	delete E;

}

Lo emplearemos para encontrar el polinomio de interpolación de Lagrange (de grado 2) a partir de tres pares ordenados (x, f(x)) que generé usando la ecuación de la parábola:

f(x) = x² + 2x -5

Los valores utilizados fueron los siguientes:

x  y
0 -5
1 -2
2  3  

Estos valores, conjuntamente con el número de pares ordenados (colocado al inicio), se incluyeron en el archivo datos/datos.in y se introducen desde el programa principal. El algoritmo usa la función monomio2, que incluye unas ligeras modificaciones en la ya conocida monomio, para permitir en este caso determinar los desarrollos de (x-x1).( x-x2) … (x-xn), en lugar de (x+x1).( x+x2) … (x+xn), para cada i diferente de k. Después de copiar el código y grabarlo como lagrange.c++, abrimos una cónsola y nos vamos hacia la carpeta LAGRANGE y lo compilamos con:

g++ lagrange.c++ -o lagrange

La ejecución del mismo (./lagrange [Enter]) produce esta salida:

Polinomio de interpolacion de Lagrange

Numero de pares ordenados = 3

Valores de x y f(x)
x f(x)
0 -5
1 -2
2 3

Coeficientes de interpolacion de Lagrange

L(2) = -2.5
L(1) = 2.0
L(0) = 1.5

Polinomio de interpolacion de Lagrange (grado 2)

P(2) = 1.0
P(1) = 2.0
P(0) = -5.0

donde se verifica que los coeficientes del polinomio cuadrático obtenido ajustan perfectamente a la ecuación f(x) = x² + 2x -5 ya referida. Sin embargo, no nos emocionemos tanto por este éxito. Si incluimos, por ejemplo, 10 pares ordenados para los valores experimentales de una variable como función de otra, el polinomio que obtendríamos sería de grado 9 y, aunque pasará por todos los puntos considerados, la gráfica que obtendremos posiblemente se comportará de una manera no deseada fuera de ellos. Existen técnicas de interpolación más poderosas como la de los splines pero, este caso, es válido e interesante como ejercicio de programación.