Coś jeszcze o interpolacji

Na podstawie Chebyshev interpolation

In [1]:
import matplotlib.pyplot as plt
from scipy import interpolate
from numpy import cos, pi, array, linspace
In [2]:
def cauchy(x):
    return (1 + x**2)**-1
In [3]:
n = 25
x = linspace(-5, 5, n) # points to interpolate at
In [4]:
y = cauchy(x)
f = interpolate.BarycentricInterpolator(x, y)

Najpierw wykres interpolowanej funkcji

In [5]:
plt.plot(x,y)
Out[5]:
[<matplotlib.lines.Line2D at 0x7fe578474100>]
In [6]:
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.

In [7]:
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$ .

In [8]:
y = cauchy(x)
f = interpolate.BarycentricInterpolator(x, y)
In [9]:
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.

In [10]:
plt.plot(xnew,ynew-yy,label='bład interpolacji')
plt.legend()
plt.show()
In [ ]: