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:
$$x$$ | $$f(x)$$ | </tr>
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')
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]:
In [3]:
sym.expand(pol)
Out[3]:
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()
In [5]:
Interp1 = scipy.interpolate.lagrange(x1, y1)
Interp1
Out[5]:
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 sí 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:
$$x$$ | $$f(x)$$ | </tr>
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])
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)
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
In [10]:
sym.expand(pol2)
Out[10]:
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]:
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).