domingo, 23 de abril de 2017

Interpolación por Polinomios de Newton y Lagrange


Justificación y propósito
La interpolación es utilizada para hacer predicciones de un valor que se encuentra entre una serie de puntos obtenidos de observaciones experimentales o empíricas. En otras palabras, valores que no pudieron ser medidos u observados en alguna clase de experimento.
El algoritmo de interpolación encuentra una función que cruza todos los puntos dados, y de esta manera, es posible evaluar cualquier valor dentro del valor mínimo y el máximo de dichos puntos para estimar cualquier otro valor. Sin embargo, no significa que no tenga error, de hecho, la función creada no describe el comportamiento real de la información utilizada para la interpolación. Por ejemplo, si se mide la temperatura de cierta región en un lapso de tiempo y se crea una función que, para estos datos, dicha función podrá describir de manera aproximada los datos intermedios, pero los datos fuera del rango pueden no tener sentido o tener un incremento o decremento descontrolado.
Volviendo a algo más específico, la razón por la que se utiliza la interpolación por Newton o Lagrange, es porque estos pueden crear una función con n puntos y, a diferencia de una línea recta (interpolación linear) o una simple parábola (interpolación cuadrática), tiende a generar una mejor precisión en las predicciones de valores intermedios.

Antecedentes
Imaginemos que tenemos 2 puntos:
Podemos dibujar una línea recta entre ellos:
Y elegir un punto arbitrario entre estos puntos:
Así, podemos calcular varias distancias:
photo_2017-04-20_19-29-10.jpg
Entonces (tomando en cuenta los 2 triángulos que se forman), podemos hacer la siguiente igualdad:
Y resolviendo para f(x):


 Explicación gráfica de la interpolación
La interpolación consiste en encontrar una función que cruce por una cierta cantidad de puntos.

En el caso de la interpolación por polinomios de Newton y Lagrange, se puede utilizar cualquier cantidad de puntos. Esto, debido a que estos métodos buscan un polinomio cuyo grado depende de la cantidad de valores conocidos, en otras palabras, utiliza una fórmula generalizada que se puede aplicar sin importar la cantidad de puntos y genera un polinomio que se ajusta a dichos puntos.

El funcionamiento gráfico es exactamente el mismo, una vez que se obtiene el polinomio, este ya se puede utilizar para evaluar el valor intermedio deseado.

Como se muestra en la siguiente imagen, con tres puntos se creó un polinomio de tercer grado que se ajusta a los puntos. Posteriormente, se evalúa en un punto intermedio y se da un valor cercano al real.



Codigo en C++

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
#include <iostream>
#include <cmath>
#include <vector>

using namespace std;
double EPSILON = 0.00001;       //Error


void interpolacionLineal( double y1, double y2, double x1, double x2, double x ) {

    double y = 0.0;

    y = y1 + ( y2 - y1 ) * ( x - x1 ) / ( x2 - x1 );

    cout << "Valor interpolado es (y): " << y << endl;
}

void interpolacionCuadrtica (  double y0, double y1, double y2, double x0, double x1, double x2 ) {

    double b0 = y0;

    // Diferencia dividida discreta
    double b1 = ( y1 - y0 ) / ( x1 - x0 );

    double b2 = ( ( ( y2 - y1 ) / ( x2 - x1 ) ) - b1 ) / ( x2 - x0);

    double a = b2;
    double b = ( b1 - b2 * x1 - b2 * x0 );
    double c = ( b0 - b1 * x0 + b2 * x0 * x1 );

    cout << "\n++++++++++++++++++++++++++++++" << endl;
    cout << "b0: " << b0 << "\tb1: " << b1 << "\tb2: " << b2 << endl;
    cout << "Funcion cuadratica: " << a << "x^2 + " << b << "x + " << c << endl << endl;
}

/**
    Metodo utilizado en la interpolacion Polinomial de Newton
*/
double divisionDiscreta( double *y, double *x, int n, int i, int j, vector<double> *B ) {

    if ( n == 1 ) {

        if ( i == 0 ) {

            B->push_back( y[i] );
            cout << "f(x0): " << y[i] << endl;
        }

        return y[i];
    } else {
        //cout << "xi: " << i << endl;
        //cout << "xj: " << j << endl;
        double val = ( divisionDiscreta( y, x, n - 1, i, j + 1, B ) - divisionDiscreta( y, x, n - 1, i - 1, j, B ) ) / ( x[i] - x[j] );
        if ( j == 0 ) {
            B->push_back( val );
        }
        cout << "f[ x" << i << ", x" << j << " ]: " << val << endl;
        return val;
    }
}

