Este es el código completo, incluyendo las figuras, del artículo http://pybonacci.org/2013/08/15/ajuste-e-interpolacion-unidimensionales-basicos-en-python-con-scipy/.
import numpy as np
from scipy.interpolate import barycentric_interpolate
def runge(x):
"""Función de Runge."""
return 1 / (1 + x ** 2)
N = 11 # Nodos de interpolación
xp = np.arange(11) - 5 # -5, -4, -3, ..., 3, 4, 5
fp = runge(xp)
x = np.linspace(-5, 5, num=1000)
y = barycentric_interpolate(xp, fp, x)
%matplotlib inline
import matplotlib.pyplot as plt
l, = plt.plot(x, y)
plt.plot(x, runge(x), '--', c=l.get_color())
plt.plot(xp, fp, 'o', c=l.get_color())
leg = plt.legend(['Interpolación', 'Real'])
leg.get_frame().set_facecolor('#fafafa')
from numpy.polynomial import chebyshev
coeffs_cheb = [0] * 11 + [1]
T11 = chebyshev.Chebyshev(coeffs_cheb, [-5, 5])
xp_ch = T11.roots()
fp_ch = runge(xp_ch)
y_ch = barycentric_interpolate(xp_ch, fp_ch, x)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
l1, = ax1.plot(x, y)
ax1.plot(x, runge(x), '--', c=l1.get_color())
ax1.plot(xp, fp, 'o', c=l1.get_color())
ax1.set_title("Nodos equiespaciados")
leg1 = ax1.legend(['Interpolación', 'Real'])
leg1.get_frame().set_facecolor('#fafafa')
l2, = ax2.plot(x, y_ch)
ax2.plot(x, runge(x), '--', c=l2.get_color())
ax2.plot(xp_ch, fp_ch, 'o', c=l2.get_color())
ax2.set_ylim(ax1.get_ylim())
ax2.set_title("Nodos de Chebyshev")
leg2 = ax2.legend(['Interpolación', 'Real'])
leg2.get_frame().set_facecolor('#fafafa')
Pato
from scipy.interpolate import InterpolatedUnivariateSpline
# Pato
P = [(0.9, 1.3), (1.3, 1.5), (1.9, 1.8), (2.1,2.1), (2.6, 2.6), (3.0, 2.7),
(3.9, 2.3), (4.4, 2.1), (4.8, 2.0), (5.0, 2.1), (6, 2.2), (7, 2.3),
(8, 2.2), (9.1, 1.9), (10.5, 1.4), (11.2, 0.9), (11.6, 0.8), (12, 0.6),
(12.6, 0.5), (13, 0.4), (13.2, 0.2)]
xi, yi = zip(*P)
print("Número de puntos: {}".format(len(xi)))
x = np.linspace(min(xi), max(xi), num=1001)
y1d = np.interp(x, xi, yi)
ysp = InterpolatedUnivariateSpline(xi, yi)(x)
plt.plot(xi, yi, 'x', mew=2)
plt.plot(x, y1d)
plt.plot(x, ysp)
leg = plt.legend(['Puntos', 'Lineal', 'Cúbico'])
leg.get_frame().set_facecolor('#fafafa')
Número de puntos: 21
import numpy.polynomial as P
# Cargamos los datos
data = np.loadtxt("polar.dat")
C_L, C_D = data
# Descarto los datos que no me sirven
stall_idx = np.argmax(C_L)
y = C_D[:stall_idx + 1]
x = C_L[:stall_idx + 1] ** 2
# Ajuste lineal, devuelve los coeficientes en orden creciente
C_D0, k = P.polynomial.polyfit(x, y, deg=1)
print(C_D0, k)
0.0221555792339 0.0380587917177
C_L_dom = np.linspace(C_L[0], C_L[stall_idx], num=1001)
C_D_int = P.polynomial.polyval(C_L_dom ** 2, (C_D0, k))
l, = plt.plot(C_D_int, C_L_dom)
plt.plot(C_D, C_L, 's', mew=2, mec=l.get_color(), mfc="none")
plt.title("Curva polar")
plt.xlabel("$C_D$")
plt.ylabel("$C_L$")
leg = plt.legend(["Ajuste parabólico", "Datos reales"], 4)
leg.get_frame().set_facecolor('#fafafa')
def func(x, A, B, C):
"""Modelo para nuestros datos."""
return A * np.exp(-B * x ** 2) + C
from scipy.optimize import curve_fit
Adat, Bdat, Cdat = 2.5, 1.3, 0.5
xdat = np.linspace(-2, 4, 12)
ydat = func(xdat, Adat, Bdat, Cdat) + 0.2 * np.random.normal(size=len(xdat))
(A, B, C), _ = curve_fit(func, xdat, ydat)
print(A, B, C)
2.41559988702 1.53213189998 0.497195046857
x = np.linspace(-2, 4, 1001)
l, = plt.plot(x, func(x, Adat, Bdat, Cdat), '--', label="Sin ruido")
plt.plot(xdat, ydat, 'o', mew=2, mfc="none", mec=l.get_color(), label="Datos")
plt.plot(x, func(x, A, B, C), label="Ajuste")
plt.ylim(0, 3.5)
leg = plt.legend()
leg.get_frame().set_facecolor('#fafafa')