Loeng 15: Teadusarvutused Pythonis, teek matplotlib ja mitme 2D teljepaari loomine joonise kohta, jooniste vormindamine, teekide NumPy ja SciPy kasutamise valitud näiteid

Viimati uuendatud 9.12.2024.

---

Koodinäited

Matplotlib keerukam süntaks ja jooniste vormindamine

Meeldetuletus

Eelmise nädala Loengus 14 tutvustasime kuidas luua jooniseid kasutades allolevat lihtsustatud ja kompaktset süntaksit:

import matplotlib.pyplot as plt
import numpy as np
%config InlineBackend.figure_format='retina'  # Vajadusel sisesta Spyderi konsooli.

x = np.linspace(0, 2, 100)
y = np.sin(2 * np.pi * x)

plt.plot(x, y, label='Legend 1')  # Argument label on legend.

plt.legend(loc=1)  # Legendi asukoht.
plt.xlabel('Argument $x$')  # '$' on LaTeXi exit tag.
plt.ylabel('Funktsioon $y(x)$')
plt.ylim(-1, 1.5)
plt.title('Pealkiri $\sin(2 \pi x)$')

plt.savefig("test.pdf")  # Joonise faili salvestamine töökausta.

plt.show()

Teegi Matplotlib keerukam süntaks

Üks teljepaar

Lisaks eelmisel nädalal tutvustatud süntaksile saame ühe ja mitme teljepaariga 2D graafikuid luua kasutades järgmist töövoogu ja süntaksit:

• Joonise tekitamine
• Telgede lisamine joonisele, telgede konfigureerimine
• Graafikute lisamine telgedele
• Joonise salvestamine ja kuvamine

x = np.linspace(0, 2)  # Vaikimisi, num = 50.
y = 1 + np.sin(2* np.pi * x)

fig, ax = plt.subplots()  # Joonise fig ja teljepaari ax loomine.

ax.plot(x, y, label='legend 1')  # Graafiku loomine ja lisamine teljepaarile.
ax.legend(loc=9)  # Legendi asukoht.
ax.grid()  # Ruudustiku lisamine.
ax.set(xlabel='Argument $x$',  # Telgede konfigureerimine.
       ylabel='Funktsioon $y(x)$',
       ylim=(0, 2.5),
       aspect=1,
       title='Pealkiri $y(x) = \int_{-a}^{x}(t-b)^2 dt$')  # Pealkiri, '$' on LaTeXi exit tag.

fig.savefig("test.pdf")  # Joonise salvestamine.

plt.show()  # Joonise kuvamine konsoolis või mujal.

Mitu teljepaari joonise kohta

Kui soovid joonisele lisada kaks 2D graafikut, s.t. kaks teljepaari:

x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)  # Joonise ja telgede loomine.
                            # 2 rida 1 veerg, jagatud x-telg.
ax1.plot(x, y, label='Legend 1')
ax1.legend(loc=1)
ax1.set(xlabel='x',
        ylabel='Funktsioon $y(x)$',
        title='Graafiku pealkiri')

ax2.scatter(x, y, label='Legend 2')  # Graafiku tüüp, saab asendada plot'ga.
ax2.legend(loc=3)
ax2.set_xlabel('x')  # Nii saab ka, vrd. ax1.
ax2.set_ylabel('y')  # Nii saab ka.

fig.suptitle('Joonise pealkiri')  # Kogu joonise pealkiri.
fig.savefig('test.pdf')

plt.show()

Kaks teljepaari paigutatuna ühte ritta:

x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)

fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True)  # 1 rida 2 veergu, jagavad y-telgi.

ax1.plot(x, y, label='Legend 1')
ax1.legend(loc=0)  # Automaatne parim paigutus.
ax1.set(xlabel='x',
        ylabel='Funktsioon $y(x)$',
        title='Graafiku pealkiri')

ax2.scatter(x, y, label='Legend 2')
ax2.legend(loc=3)
ax2.set_xlabel('x')  # Nii saab ka, vrd. ax1.

fig.suptitle('Joonise pealkiri')
fig.savefig('test.pdf')

plt.show()

Kolm teljepaari paigutatuna ühte veergu:

x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)  # 3 rida 1 veerg, jagavad x-telgi.

ax1.plot(x, y, label='Legend 1')
ax1.legend(loc=0)  # Automaatne parim paigutus.
ax1.set(xlabel='x',
        ylabel='y',
        title='Graafiku pealkiri')

ax2.plot(x, y, label='Legend 2')
ax2.legend(loc=2)
ax2.set(xlabel='x',
        ylabel='y')

ax3.scatter(x, y, label='Legend 3')   
ax3.legend(loc=3)
ax3.set_xlabel('x')

fig.suptitle('Joonise pealkiri')
fig.savefig('test.pdf')

plt.show()

Teljepaari objektile viitamine

Lisaks nimepidi muutujale viitamisele, nii nagu on tehtud eelmises koodirakus, saame teljepaarile viidata kasutades kandilisi sulge ja indekseid. Loodavatele teljepaaride objektidele saab viidata kasutades indekseid järgnevalt:

x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)

fig, axes = plt.subplots(1, 2, sharey=True)

axes[0].plot(x, y)
axes[0].set(ylabel='y',
            title='Graafiku pealkiri')

axes[1].plot(x, y)
axes[1].set(title='Jagame y-telgi')

fig.suptitle('Joonise pealkiri')
fig.savefig('test.pdf')

plt.show()

x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)

fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)  # 2 rida ja 2 veergu.

axes[0, 0].plot(x, y)
axes[0, 0].set(ylabel='y',
               title='Graafiku pealkiri')

axes[0, 1].plot(x, y)
axes[0, 1].set(title='Jagame x ja y-telgi')

axes[1, 0].remove()  # Eemaldab teljepaari.

axes[1, 1].scatter(x, 2*y)
axes[1, 1].set(xlabel='x')

fig.suptitle('Joonise pealkiri')
fig.savefig('test.pdf')

plt.show()

Erinevad graafikute tüübid

Matplotlib toetab erinevaid graafikute tüüpe. Teegi Matplotlib graafikute põhitüübid: https://matplotlib.org/stable/plot_types/index.html
Lisainfo ja graafikute näiteid Matplotlib kodulehelt: https://matplotlib.org/stable/gallery/index.html

Näiteks tutvustame graafikutüüpi mis võimaldab lisada määramatusi kujutatud andmepunktidele:

x = np.arange(10)
y = np.sqrt(x)

fig, axes = plt.subplots(1, 3, sharey=True)

axes[0].plot(x, y, 'o-')
axes[0].set(xlabel='x',
            ylabel='Päkapiku jume y(x)',
            aspect=2)

axes[1].errorbar(x, y,
                 xerr=0.6,
                 yerr=0.5 * y,
                 fmt=' ',
                 ecolor='k',
                 elinewidth=0.65,
                 capsize=3)
axes[1].set_xlabel('Juhuu $x$')
axes[1].set_aspect(4)

axes[2].errorbar(x, y,
                 xerr=0.6,
                 yerr=0.5 * y,
                 fmt=' ',
                 ecolor='k',
                 elinewidth=0.65,
                 capsize=3)
axes[2].plot(x, y, 'ko-')
axes[2].set_xlabel('x')

fig.savefig('test.pdf')

plt.show()

Kaks vertikaalset skaalat ühel graafikul

Kui on vajadus kasutada kahte skaalat korraga:

t = np.arange(0.01, 10, 0.01)

data1 = np.exp(t)
data2 = np.sin(t)

fig, ax1 = plt.subplots()

ax1.set_xlabel('Jagatud argumendid, Aeg $t$')
ax1.set_ylabel('Funktsioon $\exp(t)$', c='r')  # Või color='red'.
ax1.plot(t, data1, color='red')
ax1.tick_params(axis='y', labelcolor='red')

ax2 = ax1.twinx()  # Teine telg mis jagab argumentide-telje väärtusi.
ax2.set_ylabel('Funktsioon $\sin(t)$', color='blue')  # Teljepaar ax1 on konfigureeritud ülal.
ax2.plot(t, data2, color='blue')
ax2.tick_params(axis='y', labelcolor='blue')

fig.tight_layout()  # Kompaktsem joonise osade paigutus.
fig.savefig('test.pdf')

plt.show()

Jooniste vormindamine ja sõnaraamat matplotlib.pyplot.rcParams

Teegi Matplotlib jooniste loomisel kasutatud parameetrite ja muude muutujate väärtusi hoitakse sõnaraamatus matplotlib.pyplot.rcParams:

#print(plt.rcParams)  # Vaata sisu iseseisvalt.

Nii nagu ikka sõnaraamatutega saame sõnaraamatu matplotlib.pyplot.rcParams võtmete väärtusi meelevaldselt üle kirjutada. Muudame mitut väärtust korraga:

params = {'axes.labelsize': 15,
          'font.size': 10,
          'legend.fontsize': 9,
          'xtick.labelsize': 7,
          'ytick.labelsize': 7,
          'axes.titlepad': 1,
          'font.family': ['sans serif'],
          'figure.subplot.hspace': 0.1,
          'figure.subplot.wspace': 0.6,
          'figure.figsize': [7.9, 2.5]}  # Tollides.

plt.rcParams.update(params)  # Uuendame sõnaraamatut.

x = np.linspace(0, 2*np.pi, 300)
y = np.sin(x ** 2)

fig, axes = plt.subplots(1, 2)

axes[0].plot(x, y)
axes[0].set(ylabel='Funktsioon y',
            xlabel='Argument x',
            title='Pealkiri 1')

axes[1].plot(x, y ** 3, 'k')
axes[1].set(ylabel='y',
            xlabel='Argument x',
            title='Pealkiri 2',
            aspect=1)

fig.suptitle('Joonise pealkiri', y=0.87)  # Pealkirja asukoht y.
fig.tight_layout()
fig.savefig('test.pdf')

plt.show()

Või muuda sõnaraamtu matplotlib.pyplot.rcParams võtmete väärtusi ühekaupa, neid üle kirjutades:

plt.rcParams['figure.figsize'] = [6, 1]  # Tollides (inch).

x = np.linspace(0, 2*np.pi, 600)
y = np.sin(x ** 2.5)

plt.plot(x, y, lw=0.7)  # Eelmise nädala lühem süntaks.

plt.xlabel('Argument $x$')
plt.ylabel('$y(x)$')

plt.show()

Kui soovid parameetrid alglähtestada saad kasutada sõnaraamatut matplotlib.pyplot.rcParamsDefault:

plt.rcParams.update(plt.rcParamsDefault)

Teegi Matplotlib kasutamise näide

Histogramm ja logaritmilised teljed

Histogramm on statistiliste andmete esitamise viis mis annab ülevaate nende jaotusest sageduse järgi. Andmete väärtused on enamasti kujutatud horisontaalsel teljel ja andmete sagedused püstteljel.

Matploltib graafikutüüp matplotlib.pyplot.hist

n_bins = 100  # Histogrammi vahemike (bins) arv.
data = np.random.randn(10_000)  # Juhuslikud andmed, normaaljaotus.

fig, (ax1, ax2) = plt.subplots(1, 2)

ax1.plot(data, 'k-', lw=0.05)
ax1.set(xlabel='Massiivi indeks $i$',
        ylabel='Andmepunkti väärtus',
        title='Andmed')

ax2.hist(data, bins=n_bins, histtype='step', log=False)
ax2.set(xlabel='Andmepunkti väärtus',
        ylabel='Andmepunktide arv väärtuste vahemiku kohta',
        title='Histogramm')

fig.tight_layout()

plt.show()

Funktsioon numpy.histogram ja matplotlib.pyplot.plot

Eelmise koodirakuga sarnase tulemuse saame kui kasutame teegi NumPy algoritme:

n_bins = 100
data = np.random.randn(10_000)

hist = np.histogram(data, bins=n_bins)

fig, (ax1, ax2) = plt.subplots(1, 2)

ax1.plot(data, 'k-', lw=0.05)
ax1.set(xlabel='Massiivi indeks $i$',
        ylabel='Andmepunkti väärtus',
        title='Andmed')

ax2.plot(hist[1][1:], hist[0], '.-', lw=0.6)
ax2.set(xlabel='Andmepunkti väärtus',
        ylabel='Andmepunktide arv väärtuste vahemiku kohta',
        title='Histogramm')

fig.tight_layout()

plt.show()

Graafikutüübid matplotlib.pyplot.semilogy ja matplotlib.pyplot.loglog

Tihti on vaja andmeid esitada logaritmilistel skaaladel. Lisaks graafikutüübile matplotlib.pyplot.semilogy eksisteerib ka matplotlib.pyplot.semilogx. Kasutame eelmise koodiraku andmeid:

fig, (ax1, ax2) = plt.subplots(1, 2)

lw = 0.6
ax1.semilogy(hist[1][1:], hist[0], '.-', lw=lw)
ax1.set(xlabel='Andmepunkti väärtus',
        ylabel='Andmepunktide arv väärtuste vahemiku kohta',
        title='Histogramm 1')

ax2.loglog(hist[1][1:], hist[0], '.-', lw=lw)
ax2.set(xlabel='Andmepunkti väärtus',
        ylabel='Andmepunktide arv väärtuste vahemiku kohta',
        title='Histogramm 2')

fig.tight_layout()

plt.show()

Teekide NumPy ja SciPy kasutamise valitud näited

NumPy massiivi indekseerimine, funktsiooni tuletis ja vähimruutude meetod

Massiivi numpy.ndarray indekseerimine edasijõudnutele

Lisainfo: https://numpy.org/doc/stable/user/basics.indexing.html#advanced-indexing
Võrdlusoperaatorid indeksi kandilistes sulgudes. Alustame eelmise nädala süntaksi meeldetuletusega:

x = np.array([1, 2, 3, 4, 5, 6])

x % 2 == 0  # Element-element haaval ja edastatakse andmemassiivi numpy.ndarray.

array([False,  True, False,  True, False,  True])

NumPy andmemassiivi indekseerimine lubab rakendada tingimuslauseid massiivi elementidel järgmiselt:

x[x%2==0]  # Vali ainult need elemendid mis rahuldavad tingimust.

array([2, 4, 6])

Lisame täisarvu 20 ainult negatiivsetele massiivi elementidele:

x = np.array([1., -1., -2., 3])

x[x<0] += 20  # Ülekirjutatud operaator += itereerib üle massiivi elementide.
x

array([ 1., 19., 18.,  3.])

Asendame kuuest väiksemad elemendid arvuda 1000:

x = np.array([4, 5, 6, 7])

x[x<6] = 1000
x

array([1000, 1000,    6,    7])

Eelmised kolm näided põhinevad booli tõeväärtuste abil andmete indekseerimisele mida arutasime eelmine nädal:

x = np.array([1., -1., -2., 3])

y = x[[False, False, True, True]]  # Valime kasutades booli tõeväärtusi.
y

array([-2.,  3.])

Numbriline tuletise võtmine numpy.gradient

Võtame tuletise funktsioonist $\sin$ mille väärtused on salvestatud massiivi numpy.ndarray. Teatavasti: $$ \frac{{\rm d}}{{\rm d}x} \sin(x) = \cos(x). $$ Allolevates näidetes kasutame Matplotlibi lühemat süntaksit.

x = np.linspace(0, 2*np.pi, 300)

tuletis = np.gradient(np.sin(x), x)  # Samm dx = x[1]-x[0] jne on konstantne. Vrd järgmine koodirakk.

plt.plot(x, np.sin(x), label='Funktsioon $\sin(x)$')
plt.plot(x, tuletis, label='Tuletis $\cos(x)$')

plt.legend(loc=0)  # Automaatne paigutus.
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioonid $\sin(x)$ ja $\cos(x)$')

plt.show()

Sama mis eelmine näide, siin kasutame diferentsiaali dx väärtust massiivi x väärtuste asemel:

dx = 0.1  # Diferentsiaal.

x = np.arange(0, 2*np.pi, dx)
func = np.sin(x)

tuletis = np.gradient(func, dx)  # Eeldab, et dx on andmete ulatuses konstantne.

plt.plot(x, np.sin(x), label='Funktsioon $\sin(x)$')
plt.plot(x, tuletis, label='Tuletis $\cos(x)$')

plt.legend(loc=0)
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioonid $\sin(x)$ ja $\cos(x)$')

plt.show()

Vähimruutude meetod numpy.polyfit

Kasutame vähimruutude meetodit, et lähendada kolmanda astme polünoomi mis in kujul: $$ f(x) = a x^3 + b x^2 + c x + d, $$ kus $a$, $b$, $c$ ja $d$ on polünoomi koefitsendid, juhuslikult paigutatud ja juhuslikult kaalutud andmepunktidele. Allolevas koodis on ringi suurus (pindala) võrdeline andmepunkti kaaluga $-$ suurem ring suurem kaal.

n = 10

x = np.random.rand(n)  # Väärtused vahemikul [0, 1].
y = np.random.rand(n)

colors = np.random.rand(n)
area = (30 * np.random.rand(n)) ** 2

# Ruudu joonistamiseks:
x_0, y_0, x_1, y_1 = 0, 0, 1, 1  # Algandmete maksimaalne võimalik vahemik [0, 1].
ruut_koordinaadid = ([x_0, x_1, x_1, x_0, x_0], [y_0, y_0, y_1, y_1, y_0])

# Andmete lähendamine polünoomile:
kaal = area / 30**2  # Kaal on võrdeline ringi pindalaga.
z = np.polyfit(x, y, 3, w=kaal)  # Argument w on kaal.
pol = np.poly1d(z)  # Polünoomi objekt.

xp = np.linspace(0, 1, 100)

plt.plot(*ruut_koordinaadid, ls='--', lw=0.8, c='k')
plt.scatter(x, y, s=area, c=colors, alpha=0.6)
plt.plot(xp, pol(xp), lw = 2, label='VRM lähendus')

plt.xlabel('Koordinaat $x$')
plt.ylabel('Koordinaat $y$')
plt.title('Vähimruutude meetod')
plt.legend(loc=0)

plt.show()

Juured ehk nullkohad

Funktsiooni juurte otsimine scipy.optimize.root

