Aproksymacja wielomianowa

Biblioteka numpy zawiera procedury do aproksymacji. W szczególności:

In [1]:
import matplotlib.pyplot as plt
import numpy as np
#from scipy.interpolate import interp1d
#from scipy.interpolate import lagrange
#from scipy.interpolate import BarycentricInterpolator

Plik AVERAGE300aprox.dat zawiera dane już oczyszczone.

In [2]:
data = np.loadtxt("AVERAGE300aprox.dat")
In [3]:
plt.plot(data[:,0],data[:,1],'.')
plt.show()

Dokonamy najprostszej aproksymacji wielomianem stopnia trzeciego, używając funkcji polyfit().

Funkcja zwraca współczynniki wielomianu (w przypadku gdy stopień, jak w przykładzie jest równy 3) $$z_0 x^3 + z_1 x^2 + z_2 x + z_3$$

In [4]:
z = np.polyfit(data[:,0],data[:,1], 3)

/home/myszka/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3418: RankWarning: Polyfit may be poorly conditioned
  exec(code_obj, self.user_global_ns, self.user_ns)

Jak widać funkcja informuje, że zadanie jest źle uwarunkowane. Wynika to (w tym przypadku) z tego, że wartości $x$ są bardzo duże (gdyż reprezentują czas bezwzględny w sekundach). Bardzo łatwo zamienić je na wartości reprezentujące czas względny (względem początku doby) w sekundach odejmując od każdej wartości odciętej minimum.

In [5]:
tmin = min(data[:,0])
data[:,0] = data[:,0] - tmin

z = np.polyfit(data[:,0],data[:,1], 3)

Jak widać jest wyraźnie lepiej (nie pojawia się komunikat).

z jest tablicą zwierającą współczynniki wielomianu. Zwracam uwagę, że ze względu na spore wartości $x$ (liczba sekund w ciągu doby) współczynniki są bardzo małe. Grozi to utratą dokładności podczas obliczeń. Można spróbować zamienić sekundy na godziny dzieląc odcięte przez 3600.

Jeszcze lepsze efekty osiągniemy normalizując odcięte danych do zakresu $[-1, 1]$. Daje to szansę, że skrajne wartości podnoszone do dużych nawet potęg będą przyjmowały sensowniejsze wartości.

In [6]:
z
Out[6]:
array([-1.50077850e-15, -4.89454863e-09,  3.67306502e-04, -2.29165068e+00])
In [7]:
data[:,0] = data[:,0]/3600.

z = np.polyfit(data[:,0],data[:,1], 3)
z
Out[7]:
array([-7.00203219e-05, -6.34333503e-02,  1.32230341e+00, -2.29165068e+00])

Współczynniki wielomianu uległy modyfikacji, za wyjątkiem ostatniego $z_3$ reprezentującego wartość funkcji w punkcie 0.

Do manipulacji na wielomianie warto używac funkcji poly1d()

In [8]:
wielomian = np.poly1d(z)
In [9]:
wielomian(0.)
Out[9]:
-2.2916506803099526
In [10]:
wielomian(24)
Out[10]:
-8.06193960023452

Obejrzyjmy na wykresie rezultat aproksymacji.

In [11]:
x=np.linspace(0, 24, 25)

plt.plot(x, wielomian(x),'g-',data[:,0], data[:,1], 'r.')
plt.show()

Przybliżenie nie jest idealne, ale (w pewnym sense) oddaje charakter dobowych zmian temperatury: w ciągu doby rosną, a na zakończenie dnia temperatura jest mniejsza niż na jego początku.

Niech N będzie to stopień wielomianu interpolującego. W poniższym przykłądzie można dowlonie(?) zmieniać N i przeliczać…

In [12]:
N = 18
plt.plot(x, np.poly1d(np.polyfit(data[:,0],data[:,1], N))(x),'g-', 
         data[:,0], data[:,1], 'r.')
plt.show()

Zobaczmy jak wyglądają różnice wartości prawdziwej i aproksymowanej.

In [13]:
err = data[:,1] - np.poly1d(np.polyfit(data[:,0],data[:,1], N))(data[:,0])
plt.plot(data[:,0], err)
plt.show()

Generalnie nie jest źle: wartości skupiają się wokół zera w sposób dosyć symetryczny. Zobaczmy histogram.

In [14]:
plt.hist(err)
plt.show()

Nie do końca przypomina on histogram rozkłądu normalnego, ale nie jest ,,najgorszy''.

Można też zdefiniować sobie funkcję (na przykład wykres) i korzystać z niej w celu zobaczenia rezultatu aproksymacji.

In [15]:
def wykres(N):
    plt.plot(x, np.poly1d(np.polyfit(data[:,0],data[:,1], N))(x),'g-', 
         data[:,0], data[:,1], 'r.')
    plt.show()

wykres(5)

Można też zabawić się w wykres interakcyjny, w którym suwaczkiem zmieniamy wartość parametru. Ale to już jakaś magia…

In [16]:
%matplotlib inline
#from __future__ import print_function
from ipywidgets import interactive
import matplotlib.pyplot as plt
import numpy as np

def wykres_wykres(N):
    plt.figure(2)
    plt.plot(x, np.poly1d(np.polyfit(data[:,0],data[:,1], N))(x),'g-', 
         data[:,0], data[:,1], 'r.')
    # plt.ylim(-5, 5)
    plt.show()

interactive_plot = interactive(wykres_wykres, N=(0, 30))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

Dla $N=0$ dostajemy średnią temperaturę dobową, a dla $N=1$ trend — temperatura spada.

Również dla kolejnych wartości $N$ można znaleźć jakąś interpretację…

Niestety, dla $N >18$ funkcja informuje, że zadanie staje się źle uwarunkowane. Jest to znowu związane z tym, że wartości $x$ podniesiona do potęgi 19 czy 20 staje się baaaardzo duża:

In [17]:
24.**20
Out[17]:
4.0199887178406037e+27

Spójrzmy zatem na to co składa się na wartość wielomainu dla czasu $x$.

In [18]:
N = 19
x = 24
a=np.polyfit(data[:,0],data[:,1], N)

/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:3: RankWarning: Polyfit may be poorly conditioned
  This is separate from the ipykernel package so we can avoid doing imports until
In [19]:
for i in range(0, N):
    print(a[i]*x**(N-i))

-11330548833.759817
91343839406.0988
-323839414937.3468
649910659240.3022
-768190706710.5015
429110280561.0846
186788204073.06516
-623109225162.7025
652222465812.5778
-426061915498.8431
193788821235.2673
-63305881981.175026
14896926263.546213
-2489387080.6460786
286464930.539979
-21521206.64991253
960976.6069019116
-21266.746996815928
173.36912084905677

Jak widać sumowane są liczby o bardzo dużej rozpiętości wartości. Musi prowadzić to do problemów…

In [ ]: