# Interpolación polinomial
## Vandermonde, Lagrange, Diferencias divididas de Newton
### SCT 2020


```python
import numpy as np
import matplotlib.pyplot as plt
```

### Implementación de métodos de interpolación


```python
def vandermonde(x_i, y_i):
    n = x_i.shape[-1]
    A = np.vander(x_i, increasing=True)
    a = np.linalg.solve(A, y_i)
    V = lambda x: np.dot(np.array([x ** i for i in range(n)]), a)
    return np.vectorize(V)
```


```python
def lagrange(x_i, y_i):
    n = x_i.shape[-1]
    L = lambda x: np.dot(y_i, np.array([np.prod(x - np.delete(x_i, k)) 
                                        / np.prod(x_i[k] - np.delete(x_i, k)) for k in range(n)]))
    return np.vectorize(L)
```


```python
def barycentric(x_i, y_i):
    n = x_i.shape[-1]
    w = 1 / np.array([np.prod(x_i[i] - np.delete(x_i, i)) for i in range(n)]) 
    b1 = lambda x: y_i[np.where(np.in1d(x_i, x))] 
    b2 = lambda x: np.dot(y_i, w /(x - x_i))/ np.dot(w, 1 / (x - x_i))
    B = lambda x: b1(x) if x in x_i else b2(x)
    return np.vectorize(B)
```


```python
def newtonDD(x_i, y_i):
    n = x_i.shape[-1]
    pyramid = np.zeros((n, n)) # Create a square matrix to hold pyramid
    pyramid[:,0] = y_i # first column is y
    for j in range(1,n):
        for i in range(n-j):
            # create pyramid by updating other columns
            pyramid[i][j] = (pyramid[i+1][j-1] - pyramid[i][j-1]) / (x_i[i+j] - x_i[i])
    a = pyramid[0] # f[ ... ] coefficients
    N = lambda x: a[0] + np.dot(a[1:], np.array([np.prod(x - x_i[:i]) for i in range(1, n)]))
    return np.vectorize(N)
```


```python
f1 = lambda x: 1 / (1 + 25 * x ** 2) # Runge Function
f2 = lambda x: np.piecewise(x, [np.abs(x) <= 1e-14], [1])
```

### Gráfico de funciones


```python
xx = np.linspace(-2, 2, 1001)
plt.figure(figsize=(12, 6))
plt.plot(xx, f1(xx), label=r"$f_1(x)$")
plt.plot(xx, f2(xx), label=r"$f_2(x)$")
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.legend()
plt.grid(True)
plt.show()
```

Interpolación utilizando puntos de $f_1(x)$.


```python
f = f1
```


```python
N_i = 13
x_a, x_b = -1, 1
x_i = np.linspace(x_a, x_b, N_i)
y_i = f(x_i)
```


```python
Pv = vandermonde(x_i, y_i)
```


```python
Pl = lagrange(x_i, y_i)
```


```python
Pb = barycentric(x_i, y_i)
```


```python
Pn = newtonDD(x_i, y_i)
```


```python
N_e = 500
x_e = np.linspace(x_a, x_b, N_e)
y_e = f(x_e)
```

### Gráfico de función y el resultado de interpolación


```python
plt.figure(figsize=(14, 8))
plt.plot(x_e, y_e, 'b-', label=r"$f(x)$")
plt.plot(x_i, y_i, 'ro', label="Puntos para interpolar")
plt.plot(x_e, Pv(x_e), 'r--', label=r'$P_V(x)$')
plt.plot(x_e, Pl(x_e), 'b--', label=r'$P_L(x)$')
plt.plot(x_e, Pb(x_e), 'g--', label=r'$P_B(x)$')
plt.plot(x_e, Pn(x_e), 'k--', label=r'$P_N(x)$')
plt.grid(True)
plt.legend()
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()
```

Notamos como el interpolador efectivamente pasa por los puntos de interpolación, pero aparece el **Fenómeno de Runge**.

# Nodos de Chebyshev

Recordar que los puntos de Chebyshev se definen como:
\begin{equation}
    x_i = \cos\left(\frac{(2i - 1)\pi}{2n}\right), \quad i=1, \dots, n
\end{equation}


```python
def chebyshevNodes(n):
    i = np.arange(1, n+1)
    t = (2*i - 1) * np.pi / (2 * n)
    return np.cos(t)
```


```python
N_c = 13
```


```python
xc_i = chebyshevNodes(N_c)
```


```python
yc_i = f(xc_i)
```

### Puntos de Chebyshev

Visualización de algunos puntos...


```python
t = np.linspace(0, np.pi, 100)
plt.figure(figsize=(12, 6))
plt.plot(np.cos(t), np.sin(t))
plt.plot(xc_i, np.zeros(N_c), 'bo')
plt.plot(xc_i, xc_i * np.tan((2 * np.arange(1, N_c+1) - 1) * np.pi / (2 * N_c)), 'ro')
plt.axhline(y=0, color='k')
plt.axis('scaled')
plt.grid(True)
plt.xlabel(r"$x$")
plt.show()
```

### Obtención de polinomios utilizando los puntos de Chebyshev


```python
Pv_c = vandermonde(xc_i, yc_i)
```


```python
Pl_c = lagrange(xc_i, yc_i)
```


```python
Pb_c = barycentric(xc_i, yc_i)
```


```python
Pn_c = newtonDD(xc_i, yc_i)
```


```python
xc_e = chebyshevNodes(500)
```

### Gráfico utilizando los puntos de Chebysev


```python
plt.figure(figsize=(12, 6))
plt.plot(xc_e, f(xc_e), 'b-', label=r"$f(x)$")
plt.plot(xc_i, yc_i, 'ro', label="Puntos para interpolar")
plt.plot(xc_e, Pv_c(xc_e), 'r--', label=r'$P_V(x)$')
plt.plot(xc_e, Pl_c(xc_e), 'b--', label=r'$P_L(x)$')
plt.plot(xc_e, Pb_c(xc_e), 'g--', label=r'$P_B(x)$')
plt.plot(xc_e, Pn_c(xc_e), 'k--', label=r'$P_N(x)$')
plt.grid(True)
plt.legend()
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()
```

### Comparación interpoladores

Comparamos un error relativo entre el interpolador utilizando puntos equiespaciados y los de Chebyshev.


```python
error = lambda y, p: np.linalg.norm(y - p, np.inf) / np.linalg.norm(y, np.inf)
```


```python
print("Error puntos equiespaciados:", error(f(x_e), Pb(x_e))) # Equispaced points interpolation comparison
print("Error puntos Chebyshev:", error(f(x_e), Pb_c(x_e))) # Chebyshev points interpolation comparison
```

    Error puntos equiespaciados: 3.663235759944772
    Error puntos Chebyshev: 0.06919634037617944


# Conclusión

En este ejemplo/ejercicio se puede ver como afecta la calidad de la interpolación con el uso de puntos de Chebyshev en comparación a puntos equiespaciados.

# Referencias
* https://medium.com/@sddkal/newtons-divided-difference-method-for-polynomial-interpolation-4bc094ba90d7

7016ff58c10e49c4ca814f46c70dfa934556d743
ipynb
Jupyter Notebook
material/05_interpolacion_1D/interpolacion.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
2020-04-24T01:25:24.000Z
2021-05-10T01:08:37.000Z
material/05_interpolacion_1D/interpolacion.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
2020-04-22T00:57:29.000Z
2020-04-23T23:59:28.000Z
material/05_interpolacion_1D/interpolacion.ipynb
fitoxgs2/INF-285
f7d547c2480aea7eb1242b873ab1135a0873834f
BSD-3-Clause
2020-04-20T01:15:09.000Z
2020-07-27T22:53:59.000Z
Qwen/Qwen-72B
1. YES
2. YES