void interpolacionPolinomialNewton() {
    /* Nota: La suma de los binomios son conmutativas */
    double y[] =
    //{2.31, 3.36, 4.59, 6};
    //{0, 1.386294, 1.791759, 1.609438};
    {0.1353, 7.3890, 20.0855};

    double x[] =
    //{0.1, 0.4, 0.7, 1};
    //{1, 4, 6, 5};
    {-2, 2, 3};

    double x_inc = -1;

    vector<double> B;

    int n = sizeof( y ) / sizeof( y[0] );

    divisionDiscreta( y, x, n, n - 1, 0, &B );

    //
    int j = B.size();
    double v = 0;

    cout << "Polinomio: ";                      //Print Function
    while ( j --> 0 ) {

        cout << B[j];                           //Print Function

        double temp = B[j];
        int k = j;
        while ( k --> 0 ) {

            cout << "(x - " << x[k] << ")";     //Print Function
            temp *= ( x_inc - x[k] );
        }

        if ( j != 0)                            //Print Function
            cout << " + ";                      //Print Function

        v += temp;
    }

    cout << "\nf(x) = " << v << endl;

    // El polinomio es: b0 + b1(x - x0) + b2(x - x0)(x - x1) + b3 (x - x0)(x - x1)(x - x2)
}

/**
    Formula utilizda en la iterpolacion de Lagrange
*/
double L( int i, double x_incognita, double *x, int n ) {

    int j = 0;
    double val = 1;

    while ( j < n ) {

        if ( i != j ) {

            val *= ( x_incognita - x[j] ) / ( x[i] - x[j] );
        }

        j++;
    }

    return val;
}

double f( double *y, double *x, double x_incognita, int n ) {

    int i = 0;

    double val = 0;

    while ( i < n ) {

        val += L( i, x_incognita, x, n ) * y[i];
        i++;
    }

    return val;
}

void InterpolacionLagrange() {

    double y[] =
    //{2.31, 3.36, 4.59, 6};
    {0, 1.386294, 1.791759, 1.609438};

    double x[] =
    //{0.1, 0.4, 0.7, 1};
    {1, 4, 6, 5};

    int n = sizeof( y ) / sizeof( y[0] );
    double x_incognita = 10;

    cout << "y calculada: " << f( y, x, x_incognita, n ) << endl;;
}

int main() {
    double y0, y1, y2, x0, x1, x2;

    // y1 = 86.4;
    // y2 = 92.5;

    // x1 = 1.0;
    // x2 = 1.5;

    // x = 1.25;

    // interpolacionLineal ( y1, y2, x1, x2, x );

    y0 = 0;
    y1 = 1.386294;
    y2 = 1.791759;

    x0 = 1;
    x1 = 4;
    x2 = 6;

    //interpolacionCuadrtica( y0, y1, y2, x0, x1, x2 );

    interpolacionPolinomialNewton();

    //InterpolacionLagrange();

    return 0;
}
Ejemplos

Input
Valor x deseado y f(x) verdadero
Output
f(x) = ln(x)

x
y
1
0
4
1.386294
5
1.791759
6
1.609438

X = 2
f(2) = 0.69315





X = 3
f(3) = 1.09861
Newton:
X = 2
f(x) = 0.628767

X = 3
f(x) = 1.07513
Lagrange:
X = 2
f(x) = 0.628767

X = 3
f(x) = 1.07513
f(x) = log x

x
y
8
0.90309
9
0.9542425
11
1.0413927

X = 10
f(x) = 1
Newton:
X = 10
f(x) = 1.00034

Lagrange:
X = 10
f(x) = 1.00034
f(x) = e^x

x
y
-2
0.1353
2
7.3890
3
20.0855


X = -1
f(x) = 0.3678

X = 0
f(x) = 1

X = 1
f(x) = 2.7182
Newton:
X = -1
f(x) = -4.58112

X = 0
f(x) = -4.94431

X = 1
f(x) = -0.95427

Lagrange:
X = -1
f(x) = -4.58112

X = 0
f(x) = -4.94431

X = 1
f(x) = -0.95427
Publicado por en

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)