Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

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 (xx je raztezek, yy 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$');
<Figure size 640x480 with 1 Axes>

Za konkreten primer bi bilo, glede na poznavanje fizikalnega ozadja linearne vzmeti, primerno, da bi meritve poskušali popisati z linearno funkcijo:

f(x)=a0 x+a1f(x) = a_0\,x+ a_1
(1)

Poznamo tabelo nn podatkov xi,yix_i, y_i za i=0,1,…,n−1i=0,1,\dots,n-1; teh je več, kot jih potrebujemo za določitev dveh konstant a0a_0 in a1a_1, zato imamo torej predoločen sistem linearnih enačb:

yi=a0 xi+a1zai=0,1,…,n−1.y_i=a_0\,x_i+a_1\quad\textrm{za}\quad i=0,1,\dots,n-1.
(2)

Iščemo taki vrednosti konstanti a0a_0 in a1a_1, da se bo funkcija f(x)f(x) v znanih točkah xix_i najbolje ujemala z yiy_i.

Najprej torej potrebujemo kriterij za najboljše ujemanje.

Za vrednosti iz tabele xi,yix_i, y_i bi lahko iskali vrednosti a0a_0 in a1a_1, pri katerih bi bila vsota absolutne vrednosti odstopkov PP najmanjša:

P(a0,a1)=∑i=0n−1∣yi−(a0 xi+a1)∣.P(a_0, a_1) = \sum_{i=0}^{n-1} |y_i - (a_0\,x_i+a_1)|.
(3)

Ker pa taka funkcija P(a0,a1)P(a_0, a_1) ni zvezno odvedljiva, raje uporabimo metodo najmanjših kvadratov:

S(a0,a1)=∑i=0n−1(yi−(a0 xi+a1))2.S(a_0, a_1) = \sum_{i=0}^{n-1} \left(y_i - (a_0\,x_i+a_1)\right)^2.
(4)

Takšna funkcija S(a0,a1)S(a_0, a_1) je zvezna in zvezno odvedljiva. S parcialnim odvajanjem po parametrih a0a_0 in a1a_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)
Loading...

Poiskati moramo konstanti a0a_0, a1a_1, da bo vsota kvadratov razlik med funkcijo in tabelirano vrednostjo (xi,yix_i, y_i, kjer i=0,1,…,n−1i=0,1,\dots, n-1 in je nn število tabeliranih podatkov):

S(a0,a1)=∑i=0n−1(yi−(a0 xi+a1))2S(a_0, a_1) = \sum_{i=0}^{n-1} \left(y_i - (a_0\,x_i+a_1)\right)^2
(5)

najmanjša. Vrednost bo najmanjša v stacionarni točki, ki jo določimo s parcialnim odvodom po parametrih a0a_0 in a1a_1.

Najprej izvedemo parcialni odvod po parametru a0a_0:

∂S(a0,a1)∂a0=2 ∑i=0n−1(yi−a0 xi−a1) (−xi)\frac{\partial S(a_0, a_1)}{\partial a_0} = 2\,\sum_{i=0}^{n-1} \left(y_i - a_0\,x_i-a_1\right)\,(-x_i)
(6)

Izraz uredimo:

∂S(a0,a1)∂a0=−2(∑i=0n−1yi xi−a0 ∑i=0n−1xi2−a1 ∑i=0n−1xi)\frac{\partial S(a_0, a_1)}{\partial a_0} =-2\left(\sum_{i=0}^{n-1} y_i\,x_i -a_0\,\sum_{i=0}^{n-1} x_i^2- a_1\,\sum_{i=0}^{n-1} x_i\right)
(7)

Podobno postopamo še za a1a_1:

∂S(a0,a1)∂a1=2 ∑i=0n−1(yi−a0 xi−a1) (−1)\frac{\partial S(a_0, a_1)}{\partial a_1} = 2\,\sum_{i=0}^{n-1} \left(y_i - a_0\,x_i- a_1\right)\,(-1)
(8)
∂S(a0,a1)∂a1=−2(∑i=0n−1yi−a0 ∑i=0n−1xi−a1 ∑i=0n−11)\frac{\partial S(a_0, a_1)}{\partial a_1} = -2\left(\sum_{i=0}^{n-1} y_i -a_0\,\sum_{i=0}^{n-1} x_i- a_1\,\sum_{i=0}^{n-1}1\right)
(9)

Ker v stacionarni točki velja ∂S(a0,a1)/∂a0=0\partial S(a_0, a_1)/\partial a_0=0 in ∂S(a0,a1)/∂a1=0\partial S(a_0, a_1)/\partial a_1=0, iz zgornjih izrazov izpeljemo:

