Laboratorio 2 - Métodos Numéricos

Sergio Agudelo Bernal
Código: 285675


1. Programe y explique el algoritmo de interpolación de Lagrange y úselo para interpolar la siguiente tabla de valores, y compárelo con los resultados obtenidos del módulo interpolate de SciPy:

</tr> </table> </center>

Primero, se inicializan los arreglos de entrada y salida:

$$x$$

</td>

$$f(x)$$

1

7.93

2

10.05

3

12.66

4

15.79

5

19.47

6

23.73

7

28.60

8

34.11


In [1]:
%pylab inline
%load_ext sympy.interactive.ipythonprinting

from IPython.display import display, Math, Latex

import scipy.interpolate
from pylab import *
import sympy as sym

x1 = array([1, 2, 3, 4, 5, 6, 7, 8])
y1 = array([7.93, 10.05, 12.66, 15.79, 19.47, 23.73, 28.60, 34.11])

x = sym.Symbol('x')


Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

El polinomio de Lagrange es algo de este estilo:

$$P_n(x) = \sum\limits_{k = 0}^n y_k L_{n,k}$$ En donde:

$$L_{n,k} = \frac{\prod\limits_{j = 0}^n (x-x_j)}{\prod\limits_{j = 0}^n (x_k-x_j)}\ \ \ \ ,\ j \neq k$$ Debido a la condición $j \neq k$, el algoritmo debe excluir el producto $(x_j-x_k)$ en cada término de la sumatoria, cada vez que no se cumpla.

El siguiente algoritmo se compone de dos ciclos for, el segundo arma cada término $L_{n,k}$ y el primero suma paulatinamente cada uno de ellos, además de repetir el segundo para cada valor de $x$:


In [2]:
pol = x**0 - 1
for i in range(0, len(x1)):
    eq = x**0
    coeff = 1
    for j in range(0, len(x1)):
        if (j != i):
            eq = eq * (x - x1[j])
            coeff = coeff * (x1[i] - x1[j])
    pol = pol + ((y1[i] / coeff)*eq)
pol


Out[2]:
$$- 0.0015734126984127 \left(x - 8.0\right) \left(x - 7.0\right) \left(x - 6.0\right) \left(x - 5.0\right) \left(x - 4.0\right) \left(x - 3.0\right) \left(x - 2.0\right) + 0.0139583333333333 \left(x - 8.0\right) \left(x - 7.0\right) \left(x - 6.0\right) \left(x - 5.0\right) \left(x - 4.0\right) \left(x - 3.0\right) \left(x - 1.0\right) - 0.05275 \left(x - 8.0\right) \left(x - 7.0\right) \left(x - 6.0\right) \left(x - 5.0\right) \left(x - 4.0\right) \left(x - 2.0\right) \left(x - 1.0\right) + 0.109652777777778 \left(x - 8.0\right) \left(x - 7.0\right) \left(x - 6.0\right) \left(x - 5.0\right) \left(x - 3.0\right) \left(x - 2.0\right) \left(x - 1.0\right) - 0.135208333333333 \left(x - 8.0\right) \left(x - 7.0\right) \left(x - 6.0\right) \left(x - 4.0\right) \left(x - 3.0\right) \left(x - 2.0\right) \left(x - 1.0\right) + 0.098875 \left(x - 8.0\right) \left(x - 7.0\right) \left(x - 5.0\right) \left(x - 4.0\right) \left(x - 3.0\right) \left(x - 2.0\right) \left(x - 1.0\right) - 0.0397222222222222 \left(x - 8.0\right) \left(x - 6.0\right) \left(x - 5.0\right) \left(x - 4.0\right) \left(x - 3.0\right) \left(x - 2.0\right) \left(x - 1.0\right) + 0.00676785714285714 \left(x - 7.0\right) \left(x - 6.0\right) \left(x - 5.0\right) \left(x - 4.0\right) \left(x - 3.0\right) \left(x - 2.0\right) \left(x - 1.0\right)$$

La ecuación resultante es la siguiente:


In [3]:
sym.expand(pol)


Out[3]:
$$2.08166817117217 \times 10^{-17} x^{7} + 2.22044604925031 \times 10^{-16} x^{6} + 1.06581410364015 \times 10^{-14} x^{5} - 2.8421709430404 \times 10^{-14} x^{4} + 0.00499999999976808 x^{3} + 0.214999999999691 x^{2} + 1.4400000000004 x + 6.27000000000001$$

Ahora, se comparan los datos con la ecuación


In [4]:
C1 = array([2.081668e-17, 2.220446e-16, 1.065814104e-14, -2.84217e-14, 0.005, 0.215, 1.44, 6.27])

def f1(X):
    return C1[0]*X**7 + C1[1]*X**6 + C1[2]*X**5 + C1[3]*X**4 + C1[4]*X**3 + C1[5]*X**2 + C1[6]*X + C1[7]

xx1 = linspace(1, 8, num=1000)
p1, = plot(x1, y1, 'ko-')
p2, = plot(xx1, f1(xx1), 'r')

xlim([1,8])
ylim([7,35])
legend([p1, p2], ["Datos", "Interpolación"], loc=2)
show()