Leiame avaldise juure. Avaldis on antud kujul: $$ x + 2 \cos(x). $$ Teatavasti juure leidmiseks peame lahendama võrrandi kujul: $$ x + 2 \cos(x) = 0. $$ Märkus: Sellel probleemil puudub lahend algebralistes funktsioonides (transtsendentsus).

from scipy.optimize import root

def func(x):
    return x + 2 * np.cos(x)

sol = root(func, 0.3)  # Väärtus 0.3 on algne oletus.

print(sol.x)        # Lahend.
print(sol.success)  # Kas õnnestus leida.

x = np.linspace(-3, 5, 300)

plt.axhline(y=0, c='k', lw=0.8)
plt.plot(x, func(x))  # Etteantud funktsiooni graafik.
plt.plot(sol.x, func(sol.x), 'o')  # Lahend = juur, üks punkt.

plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon $f(x)$')

plt.show()

[-1.02986653]
True

Kahe kõvera lõikepunktide leidmine scipy.optimize.root

Leia kahe polünoomiga defineeritud kõvera lõikepunktid vahemikul $x \in [-6, 6]$. Polünoomid on esitatud kujul: $$ y_1(x) = -2 x^2 + 2.56 x + 46, $$ $$ y_2(x) = 0.2 x^3 + 2.9 x^2 + 2.56x - 6. $$ Probleemi lahendamiseks tuleb lahendada võrrand kujul: $$ y_1(x) = y_2(x), $$ muutuja $x$ suhtes. Leitud muutujate $x$ väärtustele vastavad funktsiooni väärtused saab leida kasutades polünoomi $y_1(x)$ või $y_2(x)$.

from scipy.optimize import root

def y1(x): 
    return -2*x**2 + 2.56*x + 46

def y2(x): 
    return 0.2*x**3 + 2.9 * x**2 + 2.56*x - 6

def find_intersection(fun1, fun2, x0):
    probleem = lambda x: fun1(x) - fun2(x)
    tulem = root(probleem, x0)  # Algne oletus x0.
    return tulem.x

tulem = np.array([])  # Tühi massiiv.

for x0 in np.arange(-6, 6, 0.5):
    result = find_intersection(y1, y2, x0)
    tulem = np.append(tulem, result)

tulem = np.unique(tulem.round(decimals=10))

print('Lõikepunktide\nx koordinaadid on:', tulem)
print('y koordinaadid on:', y1(tulem))

x = np.linspace(-6, 6, 50)

plt.plot(x, y1(x), label='$y_1(x)$')
plt.plot(x, y2(x), label='$y_2(x)$')
plt.plot(tulem, y1(tulem), 'ro')  # Või y2(tulem)

plt.legend(loc=0)
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioonid $y_1(x)$ ja $y_2(x)$')

plt.show()

Lõikepunktide
x koordinaadid on: [-3.52036514  3.07087005]
y koordinaadid on: [12.20192387 35.00094163]

Kahe kõvera lõikepunktide leidmine scipy.optimize.fsolve

Lahendame eelmise näite kasutades teist lähenemist ja funktsiooni scipy.optimize.fsolve:

from scipy.optimize import fsolve

def y1(x):
    return -2*x**2 + 2.56*x + 46

def y2(x): 
    return 0.2*x**3 + 2.9*x**2 + 2.56*x - 6

def find_intersection(fun1, fun2, x0):
    probleem = lambda x: fun1(x) - fun2(x)
    return fsolve(probleem, x0)

tulem = np.array([])  # Tühi massiiv.

for x0 in np.arange(-6, 6, 0.5):
    result = find_intersection(y1, y2, x0)
    tulem = np.append(tulem, result)

tulem = np.unique(tulem.round(decimals=10))

print('Lõikepunktide\nx koordinaadid on:', tulem)
print('y koordinaadid on:', y1(tulem))

x = np.linspace(-6, 6, 50)

plt.plot(x, y1(x), label='$y_1(x)$')
plt.plot(x, y2(x), label='$y_2(x)$')
plt.plot(tulem, y1(tulem), 'ro')  # Või y2(tulem)

plt.legend(loc=0)
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioonid $y_1(x)$ ja $y_2(x)$')

plt.show()

Lõikepunktide
x koordinaadid on: [-3.52036514  3.07087005]
y koordinaadid on: [12.20192387 35.00094163]

Globaalne miinimum

Funktsiooni globaalse miinimumi leidmine scipy.optimize.minimize_scalar

Leiame vahemikus $x \in [-1, 3]$ plünoomi mis on esitatud kujul: $$ f(x) = (x - 2)(x + 1)^2 = x^3 - 3 x - 2, $$ miinimumile vastava muutuja $x$ väärtuse.