a0 ∑i=0n−1xi2+a1 ∑i=0n−1xi=∑i=0n−1yi xia_0\,\sum_{i=0}^{n-1} x_i^2 + a_1\,\sum_{i=0}^{n-1} x_i=\sum_{i=0}^{n-1} y_i\,x_i
(10)

in

a0 ∑i=0n−1xi+a1 n=∑i=0n−1yi.a_0\,\sum_{i=0}^{n-1} x_i+ a_1\,n=\sum_{i=0}^{n-1} y_i.
(11)

Dobili smo sistem dveh linearnih enačb za neznanki a0a_0 in a1a_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:

A a=b,\mathbf{A}\,\mathbf{a}=\mathbf{b},
(12)

Opomba, tukaj smo vektor neznank zapisali kot a=(a0,a1)\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(xi)f(x_i) od tabeliranih vrednosti yiy_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));
Loading...

Uporaba psevdo inverzne matrike

Do enakega rezultata lahko pridemo z uporabo psevdo inverzne matrike. Iščemo y(x)=a0 x+a1y(x)=a_0\,x+a_1 in nastavimo predoločen sistem A a=y\mathbf{A}\,\mathbf{a}=\mathbf{y}, kjer je matrika koeficientov A\mathbf{A} definirana glede na vrednosti xix_i (i=0,1,…,n−1i=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 b\mathbf{b}, v našem primeru pa je to kar vektor vrednosti y\mathbf{y} z elementi yiy_i (i=0,1,2,…,n−1i=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 mm linearnih enačb z nn neznankami ter velja m>nm>n in je rang razširjene matrike n+1n+1, imamo predoločeni sistem.

Vektor neznanih parametrov a\mathbf{a} določimo z uporabo psevdo inverzne matrike:

a=A+ y\mathbf{a}=\mathbf{A}^{+}\,\mathbf{y}
(13)
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 mm:

f(a0,a1,…,am,x)=∑s=0mas xm−s⏟fs(x),f(a_0,a_1,\dots,a_{m}, x) = \sum_{s=0}^{m}a_s\,\underbrace{x^{m-s}}_{f_s(x)},
(14)

kjer fs(x)=xm−sf_s(x)=x^{m-s} imenujemo bazna funkcija (s=0,1,2,…,ms=0,1,2,\dots,m).

Tabela podatkov naj bo definirana z xi,yix_i, y_i, kjer je i=0,1,2,…,n−1i=0,1,2,\dots,n-1.

Opomba: zaradi kompaktnosti zapisa bomo konstante aa zapisali v vektorski obliki a=[a0,a1,…,am]\mathbf{a}=[a_0, a_1,\dots,a_m].

Uporabimo metodo najmanjših kvadratov:

S(a)=∑i=0n−1(yi−f(a,xi))2=∑i=0n−1(yi−∑s=0mas xim−s)2.S(\mathbf{a}) = \sum_{i=0}^{n-1} \left(y_i - f(\mathbf{a}, x_i)\right)^2= \sum_{i=0}^{n-1} \left(y_i - \sum_{s=0}^{m}a_s\,x_i^{m-s}\right)^2.
(15)

Potreben pogoj za nastop ekstrema funkcije m+1m+1 neodvisnih spremenljivk je, da najdemo stacionarno točko za vsak ava_v, iščemo torej ∂S(a)/∂av=0\partial S(\mathbf{a})/\partial a_v=0 (namesto ss smo uporabili indeks vv).

Najprej določimo parcialni odvod za izbrani ava_v:

∂S(a)∂av=∑i=0n−1−2 (yi−∑s=0mas xim−s) xim−v\frac{\partial S(\mathbf{a})}{\partial a_v} = \sum_{i=0}^{n-1} - 2\,\left(y_i - \sum_{s=0}^{m}a_s\,x_i^{m-s}\right)\,x_i^{m-v}
(16)

Opomba: ∂∂av(∑s=0mas xim−s)=xim−v\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:

∑i=0n−1(∑s=0mas xim−s) xim−v=∑i=0n−1yi xim−v\sum_{i=0}^{n-1} \left(\sum_{s=0}^{m}a_s\,x_i^{m-s}\right)\,x_i^{m-v}=\sum_{i=0}^{n-1} y_i\,x_i^{m-v}
(17)

Izraz uredimo:

∑i=0n−1∑s=0mas xi2m−s−v=∑i=0n−1yi xim−v\sum_{i=0}^{n-1} \sum_{s=0}^m a_s\,x_i^{2m-s-v}=\sum_{i=0}^{n-1} y_i\,x_i^{m-v}
(18)

Zamenjamo vrstni red seštevanja ter izpeljemo:

∑s=0m(as∑i=0n−1 xi2m−s−v)=∑i=0n−1yi xim−vza:v=0,1,…,m\sum_{s=0}^m \left(a_s \sum_{i=0}^{n-1} \,x_i^{2m-s-v}\right)=\sum_{i=0}^{n-1} y_i\,x_i^{m-v}\quad\textrm{za:}\quad v=0,1,\dots,m
(19)

Izpeljali smo enačbo vv sistema m+1m+1 linearnih enačb:

A a=b\mathbf{A}\,{\mathbf{a}}=\mathbf{b}
(20)

Element Av,sA_{v,s} matrike koeficientov je:

Av,s=∑i=0n−1xi2m−v−s,A_{v,s}= \sum_{i=0}^{n-1} x_i^{2m-v-s},
(21)

Element vektorja konstant je:

bv=∑i=0n−1yi xim−vb_{v}= \sum_{i=0}^{n-1} y_i\,x_i^{m-v}
(22)

Numerični zgled

Uporabimo podatke iz prve naloge in poskusimo aproksimirati s polinomom 2. stopnje (m=2m=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:

Av,s=∑i=0n−1xi2m−v−sA_{v,s}= \sum_{i=0}^{n-1} x_i^{2m-v-s}
(23)
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:

bv=∑i=0n−1yi xim−vb_{v}= \sum_{i=0}^{n-1} y_i\,x_i^{m-v}
(24)
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:

f(a0,a1,…,am,x)=∑v=0mav xm−vf(a_0,a_1,\dots,a_{m}, x) = \sum_{v=0}^{m}a_v\,x^{m-v}
(25)

Kar v konkretnem primeru je aproksimacijski polinom:

f(x)=3.18393375 x2+14.64106847 x−1.22621065f(x)=3.18393375\,x^2 + 14.64106847\,x -1.22621065
(26)

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();
<Figure size 640x480 with 1 Axes>

Poglejmo še napako aproksimacije:

ei=yi−f(xi)e_i=y_i - f(x_i)
(27)

za i=0,1,2,…,n−1i=0,1,2,\dots,n-1.

Pri pravilno izvedeni aproksimaciji je nekaj eie_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 aproksimacijo 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).

Opomba: numpy.polyfit in numpy.poly1d spadata v starejši vmesnik za polinome; numpy za nove programe priporoča razred numpy.polynomial.Polynomial (dokumentacija): p = Polynomial.fit(x, y, deg=2) vrne objekt, ki ga kličemo kot funkcijo (p(x)), z p.convert().coef pa dobimo koeficiente v naraščajočem vrstnem redu potenc (ravno obratno kot polyfit).

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();
<Figure size 640x480 with 1 Axes>

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. a0 sin⁡(a1 x+a2)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 funkcijo, 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();
<Figure size 640x480 with 1 Axes>

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 amplituda, ω 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.01461945, 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();
<Figure size 640x480 with 1 Axes>

Dodatno

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.

Opomba: UnivariateSpline je od SciPy 1.15 označen kot legacy; sodobnejši ekvivalent za aproksimacijo (glajenje) z zlepki je scipy.interpolate.make_smoothing_spline (dokumentacija) oz. make_splrep. Objektni vmesnik spodaj še deluje in je preglednejši za učenje.

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();
<Figure size 640x480 with 1 Axes>

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=0.1)');
plt.plot(spl_a.get_knots(), spl_a(spl_a.get_knots()), 'o', label='Vozli pri aproksimaciji');
plt.legend();
<Figure size 640x480 with 1 Axes>

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:

Ft=μFnF_t = \mu F_n
(28)

Na podlagi podanih meritev želite določiti vrednost koeficienta trenja μ\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 CDC_D modelčka letala opravljate meritve v vetrovniku. Pri različnih hitrostih toka zraka večkrat pomerite silo zračnega upora FDF_D. Veste, da je odvisnost sile FDF_D od hitrosti vv v vašem primeru podana z enačbo:

FD=0.005 CD⋅v2F_D = 0.005~C_D \cdot v^2
(29)

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

FD=0.005 CD⋅v2F_D = 0.005~C_D \cdot v^2
(30)

iz prejšnje naloge.

Pojasnite zakaj prihaja do razlike med izrisanima krivuljama, kako to vpliva na določitev koeficienta CDC_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 CDC_D tako, da pomerjene točke v, F_d aproksimirate s funkcijo:

FD=0.005 CD⋅v2F_D = 0.005~C_D \cdot v^2
(31)

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.

y=sin⁡(2π⋅f t+φ)y = \sin(2\pi\cdot f~t + \varphi)
(32)

kjer je ff število obratov propelerja v sekundi.

Določite ff tako, da podatke aproksimirate s podano funkcijo in določite optimalne vrednosti parametrov ff 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]
Numerične metode v ekosistemu Pythona
Interpolacija
Numerične metode v ekosistemu Pythona
Reševanje enačb