Según interpolate, la ecuación es:


In [5]:
Interp1 = scipy.interpolate.lagrange(x1, y1)
Interp1


Out[5]:
poly1d([  4.33680869e-18,  -4.16333634e-16,   0.00000000e+00,
        -2.84217094e-14,   5.00000000e-03,   2.15000000e-01,
         1.44000000e+00,   6.27000000e+00])

Ahora, se grafican todas las funciones:


In [6]:
def f2(X):
    return Interp1[7]*X**7 + Interp1[6]*X**6 + Interp1[5]*X**5 + Interp1[4]*X**4 + Interp1[3]*X**3 + Interp1[2]*X**2 + Interp1[1]*X + Interp1[0]

p1, = plot(x1, y1, 'ko-')
p2, = plot(xx1, f1(xx1), 'r')
p3, = plot(xx1, f2(xx1), 'g')

xlim([1,8])
ylim([7,35])
legend([p1, p2, p3], ["Datos", "Algoritmo", "Scipy"], loc=2)
show()


Los resultados son bastante similares, tanto numéricamente como gráficamente. Nótese que en este caso SciPy redondea algunos términos, pero como suele pasar en Python se entregan números con potencias muy bajas (~$\times 10^{-16}$) que deberían ser tomados como nulos.


2. Programe y explique el algoritmo para hallar el polinomio de interpolación de Newton y úselo para interpolar la tabla de valores siguiente, y compárelo con los resultados del módulo interpolate de SciPy:

</tr> </table> </center>

En este caso se deben armar varias ecuaciones, y se van sumando iterativamente. Antes de iniciar, se ingresa la tabla de datos:

$$x$$

</td>

$$f(x)$$

1.5

0

2.5

2.3

3.5

4.2

4.5

2.7

5.5

3.2

6.5

3.7

7.5

3.0

8.5

4.3

9.5

4.5

10.5

4.7

11.5

3.9

12.5

4.1


In [7]:
x2 = array([1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5])
y2 = array([0, 2.3, 4.2, 2.7, 3.2, 3.7, 3.0, 4.3, 4.5, 4.7, 3.9, 4.1])

El algoritmo consiste en hacer cálculos de diferencias divididas (e.g. pendiente de una recta), teniendo como parámetros de entrada los coeficientes anteriormente calculados (es decir, es un algoritmo recursivo).


In [8]:
Diff1 = zeros(len(x2)-1)
Diff2 = zeros(len(x2)-2)
Diff3 = zeros(len(x2)-3)
Diff4 = zeros(len(x2)-4)
Diff5 = zeros(len(x2)-5)
Diff6 = zeros(len(x2)-6)
Diff7 = zeros(len(x2)-7)
Diff8 = zeros(len(x2)-8)
Diff9 = zeros(len(x2)-9)
Diff10 = zeros(len(x2)-10)
Diff11 = zeros(len(x2)-11)

eq2 = x**0
eq3 = x**0
eq4 = x**0
eq5 = x**0
eq6 = x**0
eq7 = x**0
eq8 = x**0
eq9 = x**0
eq10 = x**0
eq11 = x**0
eq12 = x**0

#Primera Iteración
for i in range(1, len(x2)):
    Diff1[i-1] = ((y2[i] - y2[i-1])/(x2[i] - x2[i-1]))
    eq2 = eq2 * (x - x2[i-1])
    
#Segunda Iteración
for i in range(2,len(x2)):
    Diff2[i-2] = ((Diff1[i-1] - Diff1[i-2])/(x2[i] - x2[i-2]))
    eq3 = eq3 * (x - x2[i-2])

#Tercera Iteración
for i in range(3,len(x2)):
    Diff3[i-3] = ((Diff2[i-2] - Diff2[i-3])/(x2[i] - x2[i-3]))
    eq4 = eq4 * (x - x2[i-3])
    
#Cuarta Iteración
for i in range(4, len(x2)):
    Diff4[i-4] = ((Diff3[i-3] - Diff3[i-4])/(x2[i] - x2[i-4]))
    eq5 = eq5 * (x - x2[i-4])
    
#Quinta Iteración
for i in range(5,len(x2)):
    Diff5[i-5] = ((Diff4[i-4] - Diff4[i-5])/(x2[i] - x2[i-5]))
    eq6 = eq6 * (x - x2[i-5])

#Sexta Iteración
for i in range(6,len(x2)):
    Diff6[i-6] = ((Diff5[i-5] - Diff5[i-6])/(x2[i] - x2[i-6]))
    eq7 = eq7 * (x - x2[i-6])
    
#Séptima Iteración
for i in range(7, len(x2)):
    Diff7[i-7] = ((Diff6[i-6] - Diff6[i-7])/(x2[i] - x2[i-7]))
    eq8 = eq8 * (x - x2[i-7])
    
#Octava Iteración
for i in range(8,len(x2)):
    Diff8[i-8] = ((Diff7[i-7] - Diff7[i-8])/(x2[i] - x2[i-8]))
    eq9 = eq9 * (x - x2[i-8])