from scipy.optimize import minimize_scalar

f = lambda x: x**3 - 3*x - 2

res = minimize_scalar(f, method='brent')
print('Miinimumi asukoht: x =', res.x)

x = np.linspace(-1, 3, 300)

plt.plot(x, f(x))
plt.axhline(y=f(res.x), c='k', lw=0.8)
plt.plot(res.x, f(res.x), 'o')  # Lahend.

plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon $f(x)$')

plt.show()

Miinimumi asukoht: x = 1.0

Numbriline integreerimine

Määratud integraali leidmine scipy.integrate.quad

Leiame määratud integraali mis on esitatud kujul: $$ \int_0^4 \!\! x^2\, {\rm d}x. $$ Lisaks võrdleme tulemust analüütilise lahendiga kujul: $$ \int_0^4 \!\! x^2 \,{\rm d}x = \frac{~x^3}{3} \bigg\vert_0^4 = \frac{~4^3}{3} - 0 = \frac{~4^3}{3} = \frac{64}{3} = 21\frac{1}{3} = 21.(3). $$ Märkus: Nimi quad tuleneb programmeerimiskeele Fortran teegi QUADPACK nimest. Teek SciPy kasutab siin ajaloolisi QUADPACK algoritme.

from scipy.integrate import quad

x_ruudus = lambda x: x**2

lahend, err = quad(x_ruudus, 0, 4)  # Numbriline lahend.

print("Numbriline lahend on", lahend)
print("Analüütliline lahend on", 21 + 1 / 3)  # Analüütliline lahend.

x = np.linspace(-1, 4.1, 300)
x_area = np.linspace(0, 4, 300)

plt.fill_between(x_area, 0, x_ruudus(x_area), alpha=0.25)
plt.plot(x, x_ruudus(x))

plt.ylim(0, 17)
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon $x^2$')

plt.show()

Numbriline lahend on 21.333333333333332
Analüütliline lahend on 21.333333333333332

Määratud integraali leidmine scipy.integrate.simpson

Kasutame eelmise näite integraali ja selle lahendit. Võtame integraali kasutades teist algoritmi (Simpsoni valem) ja funktsiooni scipy.integrate.simpson:

from scipy.integrate import simpson

x_ruudus = lambda x: x**2

x = np.linspace(0, 4, 16)
fx_data = x_ruudus(x)

lahend = simpson(fx_data, x=x)  # Numbriline lahend mis põhineb etteantud andmetel mitte funktsioonil.

print("Numbriline lahend on", lahend)  # Numbriline lahend.
print("Analüütliline lahend on", 21 + 1 / 3)  # Analüütliline lahend.

x_area = np.linspace(0, 4, 300)

plt.fill_between(x_area, 0, x_ruudus(x_area), alpha=0.25)
plt.plot(x, fx_data, 'o')

plt.ylim(0, 17)
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon $x^2$')

plt.show()

Numbriline lahend on 21.333333333333332
Analüütliline lahend on 21.333333333333332

Päratud integraali leidmine scipy.integrate.quad

Võtame numbriliselt päratu integraali mis on antud kujul: $$ \int_0^\infty \!\!\!{\rm e}^{-x}\, {\rm d}x \equiv \int_0^\infty \!\!\!\exp(-x)\, {\rm d}x, $$ kusjuures anlüütline lahend on kujul: $$ \int_0^\infty \!\!\!{\rm e}^{-x} \,{\rm d}x \equiv \int_0^\infty \!\!\!\exp(-x)\, {\rm d}x = 1. $$

from scipy.integrate import quad

invexp = lambda x: np.exp(-x)

lahend, err = quad(invexp, 0, np.inf)

print("Numbriline lahend on", lahend)

x = np.linspace(0, 8, 300)

plt.fill_between(x, 0, invexp(x), alpha=0.25)
plt.plot(x, invexp(x))

plt.ylim(0, 1)
plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon e$^{-x}$')

plt.show()

Numbriline lahend on 1.0

Kumulatiivne integreerimine scipy.integrate.cumulative_trapezoid

Uurime sirgjoone $f(x) = x$ integraali lahendit $F(t)$, mis omab kuju: $$ F(x) = \int\! f(x)\, {\rm d}x = \int\! x\, {\rm d}x = \frac{~x^2}{2} + C, $$ käitumist vahemikul $x \in [-2, 2]$. Siin $C$ on integreerimiskonstant. Leiame integreerimiskonstandi väärtuse rajalt $x = -2$: $$ C = -\frac{(-2)^2}{~2} = -2. $$ Algoritm kasutab trapetsvalemit integrali võtmisel. Integraal lähendatakse trapetsite pindalade summale.

