Na podstawie Chebyshev interpolation
import matplotlib.pyplot as plt
from scipy import interpolate
from numpy import cos, pi, array, linspace
def cauchy(x):
return (1 + x**2)**-1
n = 25
x = linspace(-5, 5, n) # points to interpolate at
y = cauchy(x)
f = interpolate.BarycentricInterpolator(x, y)
Najpierw wykres interpolowanej funkcji
plt.plot(x,y)
[<matplotlib.lines.Line2D at 0x7fe578474100>]
xnew = linspace(-5, 5, 200) # points to plot at
ynew = f(xnew)
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.show()
Efekt interpolacjie nie jest zadawalający i łatwo można się przekonać, że im gęściej (wystarczy powyżej zmieniać wartość $n$), równomiernie będziemy umieszczali węzły interpolacji to wynik będzie tylko gorszy.
Remedium może być umieszczenie tych węzłów tam gdzie są miejsca zerowe wielomianu Czebyszewa odpowiedniego stopnia. Wówczas jakość interpolacji poprawia się w istotny sposób.
x = [cos(pi*(2*k-1)/(2*n)) for k in range(1, n+1)]
x = 5*array(x)
Polecenie x = 5*array(x)
„skaluje” te węzły tak, żeby były rozmieszczone na rozważanym przedziale $[-5, 5]$, a nie na przedziale $-1 \le x \le 1$ .
y = cauchy(x)
f = interpolate.BarycentricInterpolator(x, y)
xnew = linspace(-5, 5, 200) # points to plot at
ynew = f(xnew)
yy = cauchy(xnew)
plt.plot(x, y, 'o',label='wezły')
plt.plot(xnew, ynew, '-',label='interpolacja')
plt.plot(xnew, yy, '-.',label='funkcja')
plt.legend()
plt.show()
W przypadku nieparzystej liczby węzłów poprawnie będzie interpolowana wartość funkcji w maksimum.
Zobaczmy jak wygląda błąd interpolacji.
plt.plot(xnew,ynew-yy,label='bład interpolacji')
plt.legend()
plt.show()