#Novena Iteración
for i in range(9,len(x2)):
    Diff9[i-9] = ((Diff8[i-8] - Diff8[i-9])/(x2[i] - x2[i-9]))
    eq10 = eq10 * (x - x2[i-9])

#Décima Iteración
for i in range(10, len(x2)):
    Diff10[i-10] = ((Diff9[i-9] - Diff9[i-10])/(x2[i] - x2[i-10]))
    eq11 = eq11 * (x - x2[i-10])
    
#Undécima Iteración
for i in range(11,len(x2)):
    Diff11[i-11] = ((Diff10[i-10] - Diff10[i-11])/(x2[i] - x2[i-11]))
    eq12 = eq12 * (x - x2[i-11])

print (Diff1)
print (Diff2)
print (Diff3)
print (Diff4)
print (Diff5)
print (Diff6)
print (Diff7)
print (Diff8)
print (Diff9)
print (Diff10)
print (Diff11)


[ 2.3  1.9 -1.5  0.5  0.5 -0.7  1.3  0.2  0.2 -0.8  0.2]
[-0.2  -1.7   1.    0.   -0.6   1.   -0.55  0.   -0.5   0.5 ]
[-0.5         0.9        -0.33333333 -0.2         0.53333333 -0.51666667
  0.18333333 -0.16666667  0.33333333]
[ 0.35       -0.30833333  0.03333333  0.18333333 -0.2625      0.175      -0.0875
  0.125     ]
[-0.13166667  0.06833333  0.03       -0.08916667  0.0875     -0.0525
  0.0425    ]
[ 0.03333333 -0.00638889 -0.01986111  0.02944444 -0.02333333  0.01583333]
[-0.0056746  -0.0019246   0.00704365 -0.00753968  0.00559524]
[ 0.00046875  0.00112103 -0.00182292  0.00164187]
[  7.24757496e-05  -3.27105379e-04   3.84975750e-04]
[ -3.99581129e-05   7.12081129e-05]
[  1.01060205e-05]

In [9]:
pol2 = y2[0] + Diff1[0]*eq12 + Diff2[0]*eq11 + Diff3[0]*eq10 + Diff4[0]*eq9 + Diff5[0]*eq8 + Diff6[0]*eq7 + Diff7[0]*eq6 + Diff8[0]*eq5 + Diff9[0]*eq4 + Diff10[0]*eq3 + Diff11[0]*eq2

La ecuación resultante es:


In [10]:
sym.expand(pol2)


Out[10]:
$$1.01060205226872 \times 10^{-5} x^{11} - 0.000762538580246913 x^{10} + 0.0253979965828924 x^{9} - 0.491621920469577 x^{8} + 6.12700357556217 x^{7} - 51.4562383535879 x^{6} + 296.103871743069 x^{5} - 1163.1349817648 x^{4} + 3044.65061793133 x^{3} - 5039.05929855976 x^{2} + 4730.43501988156 x - 1900.53331413269$$

Se define la función, tomando los coeficientes resultantes; y además se prueba con interpolate:


In [11]:
def f3(X):
    c1 = 0.0000101060205226872*X**11
    c2 = 0.000762538580246913*X**10
    c3 = 0.0253979965828924*X**9
    c4 = 0.491621920469577*X**8
    c5 = 6.12700357556217*X**7
    c6 = 51.4562383535879*X**6
    c7 = 296.103871743069*X**5
    c8 = 1163.1349817648*X**4
    c9 = 3044.65061793133*X**3
    c10 = 5039.05929855976*X**2
    c11 = 4730.43501988156*X
    c12 = 1900.53331413269
    return c1 - c2 + c3 - c4 + c5 - c6 + c7 - c8 + c9 - c10 + c11 - c12

Interp2 = scipy.interpolate.lagrange(x2, y2)
Interp2


Out[11]:
poly1d([  1.01060205e-05,  -7.62538580e-04,   2.53979966e-02,
        -4.91621920e-01,   6.12700358e+00,  -5.14562384e+01,
         2.96103872e+02,  -1.16313498e+03,   3.04465062e+03,
        -5.03905930e+03,   4.73043502e+03,  -1.90053331e+03])

Finalmente se comparan gráficamente las aproximaciones y la tabla de datos:


In [12]:
def f4(X):
    return Interp2[11]*X**11 + Interp2[10]*X**10 + Interp2[9]*X**9 + Interp2[8]*X**8 + Interp2[7]*X**7 + Interp2[6]*X**6 + Interp2[5]*X**5 + Interp2[4]*X**4 + Interp2[3]*X**3 + Interp2[2]*X**2 + Interp2[1]*X + Interp2[0]

xx2 = linspace(1.5, 12.5, num=1000)

p4, = plot(x2, y2, 'ko-')
p5, = plot(xx2, f3(xx2), 'r')
p6, = plot(xx2, f4(xx2), 'g')

xlim([1.5,12.5])
ylim([0,7])
legend([p4, p5, p6], ["Datos", "Algoritmo", "Scipy"], loc=2)
show()


De nuevo, la diferencia numérica y gráfica entre el algoritmo diseñado y la función interpolate es casi nula (sólo en algunos coeficientes se nota desde el sexto decimal).