from scipy.integrate import cumulative_trapezoid

x = np.linspace(-2, 2, num=40)
fx = x  # Joone võrrand.

plt.fill_between(x, 0, fx, alpha=0.25)
plt.plot(x, fx)

plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon $f(x) = x$')

plt.show()

Fx = cumulative_trapezoid(fx, x, initial=0)  # Kumulatiivne integraal.

plt.plot(x, fx[0] + 0.5*x**2, 'b-', label='Analüütiline')
plt.plot(x, Fx, 'r.', label='Numbriline') 

plt.xlabel('Argument $x$')
plt.ylabel('Integraal $F(x) = \int_{-2}^x f(x)\,\mathrm{d} x = \int_{-2}^x x\,\mathrm{d} x$')
plt.legend(loc=0)

plt.show()

lahend, err = quad(lambda x: x, -2, 2)
print('Vastav vahemikul [-2, 2] määratud integraali väärtus on', lahend)

Vastav vahemikul [-2, 2] määratud integraali väärtus on 0.0

Kumulatiivne integreerimine scipy.integrate.cumulative_simpson

Kasutame eelmise näite integraali ja selle lahendit. Võtame kumulatiivselt integraali kasutades teist algoritmi (Simpsoni valem) ja funktsiooni scipy.integrate.cumulative_simpson:

from scipy.integrate import cumulative_simpson

x = np.linspace(-2, 2, num=40)
fx = x  # Joone võrrand.

plt.fill_between(x, 0, fx, alpha=0.25)
plt.plot(x, fx)

plt.xlabel('Argument $x$')
plt.ylabel('Funktsioon $f(x) = x$')

plt.show()

Fx = cumulative_simpson(fx, x=x, initial=0)  # Kumulatiivne integraal.

plt.plot(x, fx[0] + 0.5*x**2, 'b-', label='Analüütiline')
plt.plot(x, Fx, 'r.', label='Numbriline') 

plt.xlabel('Argument $x$')
plt.ylabel('Integraal $F(x) = \int_{-2}^x f(x)\,\mathrm{d} x = \int_{-2}^x x\,\mathrm{d} x$')
plt.legend(loc=0)

plt.show()

lahend, err = quad(lambda x: x, -2, 2)
print('Vastav vahemikul [-2, 2] määratud integraali väärtus on', lahend)

Vastav vahemikul [-2, 2] määratud integraali väärtus on 0.0

Algväärtusülesanne ehk Cauchy probleem

Algväärtusülesanne ja scipy.integrate.solve_ivp

Lahedame numbriliselt algväärtusülesande mille üldkuju on järgmine: $$ \begin{cases} \dfrac{{\rm d}}{{\rm d} t}x(t) = f(t, x(t)), \\ x(t_0) = x_0 . \end{cases} $$ Funktsioon $f$ on meelevaldne funktsioon, $x_0$ on algtingimuse väärtus, $t_0$ on alghetk, aeg $t$ on sõltumatu muutuja ja $t \in [t_0, t_\text{max}]$, kus $t_\text{max}$ on soovitud integreerimisaeg. Sõltuv muutuja on $x(t)$.

Näide: Lahenda algväärtusülesanne mis on antud autonoomsel kujul: $$ \begin{cases} \dfrac{{\rm d}x}{{\rm d} t} = \sin(x),\\ x(0) = 0.5, \end{cases} $$ sõltumatu muutuja $t$ on määratud vahemikus $t \in [0, 7]$.
Lugemist: https://en.wikipedia.org/wiki/Initial_value_problem

from scipy.integrate import solve_ivp

def func(t, x):
    return np.sin(x)  # Ei sõltu ajast ilmutatud kujul, autonoomne süsteem.

t0, t_max = 0, 7
x0 = 0.5   # Algtingimus.

sol = solve_ivp(func, [t0, t_max], [x0], dense_output=True)  # Saab ette öelda kus arvutada tulemusi, t_eval = t.

t = np.linspace(t0, t_max, 300)  # Ajahetked kus leiab väärtused.

x = sol.sol(t)

plt.plot(t, x.T, label='$x(t)$')

plt.xlabel('Sõltumatu muutuja, aeg $t$')
plt.ylabel('Lahend $x(t)$')
plt.legend(loc=0)
plt.title('Cauchy probleem')

plt.show()

Algväärtusülesanne ja scipy.integrate.odeint

Kasutame eelmise näite võrrandisüsteemi ja integreerime seda kasutades scipy.integrate.odeint:

from scipy.integrate import odeint

def func(x, t): 
    return np.sin(x)

t0, t_max = 0, 7
x0 = 0.5   # Algtingimus.

