Aproksimacija#
Uvod#
V strojniški praksi se pogosto srečamo s tabelo podatkov, ki so lahko obremenjeni z merilnimi ali numeričnimi napakami.
Oglejmo si primer meritve (linearne) vzmeti (\(x\) je raztezek, \(y\) je pomerjena sila):
import numpy as np
x = np.array([0.1, 1.1, 2.05, 3.2, 3.9 ])
y = np.array([0.22, 18.15, 44.33, 75.59, 105.63])
Poglejmo si podatke na sliki, najprej uvozimo potrebne pakete:
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(x, y, 'o');
plt.xlabel('Pomik $x$')
plt.ylabel('Sila $y$');
Za konkreten primer bi bilo, glede na poznavanje fizikalnega ozadja linearne vzmeti, primerno, da bi meritve poskušali popisati z linearno funkcijo:
Poznamo tabelo \(n\) podatkov \(x_i, y_i\) za \(i=0,1,\dots,n-1\); teh je več, kot jih potrebujemo za določitev dveh konstant \(a_0\) in \(a_1\), zato imamo torej predoločen sistem linearnih enačb:
Iščemo taki vrednosti konstanti \(a_0\) in \(a_1\), da se bo funkcija \(f(x)\) v znanih točkah \(x_i\) najbolje ujemala z \(y_i\).
Najprej torej potrebujemo kriterij za najboljše ujemanje.
Za vrednosti iz tabele \(x_i, y_i\) bi lahko iskali vrednosti \(a_0\) in \(a_1\), pri katerih bi bila vsota absolutne vrednosti odstopkov \(S\) najmanjša:
Ker pa taka funkcija \(P(a_0, a_1)\) ni zvezno odvedljiva,raje uporabimo metodo najmanjših kvadratov:
Takšna funkcija \(S(a_0, a_1)\) je zvezna in zvezno odvedljiva. S parcialnim odvajanjem po parametrih \(a_0\) in \(a_1\) lahko najdemo stacionarno točko (parcialna odvoda sta enaka 0). Postopek si bomo za linearno funkcijo pogledali v naslednjem poglavju.
Metoda najmanjših kvadratov za linearno funkcijo#
from IPython.display import YouTubeVideo
YouTubeVideo('8-oixY6GG9E', width=800, height=300)
Poiskati moramo konstanti \(a_0\), \(a_1\), da bo vsota kvadratov razlik med funkcijo in tabelirano vrednostjo (\(x_i, y_i\), kjer \(i=0,1,\dots, n-1\) in je \(n\) število tabeliranih podatkov):
najmanjša. Vrednost bo najmanjša v stacionarni točki, ki jo določimo s parcialnim odvodom po parametrih \(a_0\) in \(a_1\).
Najprej izvedemo parcialni odvod po parametru \(a_0\):
Izraz uredimo:
Podobno postopamo še za \(a_1\):
Ker v stacionarni točki velja \(\partial S(a_0, a_1)/\partial a_0=0\) in \(\partial S(a_0, a_1)/\partial a_1=0\), iz zgornjih izrazov izpeljemo:
in
Dobili smo sistem dveh linearnih enačb za neznanki \(a_0\) in \(a_1\), ki ga znamo rešiti. Imenujemo ga normalni sistem (število enačb je enako številu neznank).
Zapišimo normalni sistem v matrični obliki:
A = [[np.sum(x**2), np.sum(x)], # matrika koeficientov
[np.sum(x), len(x)]]
b = [np.dot(y,x), np.sum(y)] # vektor konstant
A = np.asarray(A)
b = np.asarray(b)
print('A:', A)
print('b:', b)
A: [[30.8725 10.35 ]
[10.35 5. ]]
b: [764.7085 243.92 ]
Sedaj moramo rešiti linearni sistem:
Opomba, tukaj smo vektor neznank zapisali kot \(\mathbf{a}=(a_0, a_1)\).
Sistem rešimo:
a0, a1 = np.linalg.solve(A, b)
a0, a1
(np.float64(27.497258679085512), np.float64(-8.135325465707014))
Preverimo še število pogojenosti:
np.linalg.cond(A)
np.float64(25.20071350870242)
Sedaj si bomo pogledali še rezultat. Najprej pripravimo sliko, ki bo vsebovala tudi informacijo o vsoti kvadratov odstopanja \(f(x_i)\) od tabeliranih vrednosti \(y_i\).
def slika(naklon=a0, premik=a1):
d=1
def linearna_f(x, a0, a1):
return a0*x+a1
def S(x, y, f):
return np.sum(np.power(y-f,2))
plt.plot(x,y,'.', label='Tabela podatkov')
linearna_f1 = linearna_f(x, naklon, premik)
linearna_f1_MNK = linearna_f(x, a0, a1)
plt.plot(x, linearna_f1, '-', label='Izbrani parametri')
plt.plot(x, linearna_f1_MNK, '-', label='Metoda najmanjših kvadratov')
napaka = S(x, y, linearna_f(x, naklon, premik))
sprememba_napake_v_smeri_a0 = (S(x, y, linearna_f(x, naklon+d, premik))-napaka)/d
sprememba_napake_v_smeri_a1 = (S(x, y, linearna_f(x, naklon, premik+d))-napaka)/d
title = f'S: {napaka:g}, \
$\\Delta S/\\Delta a_0$: {sprememba_napake_v_smeri_a0:g}, \
$\\Delta S/\\Delta a_1$: {sprememba_napake_v_smeri_a1:g}'
plt.title(title)
plt.legend()
plt.ylim(-10,110)
plt.show()
from ipywidgets import interact
interact(slika, naklon=(0, 50, 2), premik=(-10, 10, 1));
Uporaba psevdo inverzne matrike#
Do podobnega rezultata lahko pridemo z uporabo psevdo inverzne matrike. Iščemo \(y(x)=a_0\,x+a_1\) in nastavimo predoločen sistem \(\mathbf{A}\,\mathbf{a}=\mathbf{y}\), kjer je matrika koeficientov \(\mathbf{A}\) definirana glede na vrednosti \(x_i\) (\(i=0,1,\dots,n-1\)):
A=np.array([x, np.ones_like(x)]).T
A
array([[0.1 , 1. ],
[1.1 , 1. ],
[2.05, 1. ],
[3.2 , 1. ],
[3.9 , 1. ]])
Vektor konstant smo označili z \(\mathbf{b}\), v našem primeru pa je to kar vektor vrednosti \(\mathbf{y}\) z elementi \(y_i\) (\(i=0,1,2,\dots,n-1\)):
y
array([ 0.22, 18.15, 44.33, 75.59, 105.63])
Preverimo rang matrike koeficientov:
np.linalg.matrix_rank(A)
np.int64(2)
In še rang razširjene matrike:
Ab = np.hstack((A,np.array([y]).T))
Ab
array([[1.0000e-01, 1.0000e+00, 2.2000e-01],
[1.1000e+00, 1.0000e+00, 1.8150e+01],
[2.0500e+00, 1.0000e+00, 4.4330e+01],
[3.2000e+00, 1.0000e+00, 7.5590e+01],
[3.9000e+00, 1.0000e+00, 1.0563e+02]])
np.linalg.matrix_rank(Ab)
np.int64(3)
Ker rešujemo sistem \(m\) linearnih enačbami z \(n\) neznankami ter velja \(m>n\) in je rang razširjene matrike \(n+1\), imamo predoločeni sistem.
Vektor konstant \(\mathbf{a}\) določimo z uporabo psevdo inverzne matrike:
np.linalg.pinv(A).dot(y)
array([27.49725868, -8.13532547])
Metoda najmanjših kvadratov za poljubni polinom#
Linearno aproksimacijo, predstavljeno zgoraj, bomo posplošili za poljubni polinom stopnje \(m\):
kjer \(f_s(x)=x^{m-s}\) imenujemo bazna funkcija (\(s=0,1,2,\dots,m\)).
Tabela podatkov naj bo definirana z \(x_i, y_i\), kjer je \(i=0,1,2,\dots,n-1\).
Opomba: zaradi kompaktnosti zapisa bomo konstante \(a\) zapisali v vektorski obliki \(\mathbf{a}=[a_0, a_1,\dots,a_m]\).
Uporabimo metodo najmanjših kvadratov:
Potreben pogoj za nastop ekstrema funkcije \(m+1\) neodvisnih spremenljivk je, da najdemo stacionarno točko za vsak \(a_v\), iščemo torej \(\partial S(\mathbf{a})/\partial a_v=0\) (namesto \(s\) smo uporabili indeks \(v\)).
Najprej določimo parcialni odvod za izbrani \(a_v\):
Opomba: \(\frac{\partial}{\partial a_v}\left(\sum_{s=0}^{m}a_s\,x_i^{m-s}\right)=x_i^{m-v}\).
Ker je parcialni odvod v stacionarni točki enak 0, zgornji izraz preoblikujemo:
Izraz uredimo:
Zamenjamo vrstni red seštevanja ter izpeljemo:
Izpeljali smo enačbo \(v\) sistema \(m+1\) linearnih enačb:
Element \(A_{v,s}\) matrike koeficientov je:
Element vektorja konstant je:
Numerični zgled#
Uporabimo podatke iz prve naloge in poskusimo aproksimirati s polinomom 2. stopnje (\(m=2\)).
Tabela podatkov je:
x
array([0.1 , 1.1 , 2.05, 3.2 , 3.9 ])
y
array([ 0.22, 18.15, 44.33, 75.59, 105.63])
Izračunajmo matriko koeficientov:
m = 2 #stopnja
A = np.zeros((m+1,m+1))
for v in range(m+1):
for s in range(m+1):
A[v,s] = np.sum(x**(2*m-v-s))
A
array([[355.32690625, 102.034125 , 30.8725 ],
[102.034125 , 30.8725 , 10.35 ],
[ 30.8725 , 10.35 , 5. ]])
Izračunajmo še vektor konstant:
b = np.zeros(m+1)
for v in range(m+1):
b[v] = np.dot(y,x**(m-v))
b
array([2588.934425, 764.7085 , 243.92 ])
Preverimo število pogojenosti:
np.linalg.cond(A)
np.float64(963.2125856290634)
Rešimo sistem:
a = np.linalg.solve(A, b)
a
array([ 3.17760897, 14.67380042, -1.21091344])
Glede na definicijo aproksimacijskega polinoma:
Kar v konkretnem primeru je aproksimacijski polinom:
Definirajmo numerično implementacijo:
def apr_polinom(x, a):
"""Vrne vrednosti aproksimacijskega polinoma
:param x: vrednosti kjer računamo aproksimirani rezultat
:param a: koeficienti aproksimacijskega polinoma
"""
m = len(a) - 1
return np.sum(np.asarray([_*x**(m-v) for v,_ in enumerate(a)]), axis=0)
Prikažemo:
x_g = np.linspace(np.min(x), np.max(x), 100) # več točk za prikaz
plt.plot(x, y, 'o', label='Tabela podatkov')
plt.plot(x_g, apr_polinom(x_g, a), lw=5, alpha=0.5, label='Aproksimacija s kvadratno funkcijo')
plt.legend();
Poglejmo še napako aproksimacije: $\(e_i=y_i - f(x_i)\)$
za \(i=0,1,2,\dots,n-1\).
Pri pravilno izvedni aproksimaciji je nekaj \(e_i\) pozitivnih in nekaj negativnih. Poglejmo, če je to res v našem primeru:
e = y - apr_polinom(x, a)
e
array([-0.06824269, -0.62517387, 2.1057209 , -2.69396374, 1.28165939])
Opomba: višje stopnje polinoma kot uporabimo, večja je verjetnost slabe pogojenosti. Iz tega razloga s stopnjo polinoma ne pretiravamo (v praksi uporabljamo predvsem nizke stopnje)!
Uporaba numpy za aproksimacjo s polinomom#
Poglejmo si, kako uporabimo knjižnico numpy za polinomsko aproksimacijo.
Najprej uporabimo funkcijo numpy.polyfit (dokumentacija):
polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)
ki zahteva tri parametre: x in y predstavljata tabelo podatkov (lahko tudi v obliki seznamov vektorjev), deg pa stopnjo polinoma. Ostali parametri so opcijski (npr. w za uporabo uteži pri aproksimaciji).
Funkcija polyfit vrne seznam koeficientov polinoma (najprej za najvišji red); rezultat je lahko tudi seznam seznamov (če so vhodni podatki seznam vektorjev).
Poglejmo si uporabo za predhodno obravnavani primer:
koef = np.polyfit(x, y, deg=2)
koef
array([ 3.17760897, 14.67380042, -1.21091344])
a # rezultat lastne implementacije
array([ 3.17760897, 14.67380042, -1.21091344])
Ko imamo koeficiente, lahko ustvarimo objekt polinoma s klicem numpy.poly1d (dokumentacija):
poly1d(c_or_r, r=False, variable=None)
kjer c_or_r predstavlja seznam koeficientov polinoma oz. ničle polinoma v primeru, da je r=True. Funkcija vrne instanco objekta, s klicem katere lahko izračunamo vrednosti aproksimacijskega polinoma pri x, lahko pa izračunamo tudi druge stvari, kot na primer ničle polinoma.
Poglejmo si primer:
p = np.poly1d(koef) # p je instanca objekta poly1d. Aproksimacijo sedaj dobimo z y = p(x).
Izrišimo vrednosti:
plt.plot(x, y, 'o', label='Tabela podatkov')
plt.plot(x_g, p(x_g), lw=5, alpha=0.5, label='Aproksimacija s kvadratno funkcijo')
plt.legend();
Izračunajmo ničle polinoma:
p.roots
array([-4.69897273, 0.08109792])
Aproksimacija s poljubno funkcijo#
Pri aproksimaciji nismo omejeni zgolj na polinome. Tabele podatkov lahko aproksimiramo:
z linearno kombinacijo linearno neodvisnih baznih funkcij ali
s funkcijo, v kateri nastopajo parametri v nelinearni zvezi (npr. \(a_0\,\sin(a_1\,x+a_2)\)).
Za podrobnosti glejte vir J. Petrišič: Uvod v Matlab za inženirje, Fakulteta za strojništvo 2013, str 145.
Osredotočili se bomo na uporabo scipy paketa za aproksimacijo z nelinearno fukcijo, ki temelji na metodi najmanjših kvadratov.
Aproksimacija s harmonsko funkcijo#
Tabela podatkov je definirana kot:
x = np.array([ 0.1, 0.8, 1.7, 2.5, 3.4, 4.2, 5.1])
y = np.array([ 0.01, -1.13, 0.02, 0.92, -0.01, -0.98, 0.1])
Prikažimo tabelo podatkov:
plt.plot(x, y, 'o', label='Tabela podatkov')
plt.legend();
Aproksimacijo z nelinearno funkcijo bomo izvedli s pomočjo scipy.optimize.curve_fit (dokumentacija):
curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(-inf, inf), method=None, jac=None, **kwargs)
katera zahteva tri parametre: f predstavlja definicijo Python funkcije, s katero želimo aproksimirati, in katere parametre spreminjamo z uporabo metode najmanjših kvadratov. xdata in ydata predstavljata tabelo podatkov. Priporočeno je tudi, da definiramo približek iskanih parametrov p0. Ostali parametri so opcijski.
Funkcija vrne dve numerični polji: popt, ki predstavlja najdene parametre ter pcov, ki predstavljajo ocenjeno kovarianco popt.
Definirajmo najprej Python funkcijo, katere prvi parameter je neodvisna spremenljivka x, nato pa sledijo parametri, ki jih želimo določiti:
def func(x, A, ω, ϕ):
return A*np.sin(ω*x+ϕ)
kjer je A amplitua, ω krožna frekvenca in ϕ faza harmonske funkcije. S pomočjo slike lahko ugibamo prve približke: A=1, ω=1, ϕ=0
Sedaj uvozimo curve_fit in izvedemo optimizacijski postopek:
from scipy.optimize import curve_fit
popt, pcov = curve_fit(func, x, y, p0=[1, 1, 0])
popt
array([ 1.01461946, 1.89652046, -3.27943097])
Izračunali smo pričakovane vrednosti (glejte zgoraj).
func(x, *popt)
array([-0.0525477 , -0.99608783, -0.05612665, 1.00860622, -0.02753942,
-1.01426499, 0.1110182 ])
x_g = np.linspace(np.min(x), np.max(x), 50)
y_g = func(x_g, *popt) # bodi pozorni kako smo v funkcijo posredovali parametre
plt.plot(x, y, 'o', label='Tabela podatkov')
plt.plot(x_g, y_g, label='Aproksimacija')
plt.legend();
Dodatno#
Naredite .exe svojega programa:
Poglejte pandas paket.
Aproksimacija z zlepki in uporabo SciPy#
Tabela podatkov naj bo:
x = np.linspace(-3, 3, 20)
x
array([-3. , -2.68421053, -2.36842105, -2.05263158, -1.73684211,
-1.42105263, -1.10526316, -0.78947368, -0.47368421, -0.15789474,
0.15789474, 0.47368421, 0.78947368, 1.10526316, 1.42105263,
1.73684211, 2.05263158, 2.36842105, 2.68421053, 3. ])
np.random.seed(0) # seme generatorja naključnih števil
y = np.exp(-x**2) + 0.1 * np.random.normal(scale=.5, size=len(x))
y
array([ 0.08832603, 0.02075073, 0.0526001 , 0.12684217, 0.14234432,
0.08387244, 0.34226064, 0.52862107, 0.79385312, 0.99590738,
0.98257964, 0.87172774, 0.57424082, 0.30083997, 0.15492949,
0.06565014, 0.08950146, -0.00659471, 0.01639626, -0.04258138])
Poglejmo si objekt scipy.interpolate.UnivariateSpline (dokumentacija), ki omogoča tako interpolacijo kot aproksimacijo z zlepki:
UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None, ext=0, check_finite=False)
Parametra x in y predstavljata tabelo podatkov.
Opcijski parameter s določa vrednost, katere vsota kvadratov razlik aproksimacijskega zlepka in aproksimacijskih točk ne sme preseči:
sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s
w so uteži posameznih točk.
Če definiramo s=0, zahtevamo interpolacijo.
Parameter k definira stopnjo polinomskega zlepka (privzeto je k=3).
Aproksimacijo z zlepki izvedemo tako, da ob tabeli podatkov x in y definiramo še parameter s. Izvedimo interpolacijo:
from scipy.interpolate import UnivariateSpline
spl = UnivariateSpline(x, y, s=0.)
x_g = np.linspace(-3, 3, 100)
plt.plot(x, y, 'o', ms=5, label='Aproksimacijske točke')
plt.plot(x_g, spl(x_g), lw=3, label='Zlepek', alpha=0.6)
plt.legend();
Izvedimo še aproksimacijo:
spl_a = UnivariateSpline(x, y, s=.1)
plt.plot(x, y, 'o', ms=5, label='Aproksimacijske točke')
plt.plot(x_g, spl(x_g), lw=3, label='Interpolacija (s=0)');
plt.plot(x_g, spl_a(x_g), lw=3, label='Aproksimacija (s=2)');
plt.plot(spl_a.get_knots(), spl_a(spl_a.get_knots()), 'o', label='Vozli pri aproksimaciji');
plt.legend();
Dejanski preostanek:
spl_a.get_residual()
0.0999987185545532
Vprašanja za vaje#
Vprašanje 1: Prihaja deževno vreme in zanima vas koeficient trenja med avtomobilsko gumo in mokro cesto. V garaži najdete odsluženo letno pnevmatiko, iz katere izrežete vzorec za testiranje. Nanj pritrdite silomer iz katerega odčitavete silo pri vlečenju gume po mokrem asfaltu. Da bi vaša meritev bila zanesljivejša, vzorec gume obtežujete z različnimi utežmi.
Veste, da je zveza med silo trenja ter normalno silo na podlago linearna:
Na podlagi podanih meritev želite določiti vrednost koeficienta tranja \(\mu\). Pomagate si z aproksimacijo (uporabite lahko poljubno metodo).
# Podatki
F_n = np.array([5, 10, 15, 20, 25, 30]) # N
F_t = np.array([ 1.432, 2.365, 3.91 , 5.166, 6.468, 7.438]) # N
Metoda najmanjših kvadratov za linearno funkcijo#
Vprašanje 2: Za podatke iz prejšnje naloge definirajte in rešite sistem linearnih enačb za izračun parametrov linearne aproksimacije. Rezultat prikažite grafično.
Vprašanje 3: Uporabite funkcijo numpy.polyfit in podane točke x, y aproksimirajte s polinomom prvega, drugega, tretjega in četrtega reda.
# Podatki, ne briši!
x = np.linspace(0, 5, 100)
y = 0.5*x**4 + 3*x**3 - 20*x**2 + 3*x + 10*np.random.randn(len(x))
Vprašanje 4: Pripravite funkcijo napaki(y, y_apr), ki izračuna in vrne vsoto kvadratov razlik ter standardno napako (numpy.std) za numerični polji podanih vrednosti y ter rezultata aproksimacije y_apr.
Uporabite jo za primerjavo aproksimacijskih polinomov različnih stopenj iz prejšnje naloge.
Vprašanje 5: Da bi določili koeficient zračnega upora \(C_D\) modelčka letala opravljate meritve v vetrovniku. Pri različnih hitrostih toka zraka večkrat pomerite silo zračnega upora \(F_D\). Veste, da je odvisnost sile \(F_D\) od hitrosti \(v\) v vašem primeru podana z enačbo:
Na podlagi podanih vrednosti v, F_d določite vrednost koeficienta zračnega upora s pomočjo aproksimacije s polinomom druge stopnje.
# Podatki
v = [6.184, 9.843, 15.582, 20.190, 23.509, 28.890, 34.753, 40.500, 45.979, 50.167]
F_d = [0.029, 0.013, 0.072, 0.054, 0.135, 0.204, 0.300, 0.420, 0.513, 0.643]
Vprašanje 6: Na istem grafu izrišite pomerjene točke, rezultat polinomske aproksimacije in rezultat izraza
iz prejšnje naloge.
Pojasnite zakaj prihaja do razlike med izrisanima krivuljama, kako to vpliva na določitev koeficienta \(C_D\) in kako bi to napako odpravili (namig: izpišite dobljene koeficiente aproksimacijskega polinoma).
Vprašanje 7: Z uporabo ustrezne scipy funkcije določite koeficient zračnega upora \(C_D\) tako, da pomerjene točke v, F_d aproksimirate s funkcijo:
Dobljeno vrednost primerjajte z rezultatom prejšnje naloge, rezultat grafično prikažite in komentirajte.
Vprašanje 8: Na počasnem posnetku smo opazovali odmik konice propelerja od horizontalne osi pri različnih časih. Na podlagi meritev želimo določiti hitrost vrtenja propelerja.
kjer je \(f\) število obratov propelerja v sekundi.
Določite \(f\) tako, da podatke aproksimirate s podano funkcijo in določite optimalne vrednosti parametrov \(f\) in \(\varphi\). Opazujte vpliv začetnih približkov p0, ki jih opcijsko lahko podamo funkciji scipy.optimize.curve_fit(fun, x, y, p0).
# Podatki
f = 50
phi = np.pi/4
t = np.linspace(0, 0.01, 100)
y = np.sin(2*np.pi*t*f + phi) + 0.1*np.random.randn(len(t))
Aproksimacija z zlepki#
Vprašanje 9: Podane vrednosti x, y aproksimirajte s kubičnimi zlepki z uporabo scipy.interpolate.UnivariateSpline.
Izrišite aproksimacijske krivulje pri treh podanih vrednostih parametra s. Za vsako izpišete število vozlov zlepka (spl.get_knots()) in komentirajte dobljeno.
# Podatki
x = np.linspace(-1, 1, 50)
y = 0.5*np.cos(1.7*np.pi*x) + 0.1*np.random.randn(len(x))
s = [0, 0.5, 2]