Interpolación y Aproximación por polinomios

Los temas cubiertos en este capítulo están basados en el libro
Numerical Analysis, de Richard L. Burden y J. Douglas Faires, PWS-Kent, 3rd Edition (1985)

El Teorema de aproximación de Weierstrass nos asegura que dada una función definida en un intervalo cerrado, podemos encontrar un polinomio que "se acerca" a esta función tanto como lo deseemos. Veremos en este capítulo diversas formas de construir estos polinomios, y sus respectivas ventajas/desventajas.

Polinomios de Taylor

Es el polinomio de grado n que pasa por un punto dado, y tiene las n derivadas iguales a la función en ese punto. El polinomio de Taylor para una función f en x0 es: 


Pn(x) = f(x0) + f'(x0)(x - x0) +(1/2!) f''(x0)(x - x0)2 + ... + (1/n!) f(n)(x0)(x - x0)n

con un error conocido (resto Rn) igual a: 


f(x) - Pn(x) = Rn(x) = (1/(n+1)!) f(n+1)(ξ)(x - x0)n+1

Ejercicios:

  • Calcular P3(x) - el Polinomio de Taylor de orden 3 - alrededor de x0=0, para f(x) = (1+x)1/2
  • Usar este polinomio para aproximar (1.1)1/2. Comparar los errores de las sucesivas aproximaciones P1, P2 y P3
  • Usar este polinomio para integrar f(x) = (1+x)1/2 entre 0 y 0.1. Encontrar el error de esta aproximación
  • Calcular Pn(3) (n=1,7) alrededor de x0=1 para aproximar una función f(x) = 1/x .
  • Graficar P1(x), P2(x) y P3(x), entre x=0 y x=4.

    Polinomios de Lagrange

    Interpolar por este método consiste en construir un polinomio Pn(x) de grado máximo n, que pasa por n+1 puntos dados. La función f(x) tabulada en n+1 puntos xk es: 

    f(xk) = Pn(xk) para todo k=0,1,2,...,n


    El polinomio Pn(x) está dado por una combinación de los respectivos Polinomios de Lagrange Ln,k(x) 

    
Pn(x) = f(x0)Ln,0(x) + f(x1)Ln,1(x) + ... + f(xn)Ln,n(x)


    y los polinomios de Lagrange se definen como: 

    
          (x-x0)(x-x1)....(x-xk-1)(x-xk+1)....(x-xn)
Ln,k(x) = -------------------------------------------
	  (xk-x0)(xk-x1)..(xk-xk-1)(xk-xk+1)..(xk-xn)


    Como ejemplo vemos los 5 polinomios de Lagrange que pasan por los puntos x= 1.0, 1.3, 1.6, 1.9 y 2.2:


    El error dado por éste método de interpolación es: 

    
f(x) - Pn(x) = Rn(x) = (1/(n+1)!) f(n+1)(ξ)(x - x0)(x - x1) ... (x - xn)

    Ejercicios:

  • Calcular f(x) = 1/x utilizando una interpolación de Lagrange de segundo orden, para x=3. La función se da tabulada en x1=2, x2=2.5, y x3=4.
  • Se prepara una tabla con valores de ex, entre 0 y 1, con d decimales. (por ejemplo: si d=5, f(1)=2.71828). Los valores de x de la tabla están separados por un paso h.
  • Utilizar el programa lagrange.for para evaluar el valor de J0(x=1.5) (función de Bessel del primer tipo y orden cero). La función está tabulada en 5 puntos, en el file table.dat. Calcular los valores de Pn(x=1.5) utilizando interpolación de Lagrange, en los diferentes ordenes, y con las diferentes combinaciones posibles de puntos adyacentes.
  • Este programa imprime tambien los polinomios de Lagrange. Comparar graficamente el ejemplo de f(x) = 1/x, interpolado con los diferentes métodos vistos.
  • Agregar la estimación del error en el programa lagrange.for.

    Métodos de Neville y Aitken
    Una de las dificultades de la interpolación de Lagrange, es que el error es difícil (o imposible) de calcular. La forma habitual de trabajar es ir incrementando el orden de los polinomios, hasta que se obtiene un valor deseado. Sin embargo, cada cálculo es independiente del previo, perdiendose contacto entre uno y otro. Los polinomios de Legendre tambien se pueden generar aprovechando los cálculos previos, en forma iterativa.
    Calcularemos, usando los valores dados en table.dat, los polinomios de Lagrange de distinto orden y con distinta combinaciones de puntos adyacentes, para x=1.5. Llamaremos Pi,j,k(x) al polinomio de Lagrange de orden 2, que pasa por los puntos adyacentes x=xi, x=xj y x=xk.
    Se puede demostrar que el polinomio de Lagrange que pasa por los puntos adyacentes x=xi, x=xj, x=xk y x=xl se obtiene haciendo: 

    
                (x-xi)Pjkl - (x-xl)Pijk
P(x) = Pijkl(x) = ---------------------------
	                    (xl-xk)


    Por ejemplo: 

    
             (r-x1)P2 - (r-x2)P1
P12(r) = -----------------------------
                  ( x2 - x1 )
    en el caso de r=1.5 :
    
               (1.5 - 1.0)0.6200860 - (1.5 - 1.3)0.7651977
P12(1.5) =    ---------------------------------------------- = 0.5233449
                                ( 1.3 - 1.0 )
    o
    
               (1.5 - 1.3)0.4554022 - (1.5 - 1.6)0.6200860
P23(1.5) =    ---------------------------------------------- = 0.5102968
                                ( 1.6 - 1.3 )
    y en orden mas elevado: 
    
               (1.5 - 1.0)0.5102968 - (1.5 - 1.6)0.5233449
P123(1.5) =    ---------------------------------------------- = 0.5124715
                                ( 1.6 - 1.0 )


    Escribiendo los polinomios Pijkl en una matriz A

    x1A11
    x2A21A22
    x3A31A32A33
    x4A41A42A43A44
    x5A51A52A53A54A55

    el algoritmo de Neville se puede simplificar en el siguiente esquema: 

    Interpolación de la función f(x) dada en n puntos x1, x2, ... ,xn
 
INPUT: xi en el array X(i) y f(xi) en la matriz A(i,1)

    •  Paso 1: For i=2,...,n
 		 for j=2,...,i
 		                    (x-xi-j+1)Ai,j-1 - (x-xi)Ai-1,j-1
		           Ai,j = ------------------------------------
		                             xi - xi-j+1

    •  Paso 2: Output A


    Ejercicios:

  • Utilizar el programa neville.for para llenar la siguiente tabla:

    x1P1
    x2P2P12
    x3P3P23P123
    x4P4P34P234P1234
    x5P5P45P345P2345P12345
    para los datos de J0 dados en table.dat, y para x=1.5. Comprobar si el programa funciona correctamente, revisando los resultados obtenidos en neville.dat.

  • Comparar el valor de P12345(x=1.5) con el obtenido interpolando sin iteración.
  • Si se sabe que J0(x=2.5)=-0.0483838, calcular el nuevo valor de J0(x=1.5).
  • Modificar el programa neville.for para interpolar la misma función utilizando el Método de Aitken. Este consiste en escribir el polinomio parcial como combinación de los elementos de la columna previa. Mientras que en Neville se utilizan el de la misma fila y la fila anterior, aquí se utiliza el de la misma fila, y el elemento de la primer fila.
    En notación matricial sería: 

              (x-xj-1)Ai,j-1 - (x-xi)Aj-1,j-1
Ai,j = ------------------------------------
                    xi - xj-1

  • Aproximar f(x)=31/2 con la función f(x)=3x tabulada en los puntos x1=-2, x2=-1, x3=0, x4=1 y x5=2. Utilizar los métodos de Neville y Aitken, y comparar los resultados.

    Método de diferencias divididas
    Los métodos iterativos explicados anteriormente tienen la desventaja que no escriben explícitamente el polinomio obtenido. Existen métodos que sí lo hacen, y se denominan métodos de diferencias divididas. Estos son muy útiles para desarrollar expresiones de derivadas e integrales.
    En este método, el polinomio de interpolación se define como 

    
Pn = f[x0] + Sumk=1n f[x0,x1,...,xk](x-x0)...(x-xk-1)

    donde los funcionales f[] se definen en forma recursiva: 

    
f[xi] = f(xi)
		    


                  ( f[xi+1] - f[xi] )
f[xi,xi+1] = --------------------------
                     ( xi+1 - xi )
		    


                              ( f[xi+1,xi+2, ... ,xi+k] - f[xi,xi+1, ... ,xi+k-1] )
f[xi,xi+1, ... ,xi+k-1,xi+k] = -----------------------------------------------------
                                                ( xi+k - xi )

    Por ejemplo, volviendo a nuestro ejemplo de la función J0, las diferencias finitas se expresan como:

    0.7651977
    -0.4837057
    0.6200860-0.1087339
    -0.54894606.5878395E-02
    0.4554022-4.9443333E-021.8251029E-03
    -0.57861206.8068519E-02
    0.28181861.1818333E-02
    -0.5715210
    0.1103623

    siendo el polinomio:

P4(x) = 0.7651977 - 0.4837057(x-1.0) - 0.1087339(x-1.0)(x-1.3) + 
                      0.0658784(x-1.0)(x-1.3)(x-1.6) + 0.0018251(x-1.0)(x-1.3)(x-1.6)(x-1.9)

    Ejercicios:

  • Escribir un programa que calcule J0, dados los datos en table.dat. Comprobar el resultado para x=1.5.
    Si no sale, se puede ver como ayuda el programa divideddiff.for.

    Polinomios de Hermite
    Los polinomios osculadores son una generalización de los polinomios de Taylor y de Lagrange. Estos interpolan la función dada, coincidiendo con ella en n+1 puntos y en sus m derivadas.
    Un caso particular son los polinomios de Hermite, que interpola la función dada, y coincide con ella en n+1, y en n puntos de la derivada primera. El polinomio de Hermite está dado por 

    
H2n+1(x) = Sumj=0n f(xj)Hn,j(x) + Sumj=0n f'(xj)HHn,j(x)

donde

Hn,j(x) = [1 - 2(x-xj)L'n,j(xj)]L2n,j(x)

y

HHn,j(x) = (x-xj)L2n,j(x)

y los Ln,j son los polinomios de Lagrange.


    Ejercicios:

  • Dibujar esquemáticamente los polinomios de Hermite (los H y los HH)



    Si bien la descripción anterior es completa, el hecho de tener que evaluar los polinomios de Lagrange y sus derivadas, lo hace un poco tedioso.
    Una forma simple de encontrar los coeficientes es utilizando diferencias finitas, pero definiendo nuevos puntos zi en la forma: 

    
z1 = z2 = x1
z3 = z4 = x2
y en general:
z2i = z2i-1 = xi


    Se hace el cálculo de diferencias finitas explicado anteriormente, pero como f[z2i,z2i-1]=f[xi,xi] y este no está definido, entonces se usa la expresión del límite 

    f[z2i,z2i-1]=f'(xi)



    Ejercicios:

  • Modificar el programa divideddiff.for para calcular los polinomios de Hermite.
  • Plotear la función y su derivada para nuestro viejo conocido caso de J0

    Splines

     Las interpolaciones con polinomios sufren de un problema básico, y es la aparición de grandes oscilaciones espúreas, especialmente si el grado del polinomio es alto.
    Una forma alternativa de obtener funciones interpoladoras, es dividiendo el intervalo en sucesivos intervalos, y generar polinomios de bajo grado en cada uno de estos intervalos. Esto se llama aproximación polinomial por piezas.
    Para que el resultado quede "bien", estos polinomios deben mantener continuidad en la función, la derivada primera y la derivada segunda, en cada uno de estos puntos límites.
    Dada una función f definida en [a,b], y un grupo de nodos a = x0 < x1 < ... < x n = b , se llama cubic spline interpolant S, a la función que satisface las siguientes condiciones para j<(n+1): 

    (a) S  es un polinomio cúbico, que en cada subintervalo [xj,xj+1] se denomina Sj
(b) S(xj) = f(xj) (para todo j)
(c) Sj+1(xj+1) = Sj(xj+1)
(d) S'j+1(xj+1) = S'j(xj+1)
(e) S''j+1(xj+1) = S''j(xj+1)
(f) S''(x0) = S''(xn) = 0 ("free boundary")
             o
   S'(x0) = f'(x0)  y S'(xn) = f'(xn) ("Clamped boundary")


    El procedimiento para encontrar S es el siguiente:
    denominamos 

    
Sj(x) = aj + bj(x - xj)  + cj(x - xj)2  + dj(x - xj)3


    y aplicamos las condiciones pedidas anteriormente. Por ejemplo, para que se cumpla (b): 

    
Sj(xj) = fj(xj) = aj


    para que se cumpla (c): 

    
aj+1 = Sj+1(xj+1) = Sj(xj+1) = aj + bj(x - xj+1)  + cj(x - xj+1)2  + dj(x - xj+1)3 

que para simplificar, haciendo hj = (xj+1 - xj) es:

aj+1 =  aj + bjhj  + cjhj2  + djhj3


    Respecto a la derivada, que se escribe 

    
S'j(x) = bj + 2cj(x - xj) + 3dj(x - xj)2


    y aplicando la condición (d) obtenemos: 

    
bj+1(x) = bj + 2cjhj + 3djhj2


    la condición (e) se cumple si 

    cj+1(x) = cj + 3djhj


    (para que se cumpla también en el último punto n se define cn = S''(xj)/2).
    Subsituyendo todas estas condiciones, nos queda un sistema de ecuaciones (matriz tridiagonal): 

    
hj-1cj-1 +  2(hj-1 + hj)cj + hjcj+1 =
                        = 3(aj+1 - aj)/hj - 3(aj - aj-1)/hj-1


    donde j=2,3, ... ,n-2. En las filas j=0 y j=n, se escriben las condiciones de contorno en la diagonal.

    Ejercicios:

  • Utilizar el programa splines.for para interpolar la siguiente figura.

    Debido a que las derivadas no son contínuas en algunos puntos, se interpolan 3 funciones diferentes, cuyos valores estan dados en curva1.dat, curva2.dat y curva3.dat.

  • Dibujar la figura junto con el polinomio interpolador.
  • Cambiar las condiciones de contorno a clamped boundary conditions, modificando el programa dado, y comparar los resultados.


    Darío Mitnik