t = np.linspace(t0, t_max, 300)  # Ajahetked kus leiab väärtused.

sol = odeint(func, x0, t)

plt.plot(t, sol, label='$x(t)$')

plt.xlabel('Sõltumatu muutuja, aeg $t$')
plt.ylabel('Lahend $x(t)$')
plt.legend(loc=0)
plt.title('Cauchy probleem')

plt.show()

Lahendite parv mis koosned 40-st lahendist mis on leitud ühtlaselt jaotatud algtingimuste $x_0 = x(0)$ väärtuste jaoks. Algtingimused $x_0 \in [0, 4 \pi]$:

def func(x, t): 
    return np.sin(x)

n_sol = 40
d_n_sol = 4*np.pi / n_sol

t0, t_max = 0, 7
t = np.linspace(t0, t_max, 300)  # Ajahetked kus leiab väärtused.

for x0 in np.arange(0, 12.6, d_n_sol):
    sol = odeint(func, x0, t)
    plt.plot(t, sol)

plt.xlabel('Sõltumatu muutuja, aeg $t$')
plt.ylabel('Lahendid $x(t)$')
plt.title('Lahendite parv, $x_0 \in [0, 4\pi]$')

plt.show()

Konvolutsioon scipy.signal.convolve

Konvolutsioon on integraalteisendus kahe funktsiooni, nt. ühemõõtmeliste funktsioonide $f(t)$ ja $g(t)$ vahel, kus $t$ on funktsioonide sõltumatu muutuja. Konvolutsiooni tulemusena tekib kolmas funktsioon $(f * g)(t)$, mis kirjeldab kuidas ühe funktsiooni kuju muudab teist. Konvolutsiooni terminit kasutatakse nii tulemusfunktsiooni jaoks kui ka protsessi kirjelduseks. Tegemist on integraaliga funktsioonide $f(t)$ ja $g(t)$ korrutisest (integraalkorrutis), kus ühe funktsiooni argumenti on peegeldatud ja progresiivselt nihutatud.

Kahe etteantud ühemõõtmelise funktsiooni $f(t)$ ja $g(t)$ konvolutsioon on defineeritud järgmiselt: $$ (f * g)(t) = \int_{-\infty}^\infty \!\!f(\tau) \, g(t - \tau) \, {\rm d} \tau. $$ Vahelduseks kasutame allpool teegi Matplotlib pikemat süntaksit.
Lugemist: https://en.wikipedia.org/wiki/Convolution

from scipy.signal import convolve

n = 200  # Andmepunktide arv.
t_max = 1  # Maksimaalne aeg t-teljel.
t_pulse = 0.6  # Signaali pulsi asukoht aja teljel.
a = 900  # Signaali laiuse parameeter.

# Signaal.
t = np.linspace(0, t_max, n)  # Ajahetked.
sig = 0.5 * np.exp(-a * (t - t_pulse * t_max)**2)  # Signaal.
noise = np.random.rand(n) - 0.5
noisy_sig = sig + noise  # Müraga signaal.

# Tuum.
n_ker = n // 3
t_max_ker = t_max / 3
t_ker = np.linspace(0, t_max_ker, n_ker)
ker = np.exp(-a * (t_ker - 0.5 * t_max_ker)**2)

conv = convolve(noisy_sig, ker, mode='same')  # Konvolutsioon.
max_conv_i = np.argmax(conv)  # Asukoht vastavalt konvolutsiooni funktsioonile.
delta_t = np.abs(t_max*t_pulse - t[max_conv_i])  # Erinevus tegelikust asukohast.
print('Leitud asukoha erinevus tegelikust:', delta_t, 'aja ühikut.')

# Joonise loomine.
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)

ax1.plot(t, noisy_sig)
ax1.plot(t, sig)
ax1.axvline(t_max * t_pulse, lw=0.8, c='r', label='Asukoht.')
ax1.legend(loc=0)
ax1.set(ylabel='$f(t)$',
        title='Signaal ja signaal koos müraga')

ax2.plot(t_ker, ker)
ax2.set(ylabel='$g(t)$',
        title='Tuum')

ax3.plot(t, conv)
ax3.axvline(t_max * t_pulse, lw=0.8, c='r', label='Tegelik asukoht.')
ax3.axvline(t[max_conv_i], ls='--', c='g', label='Leitud asukoht.')
ax3.legend(loc=0)
ax3.set(xlabel='Aeg $t$',
        ylabel='$(f \\ast g)(t)$',
        title='Müraga signaali ja tuuma konvolutsioon')

fig.tight_layout()

plt.show()

Leitud asukoha erinevus tegelikust: 0.0030150753768845018 aja ühikut.

☻   ☻   ☻