Töölehe optimaalse kujunduse huvides veendu, et tööriistaribal on töölehe staatuseks märgitud **Trusted**. Kui nupukese pealkiri on **Not Trusted**, siis vajuta sellele.

Python ja Jupyter Notebook¶

Valter Kiisk
TÜ Füüsika Instituut
Viimati muudetud: 5.03.2024
$\newcommand{\eps}{\varepsilon}$

Sisukord

  • 1  Algteadmised
    • 1.1  Toimetamine Jupyteri töölehel
    • 1.2  Lihtsaimad arvutused
    • 1.3  Sõne
    • 1.4  Kommentaarid
  • 2  Tekst, joonised ja valemid
    • 2.1  Markdown
    • 2.2  LaTeX
  • 3  Andmetüübid
    • 3.1  Täisarvud
    • 3.2  Ujukomaarvud
    • 3.3  Kompleksarvud
  • 4  Programmeerimine
    • 4.1  IF-lause
    • 4.2  WHILE-tsükkel
    • 4.3  FOR-tsükkel
    • 4.4  Alamprogramm
  • 5  Jadad
  • 6  Graafikud
  • 7  NumPy
  • 8  Numbrilised algoritmid
  • 9  Failide lugemine-kirjutamine
  • 10  Sümbolarvutus
    • 10.1  Sümbolid ja avaldised
    • 10.2  Avaldiste teisendamine
    • 10.3  Matemaatiline analüüs
    • 10.4  Võrrandite lahendamine
    • 10.5  Lineaaralgebra

Jupyter Notebook on veebilehitsejas toimiv märkmeraamatu-tüüpi arvutuskeskkond. Sellises keskkonnas saab mitte ainult interaktiivselt arvutada, vaid ka arvutust dokumenteerida ning koostada mitmekülgse matemaatilise sisuga dokumente ja õppematerjale. Käesolevas juhendis eeldame, et arvutuskeeleks on Python, mis on hetkel populaarseim ja rikkaliku teadusliku ökosüsteemiga.

Jupyteri tööleht koosneb ühetaolistest kastikestest e. lahtritest (cell), mis on kõik töölehe laiused ja vertikaalselt rivistatud. Näiteks siinse kastikese sisuks on tavaline tekst, aga järgnev lahter on ette nähtud koodi käivitamiseks. Lihtsaim, mida Jupyter suudab teha, on kalkulaatori asendamine. Kirjuta järgmisse lahtrisse avaldis
4.3 + 1.75 * 2.1 / (8.6 - 5.1)
ja vajuta Ctrl+Enter!

In [ ]:
...

Kohe kastikese järel peaks tekkima arvutuse tulemus 5.35. Trükituna või paberile kirjutatuna näeks see avaldis välja nii: $$4,\!3 + \frac{1,\!75\times 2,\!1}{8,\!6 - 5,\!1}.$$ Muuda koodi lahtris ja veendu, et pärast Ctrl+Enter vajutamist tulemus vastavalt uueneb!

Algteadmised¶

Toimetamine Jupyteri töölehel¶

Märkmeraamatu ülaservas paiknevad tiitelriba, menüüriba ja tööriistariba. Tiitelribal kuvatakse faili nimi (muutmiseks kliki sellel). Ekraanipinna säästmiseks saab tiitelriba ära peita käsuga View | Show Header. (Vajadusel võib ka veebibrauseri viia täisekraani režiimi, F11).

Tööleht salvestub automaatselt iga paari minuti tagant. Töölehe koosseisus säilivad kõik lahtrid koos arvutustulemustega (sh graafikud). Kui salvestada tööleht käsuga Save and create checkpoint (esimene nupp tööriistaribal), luuakse ühtlasi ka varukoopia. Vajadusel saab sellest töölehe taastada käsuga File | Revert Notebook to Checkpoint.

Jupyteri töölehel on kaks erinevat seisundit. Redigeerimisrežiimis (edit mode) sisestatakse klahvivajutustega lihtsalt vastavaid sümboleid lahtrisse. Seda seisundit markeerib aktiveeritud koodiredaktor koos blinkiva kursoriga. Seevastu töölehe kui terviku kontrollimiseks on ette nähtud käsurežiim (command mode), kus isegi üksikutel klahvivajutustel võib olla eritähendus, näiteks lahtrite tekitamiseks, kustutamiseks või ümberpaigutamiseks (vt Help | Show Keyboard Shortcuts).

Uue lahtri loomiseks kliki esmalt hiirega selle lahtri päisel, mille ette või taha soovid uue lahtri tekitada (või, olles redigeerimisrežiimis, vajuta Esc). Tööleht läheb käsurežiimi. Nüüd vajuta klahvile A või B. Uus lahter vaikimisi ongi koodilahter (ehk tüübist Code, nagu näitab tekst tööriistaribal). Tekstilahtri saamiseks vajuta käsurežiimis klahvile M või vali tööriistaribalt lahtri tüübiks Markdown. Lahtri kustutamiseks vajuta (ikka käsurežiimis) piisavalt kiiresti kaks korda klahvile D.

Redigeerimisrežiimi aktiveerimiseks kliki hiirekursoriga koodilahtris (topelt-klõps tekstilahtris) või vajuta Enter. Tee seda näiteks kohe siinsamas, siis näed Markdown/HTML-koodi, mis määrab selle lahtri sisu.

Lihtsaimad arvutused¶

Matemaatilistest operatsioonidest on Pythonis koheselt olemas aritmeetikatehted ja astendamine. Suurima prioriteediga on astendamine (**), siis korrutamine-jagamine (*, /) ja viimaks liitmine-lahutamine (+, -). Arvutuste järjekorda saab muuta ümarsulgudega, nagu tegime eespool.

Paljudes süsteemides on astendamise operaatoriks katus (^), aga Pythonis on selleks kaks järjestikust korrutusmärki. Näiteks raadiusega 7,5 ühikut ringi pindala ($\pi r^2$) avaldub nii:

In [1]:
3.14 * 7.5**2
Out[1]:
176.625

(Siin ja edaspidi: aktiveeri lahter ja vajuta Ctrl+Enter.)

Täiendavaid matemaatilisi vahendeid leiab mitmesugustest lisamoodulitest. Elementaarfunktsioonid ja matemaatika põhikonstandid pi ja e saadakse moodulist math. Mooduli kasutamiseks tuleb see esmalt importida:

In [2]:
import math

Võtame ümberpööratud ülesande: teame ringi pindala 176,6 ja tahame selle kaudu leida ringi raadiuse. Selleks on tarvis ruutjuure arvutamise funktsiooni sqrt. Tuleb spetsiaalselt näidata, et me võtame selle funktsiooni ja konstandi $\pi$ moodulist math:

In [3]:
math.sqrt( 176.6 / math.pi )
Out[3]:
7.497567999028581

sqrt ja pi on muutujad, mis on defineeritud moodulis math. Need võib ka mooduli seest välja tuua (globaalsesse nimeruumi), et edaspidi poleks eesliidet math. enam tarvis:

In [4]:
from math import sqrt, pi
sqrt( 176.6 / pi )
Out[4]:
7.497567999028581

Vajadusel võib moodulist korraga kõik objektid välja tuua:

In [5]:
from math import *

Mitme mooduli üheaegsel kasutamisel seda sageli ei tehta, et vältida nimekonflikte. Näiteks moodulid math, cmath, mpmath, numpy ja sympy sisaldavad nime poolest peaaegu identse komplekti matemaatikafunktsioone erinevat tüüpi andmete jaoks.

Nagu näha, Jupyter automaatselt nummerdab kõik koodilahtrid. Lahtri käivitamisel saadud tulemus jäetakse meelde ja selle poole saab edaspidi pöörduda koodiga Out[n] või lihtsalt _n, kus n asemel on lahtri järjekorranumber (allkriips üksinda tagastab viimase arvutustulemuse). Arvutuse selguse ja dokumenteerimise huvides on mõistlikum anda nii algandmetele kui ka arvutustulemustele väljendusrikkad nimed:

In [6]:
pindala = 176.6
raadius = sqrt( pindala / pi )
raadius
Out[6]:
7.497567999028581

Sellega me defineerisime kaks uut muutujat pindala ja raadius. Muutuja nimi võib olla meelevaldne tähtede ja numbrite kombinatsioon (sh täpitähed ja kreeka sümbolid), kuid ei tohi sisaldada tühikuid ega alata numbriga. Pythoni võtmesõnade (if, for, def, jne) kasutamine on samuti keelatud. Kreeka tähestiku sümboleid saab tekitada otse Jupyteri keskkonnas sisestades esmalt vastava LaTeX'i koodi (nt \alpha, \delta, \lambda, \Sigma jne) ja vajutades seejärel tabulatsiooniklahvi. Pythoni kood on tõstutundlik (eristatakse suur- ja väiketähti). Pythoni stiilireeglid näevad ette, et muutujanimed võiksid olla väiketähtedega, kus sõnad eraldatakse allkriipsuga (nt raadius või ringi_raadius).

Pythonis iga lause (st üks terviklik programmikäsk) kirjutatakse omaette reale. Laused täidetakse järjest ülevalt alla liikudes. Seega esmalt antakse arvule 176.6 nimeks pindala, seejärel arvutatakse välja avaldise sqrt( pindala / pi ) väärtus ja sellele antakse omakorda nimeks raadius. Nüüd saab muutujaid raadius ja pindala kasutada edasiste arvutuste tegemiseks mistahes koodilahtris. Kui on soov alustada arvutustega otsast peale ja vabaneda kõigist globaalsetest muutujatest (sh muutujatest/moodulitest, mis on imporditud), tuleb arvutusmootorile teha alglaadimine (Kernel | Restart Kernel).

Harjutus. Alghetkel on pangakontol €1000. Iga aasta lisanduvad säästud €500. Hoiustamisest saadav aastane intress on 3,0% aasta alguses eksisteerivast rahasummast. Leia 7 aastaga teenitav intressitulu! Liitintressi $r$ korral on konto seis pärast $n$ perioodi möödumist $$K(n)=K_0(1+r)^n+S\frac{(1+r)^n-1}{r},$$ kus $K_0$ on algkapital ja $S$ on sissemakse perioodi jooksul. Arvutuse selguse huvides säilita eraldi muutujates nii algandmed kui ka vahetulemus $(1+r)^n$.

Trigonomeetrilised funktsioonid eeldavad alati argumenti radiaanides. Kraadide kasutamiseks on tarvis teisendustegurit, mis esindaks ühe kraadi väärtust radiaanides. Teades et 180 kraadi on võrdne $\pi$ radiaaniga, defineeri vastav teisendustegur (muutuja deg) ja seejärel veendu, et sin(30°) annab oodatava tulemuse 0,5:

In [ ]:
from math import sin, pi

deg = ...
sin(30*deg)

See viib meid üldisema probleemi juurde: loodusteaduslike või insenertehniliste rehkenduste juures läheb enamasti tarvis nii ühikute teisendamist kui ka mitmesuguseid füüsikalisi konstante. Vajadusel saab neid muidugi ise defineerida äsjamainitud viisil, kuid põhilised neist on juba defineeritud moodulis scipy.constants (väljendatuna SI põhiühikutes). Konstantide täielik nimekiri on siin.

In [7]:
from scipy.constants import m_e, h, eV, e, epsilon_0 as eps
m_e * e**4 / (8 * h**2 * eps**2) / eV
Out[7]:
13.605693122885842

Importimisel andsime muutujale epsilon_0 lühema nime eps. Toodud avaldis annab elektroni seoseenergia vesinikuaatomis elektronvoltides. Vastav õpikuvalem (SI-süsteemis) näeks välja nii: $$E=\frac{m_e e^4}{8h^2\varepsilon_0^2}.$$ Märgime, et siin Python opereerib vaid dimensioonitute arvudega, st sisulist ühikute teisendamist/dimensioonide kontrolli ei toimu. Seepärast on mõistlik arvutused teostada kindlate, omavahel kooskõlaliste mõõtühikute süsteemis. Näiteks, kui kõik arvandmed esitada SI põhiühikutes ja valemid panna samuti kirja SI-süsteemis, siis ilmselt ka arvutustulemused saavad olema SI põhiühikutes. Ühikute teisendamiseks tuleb kasutada lihtsalt sobivaid teisendustegureid. Viimastes näidetes muutujaga deg korrutamisel teisenesid kraadid radiaanideks ja muutujaga eV jagamisel teisenesid džaulid elektronvoltideks, st radiaane ja džaule käsitleti kui põhiühikuid.

Väga suured või väga väikesed arvud (nagu 1,6×10-19) tuleb programmikoodis kirja panna liht-tekstina kujul 1.6e-19 (kümnenderaldaja on alati punkt).

Harjutus. Arvuta raskuskiirendus Maa pinnal valemiga $g=GM/R^2$, kus Maa mass $M=5,\!97\times 10^{24}$ kg ja (keskmine) raadius $R=6371$ km. Gravitatsioonikonstandi saab moodulist scipy.constants.

Sõne¶

Kui arvutuslahtri viimane lause tagastab mingisuguse andmeobjekti (mida ei omistata muutujale), siis see trükitakse teksti kujul ekraanile (kohe lahtri järel). Seejuures reaalarvulised tulemused väljastatakse maksimaalse täpsusega (umbes 17 tüvenumbrit). Mõnikord on meil tarvis:

  • väljastada ühest lahtrist mitu (vahe)tulemust
  • väljastada peale arvude ka muud informatsiooni teksti kujul
  • kontrollida arvuliste tulemuste formaati ja täpsust

Esimese probleemi lahendab käsk print:

In [8]:
print(pindala)
print(raadius)
176.6
7.497567999028581

Teise probleemi lahendamiseks on tarvis uut andmetüüpi, mille nimi on string ehk eestipäraselt sõne. Sõne-tüüpi objekti sisuks on meelevaldne kirjasümbolite jada. Sõneliteraal piiritletakse ülakomade või jutumärkidega. print-käsule võib argumentidena anda hulga arvulisi või sõne-tüüpi andmeid (muutujaid, literaale), mis trükitakse ekraanile üksteise järel tühikutega eraldatult. Niiviisi saame konstrueerida tervikliku lause:

In [9]:
print( 'ringi raadius on', raadius, 'ning pindala', pindala )
ringi raadius on 7.497567999028581 ning pindala 176.6

Viimaks, spetsiaalse vorminduskoodi abil saame täpselt kontrollida arvude formaati ja täpsust. Vorminduskoodid kaasatakse esmalt sõne koosseisu ja operaatori % abil asendatakse samas järjekorras vastavate arvandmetega:

In [10]:
print( 'ringi raadius on %.3f ning pindala %.1f' % (raadius, pindala) )
ringi raadius on 7.498 ning pindala 176.6

Veelgi mugavam on selline süntaks:

In [11]:
print( f'ringi raadius on {raadius:.3f} ning pindala {pindala:.1f}' )
ringi raadius on 7.498 ning pindala 176.6

(Pane tähele f-tähte sõne ees!) Koodiga .3f vormindatakse reaalarv püsikomaarvuna 3 komakoha täpsusega. Seevastu koodiga .3e vormindatakse reaalarv ujukomaarvuna ehk eksponentkujul 3 komakoha täpsusega. Viimaks, .4g valib kompaktseima esituse, säilitades vähemalt 4 tüvenumbrit (v.a. nullid arvu lõpus).

In [12]:
print('%.4f' % raadius)
print('%.4e' % raadius)
print('%.4g' % raadius)
7.4976
7.4976e+00
7.498

Selliste 1-mastaabis arvude korral ilmselt eksponentkuju suurt väärtust ei oma (lisab alati aru lõppu e+00), optimaalne on püsikomaformaat. Seevastu väga suuri või väga väikeseid arve pole reeglina mõttekas esitada püsikomaarvuna (v.a. finantsarvutused):

In [13]:
x = 2.641e8
print('%.4f' % x)
print('%.4e' % x)
print('%.4g' % x)
264100000.0000
2.6410e+08
2.641e+08

Ilmselgelt arvu 2,6410×108 täpsus ei ole piisav õigustamaks vormingut 264100000.0000. Seevastu formaat %.4e võib olla vägagi sobilik, andes mõista, et kõik 4 kohta pärast koma on tähendusega ja selle täpsuse piires tuleks säilitada isegi nullid arvu lõpus.

Harjutus. Kirjuta programm, mis väljastab ruutvõrrandi (näiteks $2x^2+3x-35=0$) mõlemad lahendid 5 tüvenumbri täpsusega! Algandmeteks on ruutvõrrandi kordajad (näiteks a, b, c).

Jupyteri keskkonnas saab kergesti küsida infot mistahes (dokumenteeritud) elemendi kohta Pythoni koodis (olgu see muutuja, funktsioon, meetod, moodul, vm). Selleks tuleb kursoriga liikuda elemendi peale ja vajutada seejärel Shift+Tab. Näiteks tehes seda funktsiooninimel print, ilmub selle kohale lipik teabega

No description has been provided for this image

Siit saame teada, et sellel funktsioonil on parameetrid sep ja end. Esimene määrab sõne, millega eraldatakse ekraanile trükitavad väärtused (vaikimisi tühik). Parameetri end vaikeväärtus on '\n', mis tekitab reavahetuse, nii et kui mitu print-käsku on järjest, siis kõik väljastused on eraldi real. On hulk kirjasümboleid, mida ühel või teisel põhjusel ei saa otse sõne koosseisus näidata. Neid tähistatakse spetsiaalse koodiga ehk paojadaga (escape sequence), mis algavad langkriipsuga (\). Tuntumad näited on reavahetuse sümbol (\n), tabulatsioonisümbol (\t), jutumärgid (\"), apostroof (\') ja langkriips ise (\\).

No description has been provided for this imageKatseta ka sellist tegevust: kirjuta koodilahtrisse math. ja vajuta klahvile Tab. Tulemusena näidatakse nimekirja kõigist funktsioonidest ja muutujatest mooduli math koosseisus. Niiviisi saab teada ka objekti meetodid ja atribuudid. Kui sooritada see tegevus näiteks sõnemuutujaga, saame teada, et sõnel on hulk kasulikke meetodeid: split, replace, upper, jne. Klahvivajutusega Tab saab teada ka globaalsed muutujanimed. Näiteks, kui sisestada täht p ja vajutada Tab, näidatakse kõikide p-tähega algavate muutujate nimekirja, mille hulgas on ootuspäraselt pi, pindala, pow ja print.

Kommentaarid¶

Kui ühte lahtrisse mahutada pikem arvutuskood (terve algoritm või programm), on mõistlik ka otse koodi sees üht-teist kommenteerida. Pythonis üherealine kommentaar algab sümboliga # ja kulgeb rea lõpuni. Muidugi iga triviaalset liigutust ei ole mõtet kommenteerida — kommenteerimine peab tegema koodi selgemaks, mitte vastupidi. Kommenteerimise vajalikkus sõltub oluliselt näiteks muutujanimede valikust. Seega järgnevas koodis mõned piisavalt ilmsed operatsioonid on jäetud kommenteerimata.

In [14]:
from math import sin, cos
from scipy.constants import degree as deg, g

# algandmed
kiirus = 62 # algkiirus [m/s]
nurk = 45 # viskenurk kraadides

# arvutus
nurk = nurk * deg # radiaanideks
kiirus_x = kiirus * cos(nurk)
kiirus_y = kiirus * sin(nurk)
aeg = 2 * kiirus_y / g # lennuaeg
kaugus = kiirus_x * aeg
print(f'lennukaugus = {kaugus:.1f} m')
lennukaugus = 392.0 m

Üht ja sama muutujanime võib korduvalt kasutada, kusjuures muutuja uue väärtuse võib defineerida ka kohe selle muutuja olemasoleva väärtuse baasil (nagu lauses nurk=nurk*deg), sest programmeerimises võrdusmärk tähendab omistamist, mitte vasaku ja parema poole võrdsust. Seda tüüpi tehteid läheb sageli vaja, selle tõttu operatsioon x = x * a on lubatud kirjutada lühemalt kujul x *= a. Ülejäänud analoogsed operaatorid on +=, -= ja /=.

Tekst, joonised ja valemid¶

Markdown¶

Koodilahtrite vahele saab kirjutada rikkaliku kujundusega selgitavat teksti (nagu käesolev lõik), kasutades lihtsat märgistuskeelt nimega Markdown. Tekstiala tekitamiseks tuleb luua uus lahter (B), valida selle tüübiks Markdown (M), kirjutada sellesse Markdown-koodis märgendatud sisu ja viimaks lasta see kujundada klahvivajutusega Ctrl+Enter. Markdown'is lõigud eraldatakse tühja reaga (üksik reavahetus on ekvivalentne tühikuga). Lõigusisese reavahetuse tegemiseks tuleb rea lõppu jätta vähemalt kaks tühikut. Muud põhilised tekstielemendid on kirjeldatud järgnevas tabelis.

Kood Tulemus
*rõhutatud tekst (kaldkiri)* rõhutatud tekst (kaldkiri)
**tugevalt rõhutatud tekst (rasvane kiri)** tugevalt rõhutatud tekst (rasvane kiri)
`programmikood (masinkiri)` programmikood (masinkiri)
- esimene
- teine
- kolmas
  • esimene
  • teine
  • kolmas
1. esimene
1. teine
1. kolmas
  1. esimene
  2. teine
  3. kolmas
| päis | päis | päis |
| ---- | ---- | ---- |
| sisu | sisu | sisu |
| sisu | sisu | sisu |
päispäispäis
sisusisusisu
sisusisusisu
[Google](http://www.google.com) Google
![](logo.png) No description has been provided for this image

Siin viimases näites vaid lingitakse kõvakettal asuvale pildifailile logo.png (mis peab asuma töölehega samas kataloogis). Pildifaili saab kaasata ka töölehe koosseisu. Selleks tuleb lohistada fail Windows Explorer'ist märkmeraamatu lahtrisse (fail lisandub selle lahtri manuste hulka, vt View | Cell Toolbar | Attachments). Näiteks antud juhul manusfaili nimeks saab attachment:logo.png, mida saab selle lahtri piires kõikjal kasutada.

No description has been provided for this imageKeerulisema kujunduse saavutamiseks tuleb kasutusele võtta juba HTML ja CSS. Üks kujundusvõte, mida Markdown'i lihtsad reeglid ei hõlma, on üla- ja alaindeksite saamine. Näiteks kujundus 1,38 × 10−23 saavutatakse koodiga 1,38&times;10<sup>&minus;23</sup> (erisümbolite saamiseks on kasutatud vastavaid tähekoode. Lisaks on HTML/CSS vajalik jooniste suuruse ja paigutuse kontrollimisel. Näiteks käesoleva lõigu kõrval näha olev pilt on tekitatud koodiga <img src="logo.png" style="width: 40%; float: right; margin-left: 1em;">, mis annab mõista, et pildi laius peab olema 40% lahtri laiusest ja pilt tuleb paigutada vastu parempoolset äärist, nii et lõigu tekst voolab pildist vasakult mööda (jättes veerise, mille laius on võrdne šrifti suurusega).

Seni oleme arvutuslahtrist infot ekraanile trükkinud vaid print-käsuga, mis renderdab toorest ehk kujundamata teksti. Moodulis IPython.display on liides, mis lubab Pythoni koodist Jupyteri töölehele väljastada sellist rikkaliku kujundusega sisu, mida muidu kirjutatakse tekstilahtrisse. Käsuga Markdown või HTML saab luua objekti, mis sisaldab vastavalt Markdown- või HTML-märgendusega teksti, ja käsuga display renderdatakse see ekraanile. Oletagem näiteks, et meil on muutuja x=1.38e-23 ja soovime selle ekraanile kuvada kujul 1,38 × 10−23.

In [15]:
from IPython.display import display, Markdown
from math import floor, log10

x = 1.38e-23
astendaja = floor(log10(x))
mantiss = x / 10**astendaja
md = '**Boltzmanni konstandi** väärtus on %.2f &times; 10<sup>%d</sup> J K<sup>-1</sup>'
md = md % (mantiss, astendaja)
md = md.replace('-', '&minus;').replace('.', ',')
display(Markdown(md))

Boltzmanni konstandi väärtus on 1,38 × 10−23 J K−1

LaTeX¶

Valemid jm matemaatilised konstruktsioonid tuleb kirjeldada LaTeX'i koodis (mida veebilehel renderdab MathJax). Valemi paigutamiseks otse lause sisse tuleb vastav LaTeX'i kood kirjutada dollarimärkide vahele (näiteks $E=mc^2$). Seevastu valemid, mis lähevad omaette reale, kirjutatakse kahekordsete dollarimärkide vahele ($$E=mc^2$$). Piirdume siin vaid mõne näitega, põhjalikum käskude loetelu koos näidetega on eraldi veebilehel.

KoodTulemus
\vec F = m\vec a$$\vec F = m\vec a$$
x = \frac{-b\pm\sqrt{b^2 - 4ac}}{2a}$$x = \frac{-b\pm\sqrt{b^2 - 4ac}}{2a}$$
\sin\left(\alpha+\frac{\pi}{2}\right)=\cos\alpha$$\sin\left(\alpha+\frac{\pi}{2}\right)=\cos\alpha$$
\ln(n!)\approx n\ln n-n$$\ln(n!)\approx n\ln n-n$$
e^x = \sum_{n=0}^\infty \frac{x^n}{n!}$$e^x=\sum_{n=0}^\infty\frac{x^n}{n!}$$

Et genereerida valemeid Pythoni koodist, võib need asetada Markdowni koodi, aga võib kasutada ka eraldi klassi Latex. Tuleb silmas pidada, et LaTeX'i käsud algavad langjoonega, aga viimane võib Pythoni sõne koosseisus tähistada paojada algust. Sel juuhul tuleb kas kasutada kahekordset langkriipsu või lisada sõne ette täht r, mis annab mõista, et paojadasid ei tohi interpreteerida.

In [16]:
from IPython.display import Latex

display(
    Latex('$\vec F = m\vec a$'),
    Latex('$\\vec F = m\\vec a$'),
    Latex(r'$\vec F = m\vec a$')
)
$ ec F = m ec a$
$\vec F = m\vec a$
$\vec F = m\vec a$
Harjutus. Kirjuta koos selgitava tekstiga valem algkiirusega $v_0$ ja kiirendusega $a$ liikuva objekti läbitava teepikkuse leidmiseks!
Harjutus. Kirjuta matemaatilise pendli võnkeperioodi valem!

Andmetüübid¶

Otse programmikoodis sisalduvad andmed, nagu 62 või 1.6e-19 või '%.4g', on literaalid. Andmed võivad pärineda ka failist, internetist, vm. Pythonis on kolm põhilist andmetüüpi (ja vastavat literaali) skalaarsete arvandmete väljendamiseks.

Täisarvud¶

Vaid numbritest koosnev literaal (näiteks 62) esindab täisarvu ehk andmetüüpi int. Andmetüübi saab teada käsuga type:

In [17]:
x = 62
type(x)
Out[17]:
int

Kui kasutada aritmeetikaoperaatoreid +, -, *, // (täisosa) ja % (jääk), siis ka kõik arvutustulemused on täisarvud. Näiteks jagamise pöördoperatsioon täisarvude korral:

In [18]:
jagatav = 62
jagaja = 11

jagatis = jagatav // jagaja
jääk = jagatav % jagaja
print(f'{jagatav} / {jagaja} = {jagatis} (jääk {jääk})')

# taastame algse arvu jagatise ja jäägi kaudu
algne = jagatis * jagaja + jääk
print(f'{jagatis} * {jagaja} + {jääk} = {algne}')
62 / 11 = 5 (jääk 7)
5 * 11 + 7 = 62

Riistvaralise toetusega täisarvud on tüüpiliselt kuni 64-bitised ehk maksimaalse suurusega 263–1=9223372036854775807, kuid tarkvaras võib realiseerida algoritmid kuitahes suurte arvudega opereerimiseks. Pythonis on otsene toetus piiramatu suurusega täisarvudele ja vajadusel läheb süsteem automaatselt üle nende kasutamisele. Näiteks faktoriaali arvutamine ($n! = 1\cdot 2\cdot\ldots\cdot n$) viib väga kiiresti suurte arvudeni:

In [19]:
from math import factorial
factorial(500)
Out[19]:
1220136825991110068701238785423046926253574342803192842192413588385845373153881997605496447502203281863013616477148203584163378722078177200480785205159329285477907571939330603772960859086270429174547882424912726344305670173270769461062802310452644218878789465754777149863494367781037644274033827365397471386477878495438489595537537990423241061271326984327745715546309977202781014561081188373709531016356324432987029563896628911658974769572087926928871281780070265174507768410719624390394322536422605234945850129918571501248706961568141625359056693423813008856249246891564126775654481886506593847951775360894005745238940335798476363944905313062323749066445048824665075946735862074637925184200459369692981022263971952597190945217823331756934581508552332820762820023402626907898342451712006207714640979456116127629145951237229913340169552363850942885592018727433795173014586357570828355780158735432768888680120399882384702151467605445407663535984174430480128938313896881639487469658817504506926365338175055478128640000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
Harjutus. Taanda murd $\frac{2^{19}}{14!}$. Kahe täisarvu suurima ühisteguri leidmiseks saab kasutada funktsiooni math.gcd (greatest common divisor).

Harilike murdudega arvutamiseks on olemas moodul fractions.

Ujukomaarvud¶

Ujukomaarvud on sellised reaalarvud, mille diapasoon on võrdlemisi avar, aga suhteline täpsus (st tüvenumbrite arv) on ligikaudu konstant. Selliseid arve esindab Pythonis andmetüüp float. Literaali interpreteeritakse float-tüüpi väärtusena, kui selles sisaldub kümnendpunkt või see on eksponentkujul, näiteks 3.14 või 3. või 1e6 või 1.6e-19:

In [20]:
x = 3.14
type(x)
Out[20]:
float

Pythoni float hõivab 64 bitti mälu, selle diapasoon on umbes 10−308 kuni 10308 (samuti negatiivsed arvud) ja täpsus 15–17 tüvenumbrit. Viimane tähendab ühtlasi seda, et kahte reaalarvu, mis erinevad rohkem kui umbes 1016 korda, ei ole enam mõtet kokku liita. Või kui kahel reaalarvul 15 või rohkem tüvenumbrit on identsed, siis nende lahutamisel üksteisest saab tulemuse täpsus olema väga väike (vaid mõned tüvenumbrid).

Harjutus. Leia avaldise $$f(x)=2\frac{\sqrt{1+x}-1}{x}$$ väärtus järjest väiksemate $x$ väärtuste korral. Selgub, et mida väiksem on $x$, seda ebatäpsem on tulemus. Piisavalt väikese $x$ korral (umbes 10-16) ei anna arvutus enam üldse mõistlikku tulemust, kuigi teoreetiliselt $$\lim_{x\to 0}f(x)= 1.$$ Selgita selle põhjus ja teisenda avaldis sellisele kujule, mis annab normaalse täpsusega (~15 tüvenumbrit) tulemuse kuitahes väikese $x$ korral.

Reaalarvude piiratud täpsus (ümardusviga) avaldub muuhulgas juhul kui algandmed on kümnendmurrud, sest arvuti mälus on andmed kahendsüsteemis:

In [21]:
0.1 + 0.2 - 0.3
Out[21]:
5.551115123125783e-17

Kui tingimata on tarvis suurt täpsust kümnendmurdudega opereerimisel (nt finantsarvutused), tuleks kasutada püsikomaarve, näiteks moodulist decimal. Meelevaldse täpsusega ujukomaarve pakub moodul mpmath. Selliste spetsiifiliste või suure täpsusega andmetüüpide kasutamisel tuleks float-literaale üldse vältida ja sisestada andmed kas täisarvu või sõne kujul.

Kompleksarvud¶

Kompleksarvudeks nimetatakse matemaatilisi objekte kujul $x+jy$, kus $x$ ja $y$ on reaalarvud ning $j=\sqrt{-1}$. Seega kompleksarvud on reaalarvude laiendus, tehes võimalikuks hulga matemaatilisi operatsioone, mis reaalarvude vallas on võimatud: logaritm või ruutjuur negatiivsest arvust, arkussiinus või -koosinus arvust, mis on suurem kui 1, jms. Kompleksarve esindab Pythonis andmetüüp complex. Vastavalt näiteks kompleksarvu $3+7j$ saab Pythonis kas käsuga complex(3,7) või literaaliga 3+7j. Kompleksarvude korral on aritmeetilised operaatorid mõistagi "teadlikud" sellest, kuidas kompleksarvudega rehkendada:

In [22]:
(9+8j) / (5-2j)
Out[22]:
(1+2j)

Tõepoolest, $$\frac{9+8j}{5-2j}=\frac{9+8j}{5-2j}\times \frac{5+2j}{5+2j}=\frac{29+58j}{29}=1+2j.$$

Mooduli math funktsioonid suudavad opereerida vaid reaalarvude vallas. Näiteks math.sqrt ei suuda arvutada ruutjuurt negatiivsest arvust ega ammugi mitte kompleksarvust (programm katkestatakse veateatega). Kompleksarvude korral tuleb kasutada moodulit cmath.

In [23]:
from cmath import sqrt, exp

x = exp(3-2j)*sqrt(-4+1j)
print(f'{x:.3f}')
34.734-21.376j

Erinevaid arvutüüpe võib matemaatilises avaldises meelevaldselt segada. Kitsam andmetüüp teiseneb vajadusel automaatselt avaramasse formaati, nt täisarv reaalarvuks või reaalarv kompleksarvuks.

Harjutus. Kirjuta programm, mis väljastab ruutvõrrandi $x^2-10x+34=0$ lahendid!

Programmeerimine¶

Selliseid käske nagu kaugus = kiirus * aeg või print(kaugus), ükskõik kui palju neid ritta laduda, täidetakse vääramatult üksteise järel algusest lõpuni. Sellisel viisil ei saa ilmselt realiseerida algoritme, kus tuleb aegajalt võtta vastu otsus, millist käsku järgmisena täita. Samuti on praktiliselt võimatu opereerida andmemassiividega, millest tuleb juttu edaspidi.

Vaatame siin üle lihtsamad struktuurprogrammeerimise kontseptsioonid. Pythoni keelekonstruktsioonid on võrdlemisi lakoonilised ja nende aluseks on põhimõte, et iga (liht)lause kirjutatakse omaette reale, aga lausete grupeerimist (ehk liitlauset) näidatakse kindla taandega (mitte sulgudega, nagu muudes keeltes).

IF-lause¶

Tingimuslause ehk IF-lause on lihtsaim juhtstruktuur, mis võimaldab programmil hargneda, sõltuvalt sellest kas teatud tingimus on täidetud (tõene) või mitte. IF-lause algab võtmesõnaga if, millele järgneb kontrollitav tingimus ja koolon, ning järgmistel ridadel tuleb ühesuguse taandega kirjutada laused, mis kuuluvad täitmisele juhul kui tingimus on tõene. if-blokile võib järgneda meelevaldne hulk elif (else if) blokke ja viimaks veel võtmesõnaga else märgistatud blokk. Kui IF-lause tingimus on väär, siis kontrollitakse esimese elif bloki tingimust, jne. Kui ükski tingimus ei osutu tõeseks, täidetakse laused else-blokis (kui see eksisteerib). Ühe ja sama IF-lause koosseisu kuuluvate if, elif ja else taanded peavad olema muidugi identsed.

Võrdlemisi põhjalik näide IF-lause kohta on ruutvõrrandi $ax^2+bx+c=0$ lahendamine, kus tuleb süstemaatiliselt läbi vaadata kõik variandid. See demonstreerib ka mitut IF-lauset üksteise sees. Katseta erinevate algandmetega:

In [24]:
from math import sqrt

# ruutvõrrandi kordajad
a, b, c = 3, 7, -12

if a == 0:
    if b != 0:
        print( f'x = {-c / b :.3f}' )
    elif c != 0:
        print( 'Võrrandil puuduvad lahendid!' )
    else:
        print( 'Võrrand on rahuldatud iga x korral!' )
else:
    D = b**2 - 4 * a * c  # diskriminant
    if D > 0:
        D = sqrt(D)
        print( f'x1 = {(-b + D) / (2 * a) :.3f}, x2 = {(-b - D) / (2 * a) :.3f}' )
    elif D == 0:
        print( f'x = {-b / (2 * a) :.3f}' )
    else:
        print( 'Võrrandil puuduvad reaalarvulised lahendid!' )
x1 = 1.149, x2 = -3.482

Kui harilik võrdusmärk (=) tähistab muutujale väärtuse omistamist, siis kahekordne võrdusmärk (==) testib vasaku ja parema poole võrdsust. Ülejäänud sarnased operaatorid on < (väiksem), <= (väiksem või võrdne), != (mittevõrdne) jne. Sellisel testil tekib tõeväärtus-tüüpi (bool) tulemus True (tõene) või False (väär). Seega True ja False on tõeväärtusliteraalid.

WHILE-tsükkel¶

Tsükkel on programmi konstruktsioon, mis võimaldab komplekti lauseid korduvalt (tsükliliselt) jooksutada. Seda, kui kaua tsüklit korratakse, sätestatakse harilikult teatud tingimus(t)ega kas tsükli alguses või lõpus, aga vajadusel saab (lõputust) tsüklist ka "jõuga" välja tulla. Lihtsaim on eelkontrolliga ehk WHILE-tsükkel, mille päises, pärast võtmesõna while, kirjeldatakse kohe tingimus, mida kontrollitakse juba enne esmakordset tsüklisse sisenemist. Sellele järgnevad (pärast koolonit) tsükli sisuks olevad laused (kindla taandega, nagu ka IF-lause korral). Pärast nende täitmist läheb täitmisjärg tagasi tsükli algusesse ja kontrollitakse uuesti nimetatud tingimuse täidetust, jne. Kui WHILE-tsükli tingimus on algusest peale väär, siis tsükli keha moodustavaid lauseid ei täideta kordagi.

Eksisteerib ka järelkontrolliga tsükli kontseptsioon (mida läbitakse alati vähemalt korra), kuid Pythonis seda ei ole. Vajadusel saab järelkontrolliga tsükli realiseerida WHILE-tsükli baasil (mille tingimus on püsivalt True), kaasates IF-lause ja käsu break, millega saab tsüklist suvalisel hetkel "välja murda".

Klassikaline matemaatiline näide tsüklite kasutamise kohta on faktoriaali arvutamine:

In [25]:
n = 10

k, fact = n, 1
while k > 1:
    fact *= k
    k -= 1

print( f'arvu {n} faktoriaal on {fact}' )
arvu 10 faktoriaal on 3628800

Eespool oli teatava finantsülesande (konto kasvamine liitintressi ja perioodiliste sissemaksetega) lahend antud üsna keerulise lõppvalemi kujul. Viimase tuletamine nõuab mõningast matemaatilist pingutust. Alternatiivselt võiksime arvutada otse, käies tsükliga kõik perioodid läbi:

In [26]:
algsumma = 1000
säästud = 500
intress = 0.03
aastaid = 7

periood = 0 # mitu perioodi oleme juba läbi käinud
konto = algsumma
while periood < aastaid:
    konto = konto + intress * konto + säästud
    print( f'{periood + 1}. aasta lõpuks on kogunenud {konto:.2f} EUR' )
    periood += 1
1. aasta lõpuks on kogunenud 1530.00 EUR
2. aasta lõpuks on kogunenud 2075.90 EUR
3. aasta lõpuks on kogunenud 2638.18 EUR
4. aasta lõpuks on kogunenud 3217.32 EUR
5. aasta lõpuks on kogunenud 3813.84 EUR
6. aasta lõpuks on kogunenud 4428.26 EUR
7. aasta lõpuks on kogunenud 5061.10 EUR

Selle lähenemise eeliseks on, et tsükli sees teostatav arvutus on märksa lihtsam, paraku arvutamise aeg hakkab sõltuma algandmetest (antud juhul muutujast aastaid).

Veel tõsisema näitena lahendame võrrandi $e^{-x}=x$. See võrrand on transtsendentne, st selle lahendit ei saa avaldada lõpliku arvu elementaarfunktsioonide kombinatsioonina (nagu ruutvõrrandi korral). Seega otsest (kindla arvutustehete arvuga) algoritmi ei eksisteeri, vaid tuleb kasutada mõnda iteratiivset meetodit. Selle käigus, alustades alglähendist, korratakse tsükliliselt teatud arvutusskeemi, nii et tulemus saab seda täpsem, mida kauem arvutada.

Üks lihtsamaid algoritme antud probleemi lahendamiseks on järgmine. Võrrand $e^{-x}=x$ on ekvivalentne võrrandiga $e^{-x}-x=0$, seega algse võrrandi lahendamine tähendab avaldise $e^{-x}-x$ nullkoha leidmist. See avaldis on $x$-i funktsioon, mille me tähistame $f(x)$. Ilmselt me suudame välja pakkuda mingi küllaltki kitsa vahemiku $(a,b)$, kus nullkoht peab asuma, st $a\lt x\lt b$. Kui funktsioon $f(x)$ on monotoonne selles vahemikus, siis ilmselt $f(a)$ ja $f(b)$ on vastasmärgilised, sest nullkoht asub $a$ ja $b$ vahel. Selle vahemiku keskpunkt on $c = (a+b)/2$. Nüüd juhul kui $f(a)$ ja $f(c)$ on vastasmärgilised, siis lahend asub vahemikus $(a,c)$, vastasel korral vahemikus $(c,b)$. Nii saame iga iteratsiooniga lahendi diapasooni 2 korda vähendada, kuni see saab väiksemaks teatud piirväärtusest, mida järgnevas kirjeldab parameeter eps. Trükime välja ka kõik vahepealsed lähendid, et näha kuidas protsess koondub.

In [27]:
from math import exp

f = lambda x: exp(-x) - x
a, b, eps = 0.1, 1.0, 0.0002

while abs(a - b) > eps:
    c = (a + b) / 2
    print( f'{c:.5f}' )
    if f(a) * f(c) < 0:
        b = c
    else:
        a = c
0.55000
0.77500
0.66250
0.60625
0.57812
0.56406
0.57109
0.56758
0.56582
0.56670
0.56714
0.56736
0.56725

Siin avaldise $e^{-x}-x$ arvutamise kood on selle korduvkasutamise hõlbustamiseks kätketud uude funktsiooni nimega f, kasutades lambda-avaldist.

Harjutus. Täienda viimast programmi nii, et lahend leitakse kas täpsusega eps või teostades maksimaalselt max_iter iteratsiooni, kusjuures kõik väljastatavad lähendid on nummerdatud! Kontrolli koodi õigsust, vähendades eps ja/või max_iter väärtust.
Harjutus. Kirjuta programm, mis selgitaks välja, kuidas esitada irratsionaalarv $\pi=3,\!14\ldots$ võimalikult täpselt ratsionaalarvuna $a/b$, kus $a<10000$! Võrdle saadud tulemust ajalooliste arengutega $\pi$ arvutamisel. NB! Kahekordset tsüklit ei lähe vaja (ja see oleks ääretult ebaefektiivne)!

FOR-tsükkel¶

Määratud kordustega ehk FOR-tsükli puhul defineeritakse teatav tsüklimuutuja ehk loendur, mis omandab järgemööda kõik väärtused etteantud vahemikus, ja tsükli keha läbitakse loenduri iga väärtuse puhul (loenduri väärtust saab kasutada tsükli kehas). Näiteks tsükli päis kujul for i in range(m,n), kus m ja n on täisarvud, defineerib loenduri i, mis omandab järgemööda väärtused m, m+1, m+2, …, n-1. Kood range(m,n,s) genereerib arvude järjestuse sammuga s ehk m, m+s, m+2s, ….

Esimene näide trükib tabeli Celsiuse ja Fahrenheiti skaala vastavuse kohta sammuga 10 °C. Teine näide sisaldab koguni kaks FOR-tsüklit üksteise sees ja sel viisil trükib välja korrutustabeli. Koodiga 5d vormindatakse argument täisarvuna, mis on paremjoondatud 5 sümboli laiusel väljal.

In [28]:
print('   °C    °F')
for C in range(0, 101, 10):
    print( f'{C:5d} {9 * C // 5 + 32:5d}' )
   °C    °F
    0    32
   10    50
   20    68
   30    86
   40   104
   50   122
   60   140
   70   158
   80   176
   90   194
  100   212
In [29]:
for x in range(1,11):
    for y in range(1,11):
        print( f'{x * y:5d}', end='')
    print() # reavahetus
    1    2    3    4    5    6    7    8    9   10
    2    4    6    8   10   12   14   16   18   20
    3    6    9   12   15   18   21   24   27   30
    4    8   12   16   20   24   28   32   36   40
    5   10   15   20   25   30   35   40   45   50
    6   12   18   24   30   36   42   48   54   60
    7   14   21   28   35   42   49   56   63   70
    8   16   24   32   40   48   56   64   72   80
    9   18   27   36   45   54   63   72   81   90
   10   20   30   40   50   60   70   80   90  100
Harjutus. Koosta programm, mis trükiks ekraanile mõnede elementaarfunktsioonide väärtuste tabeli järgmisel kujul (sh näidatud täpsuse ja joondusega):
  α     sin(α)     tan(α) 
--------------------------
  0°   0.000000   0.000000
  1°   0.017452   0.017455
  2°   0.034899   0.034921
  3°   0.052336   0.052408
  4°   0.069756   0.069927
  5°   0.087156   0.087489
  6°   0.104528   0.105104
  7°   0.121869   0.122785
  8°   0.139173   0.140541
  9°   0.156434   0.158384
 10°   0.173648   0.176327
 11°   0.190809   0.194380
 12°   0.207912   0.212557
...
 81°   0.987688   6.313752
 82°   0.990268   7.115370
 83°   0.992546   8.144346
 84°   0.994522   9.514364
 85°   0.996195  11.430052
 86°   0.997564  14.300666
 87°   0.998630  19.081137
 88°   0.999391  28.636253
 89°   0.999848  57.289962
Harjutus. Arvu $\pi$ saab meelevaldse täpsusega välja arvutada mitmesuguste valemite järgi, näiteks Viète valem: $$\frac{2}{\pi}=\frac{\sqrt 2}2 \cdot \frac{\sqrt{2+\sqrt 2}}2 \cdot \frac{\sqrt{2+\sqrt{2+\sqrt 2}}}2 \cdots$$ Arvuta ja trüki ekraanile selle valemi järgi järjest täpsemad $\pi$ lähendused (näiteks esimesed 20 varianti). Vihjeks olgu öeldud, et vahetulemusi otstarbekalt kasutades saab selle realiseerida hästi lühikese programmiga.

Alamprogramm¶

Kui algoritm on vormistatud omaette alamprogrammina, saab seda hiljem rakendada mitmetele erinevatele andmetele kogu töölehe ulatuses, ilma et peaks algoritmi koodi kopeerima. Pythonis ja enamikes C-juurtega keeltes kutsutakse alamprogramme üldiselt funktsioonideks ning objekti koosseisus funktsioone nimetatakse meetoditeks. Funktsiooni defineerimiseks kirjutatakse kõigepealt võtmesõna def, millele järgneb defineeritava funktsiooni nimi, argumentide loetelu sulgudes ja koolon. Järgmistel ridadel tulevad kindla taandega funktsiooni keha moodustavad laused. Alamprogrammist saab suvalisel hetkel väljuda return-käsuga. Selle käigus võib, aga ei pea tagastama mingit (arvutus)tulemust (vaikimisi tagastatakse spetsiaalne objekt None).

Näitena vormistame alamprogrammina ruutvõrrandi lahendamise algoritmi, mis annab vastuse nii reaal- kui ka kompleksarvuliste kordajate ja lahendite korral. Seejuures eraldame üksteisest kaks selgelt eristuvat ülesannet: lahendite leidmise ja nende väljatrükkimise. Niisiis tuleb koostada kaks funktsiooni, mis mõlemad võtavad 3 parameetrit (ruutvõrrandi kordajad). Esimene funktsioon arvutab ruutvõrrandi lahendid ja tagastab 2-elemendilise ennikuna. Teine funktsioon esmalt pöördub esimese funktsiooni poole võrrandi lahendite saamiseks ja seejärel trükib lahendid välja.

In [30]:
import math, cmath

def leia_ruutlahend(a, b, c):
    D = b**2 - 4 * a * c
    D = math.sqrt(D) if D >= 0 else cmath.sqrt(D)  # lühivariant tingimuslausest
    return (-b + D) / (2 * a), (-b - D) / (2 * a)

def kuva_ruutlahend(a, b, c):
    x1, x2 = leia_ruutlahend(a, b, c)
    print( f'kordajad: a={a}, b={b}, c={c}' )
    print( f'lahendid: {x1:.3f}, {x2:.3f}' )
    print()

kuva_ruutlahend(2, 3, -14)
kuva_ruutlahend(5, 6, 2)
kordajad: a=2, b=3, c=-14
lahendid: 2.000, -3.500

kordajad: a=5, b=6, c=2
lahendid: -0.600+0.200j, -0.600-0.200j

Siin alamprogrammide sees on vahetulemuste säilitamiseks võetud kasutusele lokaalsed muutujad D ja x. Need muutujad on ajutised ja lähevad kaotsi alamprogrammist väljumisel, vastavatele andmeobjektidele eraldatud mälu vabastatakse. Kui juhtumisi eksisteerivad ka samanimelised globaalsed muutujad, siis need on alamprogrammis nähtamatud. Kui alamprogrammis pöördutakse näiteks muutuja D poole, siis esmalt otsitakse seda lokaalsete muutujate (sh funktsiooni parameetrite) seast ja seejärel globaalsete muutujate hulgast. Globaalsetele muutujatele pole lubatud alamprogrammi sees uut väärtust omistada, v.a. juhul kui globaalne muutuja on spetsiaalselt deklareeritud võtmesõnaga global. Globaalsete muutujate mõjutamine alamprogrammi seest on reeglina taunimisväärne kõrvalefekt — programmi toimimine on palju selgem, kui andmed sisenevad alamprogrammi parameetritena ja (arvutus-)tulemus tagastatakse return-käsuga. Sealhulgas näiteks matemaatilisi funktsioone on sel viisil palju mugavam kombineerida liitfunktsiooniks.

Kindel arv positsioonilisi argumente (nagu viimases näites) on kõige traditsioonilisem viis funktsioonide defineerimisel ja kasutamisel. Pythonis on ootuspäraselt märksa rohkem võimalusi funktsiooni parameetrite kirjeldamiseks. Järgmine alamprogramm nimega tuletis arvutab numbriliselt suvalise etteantud funktsiooni 1. või 2. tuletise: $$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h},\quad f''(x) \approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}.$$ See alamprogramm vajab põhimõtteliselt kuni 4 parameetrit:

  • esmalt muidugi funktsioon $f$, mille tuletist tahetakse leida
  • mitmendat järku tuletis arvutada
  • mis kohal $x$ tuletis arvutada
  • kui suurt sammu $h$ kasutada

Neist esimene on tingimata vajalik, aga ülejäänutele võib anda mõistlikud vaikeväärtused, näiteks 1-järku tuletis kohal $x=0$ sammuga $h=10^{-8}$. Vajadusel saame alamprogrammi väljakutsel kirjeldada parameetrid valikuliselt ja nimeliselt:

In [31]:
from math import *

def tuletis(f, järk=1, kohal=0.0, samm=1e-8):
    if järk == 1:
        return (f(kohal + samm) - f(kohal - samm)) / (2 * samm)
    
    if järk == 2:
        return (f(kohal + samm) - 2 * f(kohal) + f(kohal - samm)) / samm**2
    
    raise NotImplementedError('Sellist järku tuletise algoritm ei ole realiseeritud')

print( tuletis(exp) )
print( tuletis(exp, samm=1e-5) )
print( tuletis(sin, järk=2, kohal=pi/2, samm=1e-4) )
print( tuletis(lambda x: sin(x)*cos(x)) )
0.9999999994736442
1.0000000000121023
-0.999999993922529
1.0

Siin kummaski IF-lauses ei lähe ELSE-blokki tarvis, sest return-käsuga minnakse alamprogrammist kohe välja. Lisaks näeme, et funktsioon f ei pea sugugi olema nimeline (nagu on sin või exp), vaid võib olla ka anonüümne, realiseerituna lambda-avaldisena. Mistahes lihtsaid, üherealisi funktsioone võib def asemel realiseerida lambda-avaldise abil, nagu näidatud juba eespool.

Harjutus. Kodeeri alamprogramm, mis realiseeriks eespool sõnastatud finantsülesande lahenduse (parameetrid oleks seega $K_0$, $r$ ja $S$). Python lubab muuhulgas kasutada *rekursiooni*, st funktsiooni kehaks olev programmikood tohib ka sedasama funktsiooni (st iseennast) välja kutsuda. Sel juhul me võime algoritmi $K(n)$ arvutamiseks sõnastada hoopis nii: kui $n=0$, siis $K(n)=K_0$ (algkapital), vastasel korral $K(n)=(1+r)K(n-1)+S$ (vastavalt intressi definitsioonile).

Jadad¶

Kõik senikasutatud arvandmed on olnud skaalarsed ehk üksikväärtused. Vähegi tõsisemate arvutusülesannete lahendamisel läheb tarvis arvumassiive. Tüüpiliselt läheb massiive vaja, kui on tarvis läbi töötada suur hulk ühetaoliseid andmeid. Sellised andmed tekivad näiteks mitmesuguste füüsikaliste mõõtmiste käigus (voolutugevuse sõltuvus pingest, optilise signaali sõltuvus lainepikkusest, deformatsiooni sõltuvus rakendatud jõust, jne). Andmemassiivina tuleb käsitleda ka mahukaid statistilisi andmeid (aegread jpm). Iga selline andmemassiiv võib sisaldada tuhandeid katsepunkte. Eksperimendi või vaatluste tulemuseks võib olla suur hulk selliseid andmemassiive, mille töötlemine ja analüüs oleks mõistlik automatiseerida.

pinge (V)vool (mA)
0,5041
1,0083
1,50125
2,00166
2,50207
3,00248
3,50291
4,00332

Kuigi Pythonis klassikalist massiivi (array) ei eksisteeri, on olemas mitmeid järjestatud ühemõõtmelisi andmemassiive, mille üldnimetus on jada (sequence). Jada iga elemendi poole saab individuaalselt pöörduda elemendi järjekorranumbri ehk indeksi kaudu, mis näidatakse nurksulgudes jadamuutuja nime taga. Üks juba tuttav jada näide on sõne, mille elementideks on kirjasümbolid. Üldisemad jadad on järjend (list) ja ennik (tuple). Järjendit saab muuta, aga sõne ja ennik, olles korra loodud, säilitavad oma sisu kuni oma eksistentsi lõpuni.

Järjend on äärmiselt paindlik dünaamiline andmestruktuur, mis võib sisaldada mistahes tüüpi elemente (sh teisi jadasid). Meid huvitavad siin konkreetsemalt arvujadad. Oletagem näiteks, et takisti takistuse täpseks määramiseks (ja Ohmi seaduse kontrollimiseks) oleme süstemaatiliselt registreerinud voolutugevuse $I$ läbi takisti erinevatel pingetel $U$ (vt tabel). Vastavate järjendite tekitamiseks tuleb lihtsalt elemendid komadega eraldatult asetada nurksulgude vahele:

In [41]:
voldid = [0.50, 1.00, 1.50, 2.00, 2.50, 3.00, 3.50, 4.00]
milliamprid = [41, 83, 125, 166, 207, 248, 291, 332]

Ennik kirjeldatakse samamoodi, aga nurksulgude asemel on ümarsulud.

Lihtsaim ülesanne on milliamprite teisendamine ampriteks. Selleks on ilmselt vaja tsükliga läbi käia kõik elemendid jadas milliamprid:

In [42]:
amprid = []
for vool in milliamprid:
    amprid.append(0.001 * vool)

print(amprid)
[0.041, 0.083, 0.125, 0.166, 0.20700000000000002, 0.248, 0.291, 0.332]

Nagu näha, Pythonis on FOR-tsükkel seotud jada (või täpsemalt, iteraatori) mõistega: tsüklimuutuja vool omandab siin järgemööda kõik väärtused jadas milliamprid. Etteantud jadast saab uue järjendi tekitada ka veidi elegantsemal viisil kasutades matemaatilises hulgateoorias levinud notatsiooni (list comprehension):

In [34]:
amprid = [0.001 * vool for vool in milliamprid]
print(amprid)
[0.041, 0.083, 0.125, 0.166, 0.20700000000000002, 0.248, 0.291, 0.332]

Et saada võimalikult täpne hinnang takistuse $R$ väärtusele, tuleks läbi saadud katsepunktide sobitada Ohmi seadust esindav sirge $I = \frac{1}{R}U$. Kui eeldada, et põhiline mõõtemääramatus kaasneb sõltuva muutuja (antud juhul voolu $I$) mõõtmisega, siis vähimruutude mõttes parim hinnang takistuse väärtusele on $$R=\frac{\sum_k U_k^2}{\sum_k U_k I_k},$$ kus $k$ on katsepunkte nummerdav indeks.

Summade arvutamiseks on erinevaid võimalusi. $\sum_k U_k^2$ arvutamiseks saaksime kasutada lihtsaimat FOR-tsüklit. Kui aga tahame FOR-tsüklis sünkroonselt liikuda läbi mitme (sama pikkusega) jada (nagu $\sum_k U_k I_k$ arvutamisel), tuleb kasutada kas indekseid või funktsiooni zip. Esimesel juhul läheb vaja funktsiooni range. Nimelt range(n) tekitab iteraatori, mis genereerib $n$-elemendilise jada kõik indeksid $0\ldots n-1$ (Pythonis indekseerimine algab nullist). Tarvis on veel ka funktsiooni len, mis annab elementide arvu jadas ehk jada pikkuse (length).

In [36]:
lugeja, nimetaja = 0, 0
for i in range(len(voldid)):
    lugeja += voldid[i] * voldid[i]
    nimetaja += voldid[i] * amprid[i]

print(f'lugeja = {lugeja:.2f}, nimetaja = {nimetaja:.3f}')
lugeja = 51.00, nimetaja = 4.231

Veidi elegantsema lahenduse saab funktsiooniga zip. See tekitab iteraatori, mis võimaldab sünkroonselt liikuda läbi mitme jada:

In [37]:
lugeja, nimetaja = 0, 0
for pinge, vool in zip(voldid, amprid):
    lugeja += pinge * pinge
    nimetaja += pinge * vool    

print(f'lugeja = {lugeja:.2f}, nimetaja = {nimetaja:.3f}')
lugeja = 51.00, nimetaja = 4.231

Nii nagu range, tekitab ka zip iteraatori, millest vajadusel saab järjendi/enniku käsuga list/tuple:

In [38]:
from pprint import pprint #pretty print
pprint( list(zip(voldid, amprid)) )
[(0.5, 0.041),
 (1.0, 0.083),
 (1.5, 0.125),
 (2.0, 0.166),
 (2.5, 0.20700000000000002),
 (3.0, 0.248),
 (3.5, 0.291),
 (4.0, 0.332)]

Veel lihtsama lahenduse summade arvutamiseks pakub funktsioon sum, mis leiab arvujada (järjendi, enniku) või iteraatori summa. Pane tähele, et siin funktsiooni sum argument ei ole nurksulgudes. Selline konstruktsioon tekitab iteraatori (mitte järjendi).

In [43]:
lugeja = sum(pinge * pinge for pinge in voldid)
nimetaja = sum(pinge * vool for pinge, vool in zip(voldid, amprid))
print(f'lugeja = {lugeja:.2f}, nimetaja = {nimetaja:.3f}')
lugeja = 51.00, nimetaja = 4.231

On mõeldav, et juba lähteandmed on sisestatud tabelitaolise struktuurina. Sel juhul muidugi zip-funktsiooni enam tarvis ei lähe:

In [44]:
andmed = ((0.5, 41),
          (1.0, 83),
          (1.5, 125),
          (2.0, 166),
          (2.5, 207),
          (3.0, 248),
          (3.5, 291),
          (4.0, 332))

lugeja = sum(pinge * pinge for pinge, vool in andmed)
nimetaja = sum(0.001 * pinge * vool for pinge, vool in andmed)
print(f'lugeja = {lugeja:.2f}, nimetaja = {nimetaja:.3f}')
lugeja = 51.00, nimetaja = 4.231

Viimaks saame arvutada takistuse:

In [45]:
takistus = lugeja / nimetaja
print( f'takistus = {takistus:.3f} Ω' )
takistus = 12.054 Ω

Mõned muud põhioperatsioonid jadadega:

  • Kahe jada liitmisel A+B moodustatakse uus jada, mis sisaldab (samas järjestuses) esmalt jada A ja seejärel jada B elemente.
  • Koodidega x in A ja x not in A saab kontrollida, kas element x sisaldub (või ei sisaldu) jadas A.
  • A.reverse() muudab järjendi A elementide järjekorra vastupidiseks.
  • No description has been provided for this imageA[a:b] tagastab alam-jada indeksite vahemikus $a\leq i \lt b$. Seda nimetatakse jada "viilutamiseks" (slicing). Nagu näha, element indeksiga $b$ jääb parajasti välja. Sellist skeemi on hulga mugavam ette kujutada nii, et indeksid a ja b osutavad mitte elementidele endile, vaid nende vahekohtadele, nii et jada alguses on indeks 0 ja lõpus len(A). Indeksid tohivad olla ka negatiivsed, sel juhul loendamine toimub jada lõpust. Indeksi tohib ka ära jätta, kui tahetakse kopeerida teatud pikkusega lõik jada alguses või lõpus. Konstruktsiooniga A[a:b:s] tagastatakse elemendid sammuga s.

Enniku saab sageli kirjeldada ka ilma sulge kasutamata:

In [65]:
x = 1, 2, 3
type(x)
Out[65]:
tuple

Alati seda siiski teha ei saa, näiteks funktsiooni väljakutsel f(1,2,3) ja f((1,2,3)) on täiesti erineva tähendusega (esimesel juhul 3 täisarvulist parameetrit, teisel juhul üks ennik).

Igasuguse (piisavalt lühikese) jada saab lihtsasti "lahti pakkida" üksikelementideks:

In [71]:
a, (b, c) = [1, (2, 3)]
print(f'a = {a}, b = {b}, c = {c}')
a = 1, b = 2, c = 3

Muuhulgas lubab see ühe võrdusmärgiga omistada väärtuse mitmele muutujale:

In [67]:
a, b, c = 1, 2, 3

Sellisel viisil on mugav ja efektiivne ära vahetada kahe muutuja väärtused (mis traditsioonilisemas keeles nõuaks ajutise muutuja kasutamist):

In [68]:
a, b = b, a
print('a =', a)
print('b =', b)
a = 2
b = 1

Jadadega seonduvalt vaatame määramata parameetrite arvuga funktsiooni defineerimist:

In [73]:
def keskmine(*arvud):
    if len(arvud) == 0:
        raise RuntimeError('Funktsioon vajab vähemalt ühe argumendi!')
    return sum(arvud) / len(arvud)

keskmine(2, 4, 8, 3, 5)
Out[73]:
4.4

Siin funktsiooni argument kujul *arvud annab mõista, et funktsiooni väljakutsel Python peaks kõik funktsiooni argumendid kokku koguma ühte järjendisse nimega arvud, mida saaks siis funktsiooni kehas kasutada. Antud funktsioon on ilmselt mõttetu, kui parameetrite arv on 0, sel juhul programm katkestatakse veateatega (luuakse üks Pythoni standardne erind RuntimeError).

Operatsioon toimib ka vastupidi: funktsiooni väljakutsel funk(*jada) jada elemendid edastatakse funktsioonile eraldi argumentidena, st funk(*jada) on ekvivalentne koodiga funk(jada[0], jada[1], ...). Koostame näiteks funktsiooni summa, mille argumendid tohivad olla (läbisegi) nii arvud kui ka jadad. Objekti andmetüüpi saab kontrollida funktsiooniga isinstance. Seda algoritmi on mugav realiseerida rekursiooni abil.

In [74]:
def summa(*objektid):
    if len(objektid) == 0:
        raise RuntimeError('Funktsioon vajab vähemalt ühe argumendi!')
    tulemus = 0.
    for x in objektid:
        if isinstance(x, (int, float)):
            tulemus += x
        elif isinstance(x, (tuple, list)):
            tulemus += summa(*x)
        else:
            raise TypeError('Argumendi tüüp on tundmatu!')
    return tulemus

summa(2, 4, (8, 6), 3, [5, 9, (7, 1)])
Out[74]:
45.0
Harjutus. Eespool sisestati katseandmed tabelina andmed, mille iga element kirjeldas 2-elemendilise ennikuna ühte katsepunkti. Realiseeri vähemalt kolm erinevat lahendust, kuidas sellest tabelist eraldi järjenditena kätte saada pinge- ja vooluväärtuste massiivid, näiteks:
  • FOR-tsüklid ja append-meetod
  • list-comprehension
  • zip-funktsioon
Harjutus. Kodeeri funktsioon summa(*objektid) ilma rekursioonita!

Graafikud¶

Matplotlib on mahukas teek mitmesuguste graafikute kujundamiseks Pythonis. Tegemist on keerulise objekt-orienteeritud süsteemiga, kuid moodulisse matplotlib.pyplot on koondatud komplekt funktsioone, mis realiseerivad Matlab'ile sarnase lihtsa käsustiku. Enamus graafikuid saab tehtud selle käsustiku baasil.

Uus joonis luuakse käsuga figure. Viimasel on ka hulk nimelisi argumente, millest olulisimad on joonise mõõdud tollides (figsize) ja punktide (pikslite) arv tolli kohta (dpi). Joonisele paigutatakse üks või mitu teljestikku käsuga subplot või axes. Näiteks käsuga subplot(231) antakse mõista, et joonisele kavatsetakse paigutada kokku kuni 6 graafikut kahes reas ja kolmes veerus ning luuakse ja aktiveeritakse 1. teljestik (st teljestik 1. rea 1. veerus). Seejärel andmete kandmine aktiivsele teljestikule toimub käsuga plot, millele tuleb anda argumentidena andmemassiivid. Käsu plot korduva väljakutsega saab ühele teljestikule lisada rohkem kui ühe andmeseeria. Lisaparameetritena saab näidata iga andmeseeria jaoks ka nime ja kujunduse. Seeriate eristamiseks üksteisest on põhimõtteliselt hulk võimalusi: ühelt poolt joone värv, stiil ja paksus, teiselt poolt sümboli värv, kuju ja suurus. Lihtsamate juhtude jaoks on olemas spetsiaalsed lühikesed koodid. Näiteks koodiga r- tehakse punast värvi (red) pidevjoon, koodiga bo-- tehakse sinised (blue) mummud, mis on kriipsjoonega ühendatud, jne. Andmeseeria tähis ja nimi kantakse legendile, mis tekib käsuga legend. Telgede nimed/tähised ja teljestiku pealkiri kuvatakse käskudega xlabel, ylabel ja title. Käsk grid lisab ruudustikujooned (grid lines). Viimaks graafik kuvatakse käsuga show. Pärast graafiku kuvamist hakkavad järgmised käsud kujundama juba uut graafikut.

Käsud figure ja subplot (või axes) tohib ka üldse ära jätta, sel juhul teljestik luuakse automaatselt. Seejuures käsu plot korral tekib ristkoordinaadistik, käsu polar korral polaarkoordinaadistik.

Niisiis elementaarne kood lihtsa graafiku kujundamiseks on järgmine:

In [5]:
from matplotlib.pyplot import *
from math import sin, cos, pi

n = 200 # mitu andmepunkti
X = [4 * pi * i / n for i in range(n)] # vahemikus 0...4pi
Y1 = [sin(x) for x in X]
Y2 = [cos(x) for x in X]

figure(figsize=(6,3), dpi=120)
plot(X, Y1, 'r-', label='sin(x)')
plot(X, Y2, 'b--', label='cos(x)')
xlim(0, 4*pi)
title('harmooniline võnkumine')
xlabel('aeg')
ylabel('signaal')
legend(loc='lower left')
grid(color='gray', linestyle=':')
show()
No description has been provided for this image

Polaarteljestikuga graafik tekib analoogiliselt:

In [6]:
figure(figsize=(4,4), dpi=120)
n = 300
θ = [2 * pi * i / n for i in range(n)]
R = [2 - sin(6*t) - 0.5 * cos(30*t) for t in θ]
polar(θ, R, 'b-')
ylim(0, 3.5)
show()
No description has been provided for this image

Mitut graafikut sisaldava joonise konstrueerimiseks tuleb käsku subplot välja kutsuda mitu korda. Vaikimisi tekib ristkoordinaadistik, polaarkoordinaadistik luuakse parameetriga projection='polar' või polar=True.

Olgu meil eksperimentaalselt mõõdetud sõltuvus, mida teoreetiliselt peaks kirjeldama astme-, eksponent- või logaritmfunktsioon. Selle veenvaks demonstreerimiseks tuleks graafik (st selle üks või mõlemad teljed) viia logaritmilisse skaalasse, nii et see sõltuvus muutuks lineaarseks:

In [11]:
from math import exp

# "eksperimentaalsed" andmepunktid
X1 = [1, 2, 3, 4, 5]
Y1 = [2.5, 8, 20, 56, 173]

# teoreetiline sõltvus (pidevjoon) peenema sammuga
X2 = [0.1*i for i in range(55)]
Y2 = [exp(x) for x in X2]

figure(figsize=(7,3), dpi=120)
subplot(121)
plot(X1, Y1, "bo", label="katse")
plot(X2, Y2, "r-", label="teooria")
grid()
legend()
subplot(122)
plot(X1, Y1, "bo")
plot(X2, Y2, "r-")
yscale('log')
grid()
show()
No description has been provided for this image

Automaatne skaalakriipsude (ticks) valik ei pruugi olla alati optimaalne. Näiteks äsja saadud graafikutel võiks skaalakriipsud olla kohati tihedama sammuga. Kõige konkreetsem lahendus on kõikide kriipsude (ja vajadusel ka nende tähiste) asukohad jadana ette anda (käsuga xticks või yticks).

In [10]:
figure(figsize=(7,3), dpi=120)
subplot(121)
plot(X1, Y1, 'bo')
plot(X2, Y2, 'r-')
xticks( tuple(range(6)) )
grid()

subplot(122)
plot(X1, Y1, 'bo')
plot(X2, Y2, 'r-')
yscale('log')
xticks( tuple(range(6)) )
asukohad = 1, 3, 10, 30, 100
yticks( asukohad, asukohad )  # kriipsud ja nende tähised
grid()
show()
No description has been provided for this image

Vajadusel saab käskudega xticks ja yticks skaala kriipsud ja tähised ka üldse ära kaotada (andes argumendina tühjad jadad). Koordinaatteljed tervikuna (koos kõigi tähistega) saab ära peita käsuga axis('off').

Et kasutada kujundatud graafikut artiklis või esitluses, tuleks see salvestada. Lihtsaim variant on kasutada käsku savefig, millele tuleb anda failinimi (kui faili teed ei näidata, tekib fail samas kataloogis, kus asub Jupyteri töölehe fail). Failitüüp valitakse automaatselt failinime laiendi järgi. Olemas on nii vektorgraafika (näiteks PDF) kui ka rastergraafika (PNG) formaadid. Parameeter bbox_inches='tight' eemaldab tühja ala joonise ümbert ning parameetriga dpi saab muuta punktitihedust (vajaliku teravusega rasterkujutise saamiseks). Käsk savefig tuleb anda enne käsku show.

Harjutus. Eespool toodi pinge ja voolutugevuse andmed takisti mõõtmisel. Kanna katsepunktid koos lähendussirgega graafikule, vormista see korrektselt ja salvesta viimaks nii PDF kui ka PNG formaadis (kaks savefig käsku järjest). Ava failid ja võrdle nende sisu!
Harjutus. Lissajous' kõverate parameetrilised võrrandid on järgmised: \begin{align} x &= \cos(at+\delta)\\ y &= \sin(bt) \end{align} Siin $t$ muutub vahemikus $0\ldots 2\pi$. Iga kujundi jaoks on $a$, $b$ ja $\delta$ kindlad konstandid. Konstrueerige järgnev joonis kolme Lissajous' kujundiga, mille parameetrid on näidatud vastava teljestiku kohal. No description has been provided for this image

NumPy¶

Python on dünaamiline programmeerimiskeel, näiteks muutujate tüüp võib meelevaldselt muutuda programmi täitmise käigus. Sellest tulenevalt on Pythoni kood võrdlemisi aeglane, mis avaldub hiljemalt siis kui hakata matemaatilisi operatsioone sooritama suurte arvumassiividega (liikudes FOR- või WHILE-tsükliga element-haaval läbi massiivi). Igasugused standardsed operatsioonid arvumassiividega tuleks võimaluste piires vektoriseerida, st vältida ilmutatud kujul tsükli kasutamist ja realiseerida mõnes kiires madala taseme keeles, nagu Fortran või C. Liiatigi on paljud klassikalised arvutusalgoritmid juba ammu realiseeritud tuntud teekides ja neid ei ole mingit mõtet hakata Pythonis ümber kirjutama.

Pythoni teek NumPy (Numerical Python) realiseerib teadusarvutusteks sobivad andmemassiivid (numpy.ndarray) ja põhilised (vektoriseeritud) matemaatilised operatsioonid nende massiividega. Kui Pythoni järjend või ennik on sisuliselt mäluviitade massiiv (kus viidad võivad osutada suvalistele Pythoni objektidele, mis on arvuti mälus laiali), siis NumPy massiivid on homogeensed ja sidusad, st kõik elemendid on sama tüüpi ja paiknevad arvuti mälus katkematu rivina. Selle tõttu NumPy massiivid nõuavad vähem mälu ja matemaatilised operatsioonid on kiiremad.

Sõltuvalt olukorrast saab NumPy arvumassiivi tekitada mitmesuguste võtetega:

  • anda elemendid ette järjendi või ennikuna (funktsioon array)
  • tekitada lihtsamaid arvumassiive funktsioonidega arange, linspace, geomspace, zeros, jne
  • luua mitmesuguseid juhuslike arvude massiive mooduli numpy.random funktsioonidega
  • lugeda arvumassiiv tekstifailist (loadtxt)

Vaatame esmalt triviaalse näite varal (jällegi takistuse mõõtmise andmed) NumPy massiivi tegemist ja vektoriseeritud arvutusoperatsiooni tähendust. Andmed sisestatakse endiselt Pythoni jadana, mille funktsioon numpy.array konverteerib NumPy massiiviks ehk andmetüübiks numpy.ndarray. Vajadusel oleks võimalik täpsustada ka massiivi elementide andmetüüpi, aga praegu valitakse sisestatud andmete põhjal automaatselt andmetüübiks ujukomaarv (sisuliselt sama mis Pythoni float).

In [75]:
import numpy as np

voldid = np.array( (0.50, 1.00, 1.50, 2.00, 2.50, 3.00, 3.50, 4.00) )
amprid = np.array( (0.041, 0.083, 0.125, 0.166, 0.207, 0.248, 0.291, 0.332) )

voldid / amprid
Out[75]:
array([12.19512195, 12.04819277, 12.        , 12.04819277, 12.07729469,
       12.09677419, 12.02749141, 12.04819277])

Seega õnnestus arvumassiivid voldid ja amprid element-kaupa läbi jagada ilma tsüklit kasutamata.

Nende näidisandmetega jätkates, ka valemis $$R=\frac{\sum_k U_k^2}{\sum_k U_k I_k}$$ sisalduvad arvutusmustrid on hästi tuntud ja triviaalsed realiseerida NumPy vektoriseeritud operatsioonide abiga:

In [76]:
takistus = np.sum(voldid**2) / voldid.dot(amprid)  # või np.sum(voldid**2) / np.sum(voldid*amprid)
print( f'takistus = {takistus:.3f} Ω' )
takistus = 12.054 Ω

Eespool mitmesuguste graafikute tegemisel tuli üsna palju vaeva näha andmete prepareerimisega. NumPy teeb selle töö märksa lihtsamaks. Kasulikuks osutub funktsioon numpy.linspace(a,b,n), mis tekitab $n$-elemendilise arvuvektori, kus arvud vahemikus $a$ kuni $b$ paiknevad konstantse sammuga. Matplotlib aktsepteerib (ja kasutab ka sisemiselt) NumPy massiive. Näiteks siinuse ja koosinuse graafiku saab nüüd nii:

In [34]:
from matplotlib.pyplot import *
from numpy import sin, cos, pi, linspace

x = linspace(0, 4*pi, 200)  # 0 kuni 4*pi, 200 punkti

figure(figsize=(6,3), dpi=120)
plot(x, sin(x), 'r-')
plot(x, cos(x), 'b--')
xlim(0, 4*pi)
xlabel('aeg')
ylabel('signaal')
grid()
show()
No description has been provided for this image

Rõhutame, et siin funktsioonid sin ja cos on võetud moodulist numpy. Vastavad math funktsioonid ei oska NumPy massiividega midagi teha.

Tõsisema näitena joonestame graafiku soojuskiirguse spektritega erinevatel temperatuuridel (Plancki kiirgusseadus): $$B(\lambda, T)=\frac{2\pi c^2}{\lambda^5}\cdot \frac{1}{\exp\left(\frac{hc}{\lambda kT}\right)-1}.$$ Selle võrdlemisi keerulise matemaatilise avaldise saame NumPy abiga loomulikul kujul kirja panna ja tsükleid kasutamata läbi arvutada 200-elemendiliste massiividega. Kuna kvantiteedid muutuvad siin mõlema koordinaadi sihis mitmeid suurusjärke, oleme sunnitud kasutama topelt-logaritmilist skaalat. Käepäraseks osutub siin funktsioon numpy.geomspace(a,b,n), mis tekitab $n$-elemendilise arvuvektori, kus arvud vahemikus $a$ kuni $b$ paiknevad ühtlase sammuga logaritmilisel skaalal.

In [42]:
from numpy import exp, geomspace
from scipy.constants import micro, c, h, k

# λ on lainepikkus või lainepikkuste vektor [µm], T on temperatuur [K]
def planck(λ, T):
    λ = λ * micro
    a = exp(h * c / (λ * k * T)) - 1  # abimuutuja
    return 2 * pi * h * c**2 / λ**5 / a * micro

x = geomspace(0.1, 100, 200) # 0.1 kuni 100 µm, 200 punkti

figure(figsize=(6,4), dpi=120)
for temp, värv in zip( (400, 1000, 3000), ('b', 'g', 'r') ):
    plot(x, planck(x, temp), '-', color=värv, label=f'T={temp:.0f} K')

xlabel('lainepikkus (µm)')
ylabel('kiirguse intensiivsus\nW m$^{-2}$ µm$^{-1}$')
ylim(ymin=1)
xscale('log')
yscale('log')
legend()
grid()
show()
No description has been provided for this image

NumPy massiivide indekseerimine on väga paindlik. Üldjoontes toimib see samuti nagu Pythoni jadade korral, aga on ka hulk lisavõimalusi. Näiteks, kui A on 2D massiiv, siis A[2,5] tagastab 3-nda rea 6-nda veeru elemendi (indekseerimine algab nullist), A[2,:] ehk A[2] tagastab 3-nda rea tervikuna (1D massiivina), A[:,5] tagastab 6-nda veeru, A[:4,:] ehk A[:4] tagastab esimesed 4 rida, A[:,-2:] tagastab viimased 2 veergu, jne. Kui A on 1D massiiv ja I = [2,4,7], siis A[I] tagastab elemendid indeksitega 2, 4 ja 7, A[0::2] tagastab kõik paarituarvulised elemendid ja A[1::2] kõik paarisarvulised elemendid (viimane komponent nurksulgudes määrab sammu). Enamikel juhtudel nende operatsioonidega ei tekitata mitte koopiat massiivi elementidest, vaid tagastatakse vastav vaade (view) algsele massiivile. Samal viisil saab massiivi elemente või terveid alam-massiive üle kirjutada. Toimib ka massiivi "lahtipakkimine", näiteks kujul a,b,c=A, kui massiivi A pikkus esimese mõõtme sihis on 3.

Massiivi mõõtmete teadasaamiseks on hulk võimalusi. Sarnaselt jadadele len(A) tagastab massiivi A pikkuse. Mitmemõõtmelise massiivi korral on tulemuseks massiivi pikkus esimese mõõtme sihis, mida võib tõlgendada kui ridade arvu. Atribuut A.size sisaldab elementide koguarvu massiivis. Kõige detailsemat infot annab A.shape, mis sisaldab ennikuna massiivi kõiki mõõtmeid. NumPy massiivide üks omapära on see, et massiivi "kuju" ei ole kindlalt fikseeritud. Kindel on vaid elementide arv massiivis (mis on kõik arvuti mälus ühe sidusa rivina), aga mitmeks reaks, veeruks või veel kõrgemaks mõõtmeks need elemendid jaotuvad, on juba tõlgendamise küsimus. On lubatud atribuudile shape omistada mingi muu arvukombinatsioon, tingimusel et nende arvude korrutis endiselt võrdub elementide koguarvuga (size) massiivis. Sama teeb meetod reshape. Oletagem näiteks, et eespool kirjeldatud takistuse mõõtmise katseandmed on lihtsuse huvides sisestatud 1D massiivina, kus pinge ja voolu väärtused on vaheldumisi. 2D massiivi saab sellest järgmiselt:

In [86]:
andmed = np.array( (0.50, 41, 1.00, 83, 1.50, 125, 2.00, 166,
                    2.50, 207, 3.00, 248, 3.50, 291, 4.00, 332) )

andmed = andmed.reshape(-1, 2)
print( andmed )
[[   0.5   41. ]
 [   1.    83. ]
 [   1.5  125. ]
 [   2.   166. ]
 [   2.5  207. ]
 [   3.   248. ]
 [   3.5  291. ]
 [   4.   332. ]]

Siin nõutud ridade arvuks on märgitud -1, mis tähendab, et see arvutatakse automaatselt. Paaris- või paarituarvulistest elementidest saaks andmevektorid koostada ka massiivi indekseerimise teel, nagu kirjeldatud eespool.

Harjutus. Reprodutseerige ühes eelnevas ülesandes mainitud Lissajous' kujundite joonis NumPy massiivide abil.

Numbrilised algoritmid¶

Lisaks andmemassiividele sisaldab NumPy ka hulga kasulikke funktsioone, kuid enamus numbrilisi algoritme on koondatud omaette paketti SciPy (Scientific Python). Nimetatud teekidest võib leida vahendeid lineaaralgebra, statistika, interpoleerimise, optimeerimise, diferentsiaalvõrrandite, signaalitöötluse jm probleemide lahendamiseks. Vaatleme siin vaid paari näidet: regressioonsirge sobitamist läbi katsepunktide ja transtsendentse võrrandi lahendamist.

Sirge on 1. järku polünoom, ja suvalist järku polünoomi sobitab vähimruutude mõttes läbi andmepunktide funktsioon numpy.polyfit. See tagastab NumPy massiivina optimaalse polünoomi kordajad. Sirge $y=ax+b$ korral võib kordajad kohe lahti pakkida eraldi muutujatesse (tõus $a$ ja algordinaat $b$), aga üldjuhul saab polünoomi väärtuse arvutada mugavalt funktsiooniga numpy.polyval, mis tahab kõik kordajad ühe massiivina.

Mõõtemääramatuse simuleerimiseks saab juhuslikke arve mooduli numpy.random vahenditega, näiteks randn(n) tekitab standardnormaaljaotusega juhuslike arvude vektori pikkusega n. Juhul kui me tahame saadud juhuslike arvude järjestust (ja seega kogu arvutust) hiljem täpselt reprodutseerida, tuleks eelnevalt ka juhuslike arvude generaator algväärtustada mingi kindla "seemnega", kasutades funktsiooni numpy.random.seed.

In [16]:
from numpy import linspace, polyfit
from numpy.random import randn, seed
seed(42)

n = 20
X = linspace(0, 5, n)
Y = 3 * X + 7 + 2 * randn(n)
a, b = polyfit(X, Y, 1)
Y1 = a * X + b
figure(figsize=(6,4), dpi=120)
plot(X, Y, 'bo')
plot(X, Y1, 'r-')
text(0, 20, f'$y={a:.3f}x+{b:.3f}$', backgroundcolor='white')
grid()
show()
No description has been provided for this image
Harjutus. Simuleerige analoogilisel viisil müraga andmepunktid, mis sirge asemel järgiksid sõltuvust $y=-0.7x^3+5x^2-8x+6$. Veenduge, et polyfit suudab teatud täpsusega tuvastada algse polünoomi kordajad, ja joonestage sile (st peenema sammuga) lähendusfunktsioon samale graafikule.

Eespool me realiseerisime lihtsa algoritmi suvalise mittelineaarse võrrandi (näiteks $e^{-x}=x$) lahendamiseks. Moodulis scipy.optimize on hulk funktsioone, mis teevad sama töö (enamasti keerukama, aga efektiivsema algoritmiga). Igasugused mittelineaarse optimeerimise algoritmid vajavad alglähendit (või vahemikku), millest startida. Kui näiteks võrrandil on mitu lahendit, siis erinevad alglähendid võivad viia erinevate tulemusteni.

In [17]:
from numpy import exp, linspace
from scipy.optimize import brentq
from matplotlib.pyplot import *

f = lambda x: exp(-x) - x
nullkoht, info = brentq(f, 0.1, 1.0, xtol=1e-10, full_output=True)
print( f'nullkoht = {nullkoht:.10f}' )
print( f'funktsiooni väljakutseid: {info.function_calls}' )

x = linspace(0.4, 0.75, 100)
figure(figsize=(3,3), dpi=120)
plot(x, f(x), 'b-')
plot(nullkoht, 0, 'ro')
grid()
show()
nullkoht = 0.5671432904
funktsiooni väljakutseid: 7
No description has been provided for this image

Seega 10 komakoha täpsusega lahendi leidmiseks oli tarvis vaid 7 korda arvutada funktsiooni f(x) väärtust, mis on märksa efektiivsem kui eelmises lahenduses.

Harjutus. Leia vähemalt kaks erinevat lahendit võrrandile $\tan(x)=1/x$. Sobivad alglähendid saab leida graafiliselt.

Failide lugemine-kirjutamine¶

Arvutusteks vajalikud lähteandmed (nt mõne teadusaparaadiga mõõdetud tulemus) asuvad sageli failis pika arvumassiivina. Kõige lihtsam inimloetav failiformaat on selline struktureeritud tekstifail, kus iga kirje on ühes reas ja väljad (st arvutulbad) on eraldatud tühiku, tabulatsioonisümboli või komaga.

Esmalt kasutame elementaarseid Pythoni vahendeid ja loeme sellise faili sisu järjendisse. Käsk open avab faili lugemiseks ('r'), kirjutamiseks ('w') või lisamiseks ('a') ja tagastab file-tüüpi objekti. Viimasel on mitmesuguseid meetodeid teksti lugemiseks failist või kirjutamiseks faili. Näiteks meetodiga readline() loetakse järgmine rida failist (seejuures, faili lõppu jõudmisel tagastatakse tühi sõne). Pärast töö lõppu tuleks fail sulgeda meetodiga close().

Faili lugemisel saab file-tüüpi muutujat käsitleda kui tekstiridade jada. Seega WHILE-tsüklisse kätketud readline() asemel võib kirjutada lihtsalt FOR-tsükli:

In [6]:
tabel = []
f = open('spekter.txt', 'r')
for rida in f:
    sõnad = rida.split()
    arvud = tuple( float(s) for s in sõnad )
    tabel.append( arvud )
f.close()

print(f'failist loeti {len(tabel)} rida')
for rida in tabel[:5]:
    print( '%7.2f %10.6f' % rida )
failist loeti 4592 rida
 380.66   0.000020
 380.73   0.000013
 380.79   0.000009
 380.86   0.000006
 380.92   0.000004

Sõnemeetod split(s) tagastab järjendi kõigist "sõnadest" stringis, kasutades sõnade eraldamiseks stringi s. Vaikimisi on eraldajaiks tühiku või tabulatsioonisümbolid (whitespace). Seejuures ei ole vahet kas sõnad on üksteisest eraldatud ühe või mitme eraldusmärgiga.

Käsuga float(x) antakse mõista, et objekt x (ükskõik mis tüüpi see on) tuleks püüda teisendada ekvivalentseks ujukomaarvuks (st tüüpi float). Kui x on näiteks täisarv, siis operatsioon float(x) on alati edukas. Kui aga x on sõne, siis see peab sisaldama ujukomaarvu tekstina (nt 0.0025 või 2.5e-3), vastasel korral tekib programmiviga. Seega tekstirea konverteerimine arvujadaks ei pruugi tingimata õnnestuda, eriti kui failis esineb vigu, päiseridu, kommentaare, vms. Veakontrolli saab realiseerida TRY…EXCEPT-konstruktsiooniga. Funktsioon sys.stderr.write võimaldab veateate spetsiaalsel kujul välja trükkida.

In [29]:
import sys

tabel = []
with open('spekter2.txt', 'r') as f:
    for rida in f:
        sõnad = rida.split()
        try:
            arvud = tuple( float(s) for s in sõnad )
        except:
            sys.stderr.write( 'Vigases formaadis rida: ' + rida )
            continue
        tabel.append(arvud)

print(f'failist loeti {len(tabel)} rida')
for rida in tabel[:5]:
    print( '%7.2f %10.6f' % rida )
failist loeti 4592 rida
 380.66   0.000020
 380.73   0.000013
 380.79   0.000009
 380.86   0.000006
 380.92   0.000004
Vigases formaadis rida: andmed algavad siit
Vigases formaadis rida: faili l6pp

Kui TRY-plokis tekib viga, siis programmi täitmine sellel kohal katkeb ja luuakse erind (exception), st spetsiaalne objekt, mis sisaldab infot vea kohta. Erind "püütakse kinni" EXCEPT-plokis, kus saab siis otsustada, mida edasi teha (antud juhul trükitakse veateade ja käsuga continue suunatakse täitmisjärg uuesti tsükli algusesse). Kui viga ei teki, siis EXCEPT-plokki üldse ei siseneta.

WITH-lausega tagatakse, et fail suletakse automaatselt, seda isegi juhul kui WITH-plokis peaks tekkima programmiviga. Alternatiivselt võiks kasutada TRY ja FINALLY plokke, kus viimases paikneks f.close().

Andmete kirjutamine faili on täiesti analoogne. Kirjutame faili näiteks raporti eespool simuleeritud elektrimõõtmise andmetega. file-objekti meetod write kirjutab etteantud sõne faili, kuid reavahetuse sümbolit automaatselt ei lisa, selle tõttu sõne ise peab sisaldama reavahetuse sümbolit \n. Faili avamisel saab seadistada ka teksti kodeeringu. Kui on teada, et fail võib sisaldada ASCII kooditabelisse mittekuuluvaid sümboleid (kasvõi täpitähti), tuleks ASCII asemel kasutada mõnda Unicode kodeeringut (näiteks UTF-8, mis on kõige levinum).

In [6]:
andmed =   ((0.5, 41),
            (1.0, 83),
            (1.5, 125),
            (2.0, 166),
            (2.5, 207),
            (3.0, 248),
            (3.5, 291),
            (4.0, 332))

with open('andmed.txt', 'w', encoding='utf8') as f:
    f.write( '%10s%10s%14s\n' % ('pinge (V)', 'vool (A)', 'takistus (Ω)') )
    for rida in andmed:
        pinge, vool = rida
        vool *= 1e-3
        takistus = pinge / vool
        f.write( f'{pinge:10.2f}{vool:10.3f}{takistus:14.2f}\n' )

Faili andmed.txt sisu näeb nüüd välja selline:

 pinge (V)  vool (A)  takistus (Ω)
      0.50     0.041         12.20
      1.00     0.083         12.05
      1.50     0.125         12.00
      2.00     0.166         12.05
      2.50     0.207         12.08
      3.00     0.248         12.10
      3.50     0.291         12.03
      4.00     0.332         12.05
Harjutus. Täienda failist lugemise algoritmi nii, et see toimiks ka juhul kui kümnenderaldaja on koma (mitte punkt). Kasulikuks võib osutuda sõnemeetod replace(vana, uus), mis asendab kõik alamsõned vana sõnega uus.

NumPy ja Pandas loevad/kirjutavad faili sisu ühe käsuga. Lihtsaim on funktsioon numpy.loadtxt, mis vaikimisi eeldab, et andmeveerud on eraldatud tühiku- või tabulatsioonisümbolitega (whitespace). Kümnenderaldaja peab olema punkt.

In [30]:
import numpy as np

tabel = np.loadtxt('spekter.txt')
ridu, veerge = tabel.shape
print(f'failist loeti {ridu} rida')
print( tabel[:5] )
failist loeti 4592 rida
[[3.8066e+02 1.9606e-05]
 [3.8073e+02 1.2948e-05]
 [3.8079e+02 8.8751e-06]
 [3.8086e+02 6.2632e-06]
 [3.8092e+02 4.4803e-06]]

Nagu enne Pythoni järjendis, on ka siin andmed ridadekaupa. Massiivi transponeerimise teel (või lisades funktsioonile loadtxt parameetri unpack=True) saame arvutulbad eraldada, et teostada nendega arvutusi või teha graafik:

In [31]:
from matplotlib.pyplot import *

X, Y = tabel.transpose()

figure(figsize=(7,3), dpi=120)
plot(X, Y, 'b-')
xlabel('lainepikkus')
ylabel('signaal')
grid()
show()
No description has been provided for this image

Tabuleeritud andmete haldamiseks ja analüüsiks on olemas ka spetsiaalne teek Pandas (mis on ehitatud NumPy peale). Näiteks faili sisu kuvamine ilusti kujundatud andmetabelina:

In [55]:
import pandas as pd

tabel = pd.read_table('andmed.txt')
tabel.index += 1  # muidu algavad reanumbrid 0-st
tabel
Out[55]:
pinge (V) vool (mA)
1 0.5 41
2 1.0 83
3 1.5 125
4 2.0 166
5 2.5 207
6 3.0 248
7 3.5 291
8 4.0 332

Tulemuseks on pandas.DataFrame objekt. Veergude nimed võeti seekord automaatselt faili päisest. Kui anname ise tulpadele nimed, siis päiseread tuleks vahele jätta:

In [57]:
tabel = pd.read_table('andmed.txt', skiprows=1, names=('x', 'y'))
tabel.index += 1
tabel
Out[57]:
x y
1 0.5 41
2 1.0 83
3 1.5 125
4 2.0 166
5 2.5 207
6 3.0 248
7 3.5 291
8 4.0 332

Andmed tagastab NumPy massiivi kujul atribuut values.

Funktsiooni read_table korral on andmeveergude eraldaja (parameeter sep) vaikimisi tabulatsioonisümbol (\t), funktsiooni read_csv korral koma. On hulk võimalusi, mida numpy.loadtxt ei paku, näiteks kümnenderaldaja määramine (parameeter decimal), pakitud tekstifaili lugemine (parameeter compression), jms.

Sümbolarvutus¶

Pakett SymPy võimaldab Pythoni keskkonnas teostada sümbolarvutust umbes samal viisil nagu seda saab teha kommertsiaalsetes süsteemides Maple ja Mathematica. Sümbolarvutus opereerib matemaatiliste sümbolitega, mis vaid sümboliseerivad teatud tüüpi arvulisi väärtuseid, kuid ei oma mingit konkreetset arvväärtust. Näiteks samasus $\log ab=\log a+\log b$ kehtib sõltumata sellest, millised on $a$ ja $b$ konkreetsed väärtused, kuigi eeldusega et tegemist on positiivsete reaalarvudega.

SymPy realiseerib hulga spetsiaalseid andmetüüpe matemaatiliste sümbolite, arvude ja operatsioonide kirjeldamiseks. On olemas põhjalik komplekt matemaatilisi funktsioone, mis on kõik "teadlikud" sümbolobjektidest. Selle baasil saab konstrueerida kuitahes keerulisi matemaatilisi avaldisi. SymPy andmetüüpide jaoks on aritmeetilistel operaatoritel spetsiaalne tähendus, nii et avaldised saab ikkagi enam-vähem loomulikul kujul kirja panna. Seejärel saab rakendada mitmesuguseid SymPy käske, näiteks avaldise teisendamiseks sobivale kujule, võrrandi lahendamiseks, summa või integraali arvutamiseks, jne.

Sümbolarvutus on olemuslikult täpne, seega float-tüüpi ligikaudseid arve tuleks vältida. Reeglina piirdutakse täis- või ratsionaalarvudega. Kui on siiski vaja andmeid sisestada või saada tulemusi kümnendmurru kujul, kasutatakse selleks spetsiaalset, piiramata täpsusega andmetüüpi. Sellised andmetüübid vajavad muidugi rohkem arvutusressurssi.

Sümbolid ja avaldised¶

Matemaatilisi sümboleid (muutujaid) esindab andmetüüp sympy.Symbol. Suurema hulga sümbolmuutujaid korraga saab luua käsuga sympy.symbols. Mõlemal juhul tuleb sümbolite nimed sõnena ette anda.

In [61]:
from sympy import *
v,g = symbols('v,g')  # näiteks algkiirus ja raskuskiirendus
type(v)
Out[61]:
sympy.core.symbol.Symbol

Siinjuures Pythoni vastavad muutujanimed ei pruugi üldsegi olla v ja g, aga enamasti pole mingit põhjust miks nad peaksid olema erinevad. Kui matemaatilised sümbolid ja vastavad muutujanimed ühtivad, siis saab kasutada käsku var, mis ise tekitab samanimelised muutujad globaalses nimeruumis.

Sümbolmuutujaid saab nüüd matemaatiliste operatsioonidega kombineerida keerulisemateks avaldisteks. Täisarve esindab SymPy andmetüüp Integer ning ratsionaalarve andmetüüp Rational. Vajadusel tuleb igasugused täpsed arvulised tegurid esitada nende andmetüüpide kaudu.

In [62]:
t = v / g
h = v * t - Rational(1,2) * g * t**2
h
Out[62]:
v**2/(2*g)

Käsuga sympy.srepr saab uurida, millistest komponentidest (SymPy objektidest) avaldis on konstrueeritud:

In [63]:
srepr(h)
Out[63]:
"Mul(Rational(1, 2), Pow(Symbol('g'), Integer(-1)), Pow(Symbol('v'), Integer(2)))"

Nagu siit võib aimata, kõik viis aritmeetilist tehet +, -, *, /, ** taanduvad kolme SymPy klassi Add, Mul ja Pow kombineerimisele.

Tavaliselt kohe pärast sympy importimist antakse käsk init_printing(), mille järel SymPy püüab vormindada väljastatavad tulemused nii ilusti kui võimalik (Jupyteri töölehel MathJax abil):

In [64]:
init_printing()
h
Out[64]:
$$\frac{v^{2}}{2 g}$$

Avaldiste teisendamine¶

Järgnevalt loetleme mõningaid põhioperatsioone sümbolavaldiste teisendamiseks. Sulgavaldiste lahtikorrutamist jm sarnaseid operatsioone teostab sympy.expand, vastupidist operatsiooni viib aga läbi factor:

In [82]:
var('a,b')
expand( (a + b)**2 )
Out[82]:
$$a^{2} + 2 a b + b^{2}$$
In [83]:
factor(_)
Out[83]:
$$\left(a + b\right)^{2}$$

Funktsioonid factor ja expand toimivad muidugi ka murdavaldiste korral. Kui eesmärk ongi ühiste tegurite taandamine lugejast ja nimetajast, on efektiivseim funktsioon cancel:

In [84]:
cancel((a**2-b**2)/(a+b))
Out[84]:
$$a - b$$

Trigonomeetriliste avaldiste puhul tuleb expand asemel kasutada funktsiooni expand_trig ning factor asemel trigsimp:

In [88]:
expand_trig( sin(a + b) )
Out[88]:
$$\sin{\left (a \right )} \cos{\left (b \right )} + \sin{\left (b \right )} \cos{\left (a \right )}$$
In [89]:
trigsimp(_)
Out[89]:
$$\sin{\left (a + b \right )}$$

SymPy avaldised on objektid, millel on hulk meetodeid. Selle asemel, et sympy mooduli funktsioone kasutada, saab mitmeid matemaatilisi operatsioone teostada samanimeliste meetodite väljakutsumise teel:

In [90]:
( (a + b)**2 ).expand()
Out[90]:
$$a^{2} + 2 a b + b^{2}$$

Sellisel viisil on mugav kirjeldada ka pikemat teisenduste jada:

In [91]:
var('x')
( (a + b)**2 ).subs(a,x).subs(b,1).expand()
Out[91]:
$$x^{2} + 2 x + 1$$
Harjutus. Koolifüüsikast on teada lihtne seos ühtlase kiirendusega $a$ liikuva objekti algkiiruse $v_0$, lõppkiiruse $v$ ja läbitud teepikkuse $s$ vahel: $$s=\frac{v^2-v_0^2}{2a}.$$Tuletage see valem, teades et teepikkus avaldub kui keskmise kiiruse ja aja korrutis.

Vaikimisi SymPy ei tee mitte ühtegi "ligikaudu õiget" või "enamasti õiget" teisendust:

In [92]:
var('a,b')
log(a*b).expand()
Out[92]:
$$\log{\left (a b \right )}$$

Siin SymPy keeldub lahti kirjutamast korrutise logaritmi, sest valem $\log ab=\log a + \log b$ ei pruugi kehtida kui $a$ ja $b$ on negatiivsed või koguni komplekssed. Selliseid asjaolusid saab lasta ignoreerida, lisades parameetri force=True:

In [94]:
var('a,b')
log(a*b).expand(force=True)
Out[94]:
$$\log{\left (a \right )} + \log{\left (b \right )}$$

Samas juba sümbolmuutujate defineerimisel saab lisada mitmesuguseid eelduseid arvväärtuse kohta, mida need sümbolid esindavad (vaikimisi ei tehta mingeid eelduseid, st muutuja võib olla ka kompleksarvuline). Näiteks positiivsete reaalarvude korral on mainitud samasus alati kehtiv:

In [95]:
var('a,b', positive=True)
log(a*b).expand()
Out[95]:
$$\log{\left (a \right )} + \log{\left (b \right )}$$

Muud eeldused on real, integer, negative, nonnegative, nonzero jms. Mitte ainult eeldus real, vaid ka positive, negative, jne annavad mõista, et tegemist on reaalarvulise suurusega (kompleksarvude positiivsus või negatiivsus ei omagi tähendust).

Vastupidist laadi operatsioone teostab logcombine:

In [96]:
var('a,b', positive=True)
logcombine(log(a)-log(b))
Out[96]:
$$\log{\left (\frac{a}{b} \right )}$$
In [99]:
var('n', real=True)
logcombine(n*log(a))
Out[99]:
$$\log{\left (a^{n} \right )}$$

Lisades expand-funktsioonile parameetri complex=True, esitatakse kompleksavaldis standardkujul $x+iy$ ($x,y \in \mathbb{R}$). Järgmine näide tuvastab tuntud seose $e^{i\alpha}=\cos\alpha+i\sin\alpha$. Siin SymPy muutuja I esindab imaginaarühikut. Ühtlasi näeme kuidas defineerida kreeka tähestiku sümboleid.

In [100]:
var('alpha', real=True)
exp(I*alpha).expand(complex=True)
Out[100]:
$$i \sin{\left (\alpha \right )} + \cos{\left (\alpha \right )}$$

Kui me ei eelda $\alpha$ kohta midagi, siis saame veidi üldisema vastuse, kus on paratamatult olemas nii $\alpha$ reaal- kui ka imaginaarosa panus:

In [101]:
var('alpha')
exp(I*alpha).expand(complex=True)
Out[101]:
$$i e^{- \Im{\alpha}} \sin{\left (\Re{\left(\alpha\right)} \right )} + e^{- \Im{\alpha}} \cos{\left (\Re{\left(\alpha\right)} \right )}$$

Matemaatilisi konstante $\pi$ ja $e$ esindavad muutujad pi ja E. Veendume, et SymPy suudab tuvastada Euleri samasuse $e^{i\pi}+1=0$:

In [102]:
E**(I * pi) + 1
Out[102]:
$$0$$

Siinjuures väärib rõhutamist, et tulemus 0 ei ole saadud mitte numbrilise arvutamise teel. SymPy arvutused pole oma olemuselt numbrilised ega ligikaudsed. Kasutades oma matemaatilisi "teadmisi", SymPy "sai aru", et avaldise $e^{i\pi}+1$ tõeline, täpne väärtus on täisarv 0. Selline matemaatiline rangus avaldub alati kui avaldistes sisalduvad vaid täpsed numbrilised väärtused või matemaatilised konstandid. Näiteks meelevaldsete arvuliste argumentidega varustatud matemaatilised funktsioonid jäetakse kas üldse välja arvutamata või parimal juhul lihtsustakse niipalju kui võimalik:

In [103]:
from IPython.display import display

display( sqrt(17), sqrt(-8), sin(pi/7), sin(pi/3) )
$$\sqrt{17}$$
$$2 \sqrt{2} i$$
$$\sin{\left (\frac{\pi}{7} \right )}$$
$$\frac{\sqrt{3}}{2}$$

Tulemuseks on harilikult irratsionaalarv, mida ei saakski lõpliku arvu numbritega kirja panna. See muidugi ei takista vajadusel kümnendmurruna arvulise hinnangu saamist kuitahes suure täpsusega. Seda saab kas funktsiooniga N või meetodiga evalf, millele saab ette anda ka vajalike tüvenumbrite arvu:

In [104]:
print( N(sin(pi/7)) )
print( sin(pi/7).evalf(60) )
0.433883739117558
0.433883739117558120475768332848358754609990727787459876444547

Sümbolarvutuse puhul on ka lõpmatus ($\infty$) täiesti aktsepteeritav ja korrektne resultaat:

In [105]:
tan(pi/2)
Out[105]:
$$\tilde{\infty}$$

Sama arvutus numbrilisel kujul annab mingi ebamäärase suure arvu, mis on määratud riistvaraliste ujukomaarvude täpsusega:

In [106]:
import math
math.tan(math.pi/2)
Out[106]:
$$1.633123935319537e+16$$

Eespool mainiti numbrilist arvutusprobleemi, kus andmetüübi float täpsusest jäi vajaka, näiteks avaldise $\frac{\sqrt{1+2x^2}-1}{x}$ väljaarvutamisel:

In [59]:
import math
x = 1.3e-8
(math.sqrt(1 + 2*x**2) - 1)/x
Out[59]:
1.7080354225002406e-08

SymPy-s esindab meelevaldse täpsusega reaalarve andmetüüp Float. float-literaali tuleb mõistagi vältida ja anda kümnendmurd sõnena. Teise parameetrina saab näidata nõutava täpsuse, näiteks 30 tüvenumbrit.

In [65]:
x = Float('1.3e-8', 30)
(sqrt(1 + 2*x**2) - 1)/x
Out[65]:
$$0.0000000129999999999999988559731110623$$

Matemaatiline analüüs¶

Eksisteerib võrdlemisi lihtsaid avaldisi, näiteks $\sin(x)/x$, mida teatud $x$ väärtustel ei ole võimalik otseselt välja arvutada (tekib määramatus $0/0$ vms). Piirväärtuse mõttes vastus siiski eksisteerib. Piirväärtust saab leida funktsiooniga limit:

In [107]:
var('x')
limit(sin(x)/x, x, 0)
Out[107]:
$$1$$

Teatavasti piirväärtus $$\lim_{n\to\infty}\left(1+\frac{1}{n}\right)^n$$ peaks tulema täpselt võrdne naturaallogaritmide alusega. Siin lõpmatuse jaoks kasutame moodulis sympy defineeritud sümbolit oo:

In [108]:
var('n')
limit((1+1/n)**n, n, oo)
Out[108]:
$$e$$

Mistahes piisavalt pideva funktsiooni graafik lokaalselt meenutab sirget. Täpsemal vaatlemisel on kõverus siiski märgatav ja selle kirjeldamiseks tuleks kasutada juba parabooli. Niimoodi jätkates saab funktsiooni kuju lokaalselt kirjeldada kuitahes täpselt võttes piisavalt kõrget järku polünoomi. Seda esitust nimetatakse funktsiooni Taylori reaks. Vastava operatsiooni teostab SymPy funktsioon series, millele saab ka öelda, millise muutuja suhtes, millise punkti ümbruses ja kui kõrget järku liikmeteni tuleks rittaarendus teostada:

In [109]:
var('x')
series(ln(1 + x), x, 0, 2)
Out[109]:
$$x + \mathcal{O}\left(x^{2}\right)$$
In [110]:
var('n')
series((1 + x)**n, x, 0, 2)
Out[110]:
$$1 + n x + \mathcal{O}\left(x^{2}\right)$$
In [111]:
series(cos(x), x, 0, 3)
Out[111]:
$$1 - \frac{x^{2}}{2} + \mathcal{O}\left(x^{3}\right)$$

Nende näidetega oleme reprodutseerinud hulga üldtuntud ligikaudse arvutamise valemeid: $(1+x)^n\approx 1+nx$, $\cos(x)\approx 1-x^2/2$, jne, mis kehtivad tingimusel $x\ll 1$.

Tuletise arvutamine toimub funktsiooniga diff. Vaikimisi võetakse esimest järku tuletis, kuid saab arvutada kõrgemat järku või segatuletist:

In [114]:
var('x,y,n')
diff(x**n, x)
Out[114]:
$$\frac{n x^{n}}{x}$$
In [115]:
simplify(_)
Out[115]:
$$n x^{n - 1}$$
In [116]:
sqrt(x**2 + y**2).diff(x, y)
Out[116]:
$$- \frac{x y}{\left(x^{2} + y^{2}\right)^{\frac{3}{2}}}$$

Käsuga integrate saab arvutada nii määramata kui ka määratud integraale:

In [117]:
var('x')
integrate(1 / x, x)
Out[117]:
$$\log{\left (x \right )}$$
In [118]:
integrate(sqrt(1 - x**2), (x, -1, 1))
Out[118]:
$$\frac{\pi}{2}$$

Viimane integraal,$$\int_{-1}^1 \sqrt{1-x^2} dx,$$ on juhtumisi ühikulise raadiusega poolringi pindala, seega vastus pidigi tulema $\pi/2$.

SymPy on teadlik ka erifunktsioonidest, nagu veafunktsioon, Gammafunktsioon jpt. Sageli mõned integraalid või diferentsiaalvõrrandite lahendid, mis elementaarfunktsioonides ei avaldu, saab väljendada selliste erifunktsioonide kaudu:

In [119]:
display( integrate(exp(-x**2), x) )
display( diff(erf(x), x) )
$$\frac{\sqrt{\pi}}{2} \operatorname{erf}{\left (x \right )}$$
$$\frac{2}{\sqrt{\pi}} e^{- x^{2}}$$

Jupyteri töölehel arvutuse kirjeldamise huvides oleks mõnikord kasulik, kui saaks tervikliku valemina kuvada nii esialgse avaldamata integraali või tuletise kui ka selle väärtuse, mis saadakse sümbolarvutusega. Avaldamata integraal luuakse käsuga Integral ning tuletis käsuga Derivative (need on siis vastavate klasside nimed). Sellisel objektil on meetod doit, mille väljakutsega saab vastava integraali või tuletise välja arvutada. Viimaks sümbolvõrrandi (mis on klassi sympy.Equality objekt) saab tekitada käsuga Eq, mille argumentideks on võrrandi vasak ja parem pool.

In [120]:
var('C') # integreerimiskonstant
a = Integral(exp(-x**2), x)
b = Derivative(erf(x), x)
display( Eq(a, a.doit() + C) )
display( Eq(b, b.doit() + C) )
$$\int e^{- x^{2}}\, dx = C + \frac{\sqrt{\pi}}{2} \operatorname{erf}{\left (x \right )}$$
$$\frac{d}{d x} \operatorname{erf}{\left (x \right )} = C + \frac{2}{\sqrt{\pi}} e^{- x^{2}}$$

Järgmine lõpmatu rea summa arvutamine kirjeldab näiteks järgmist mõttekäiku (vana-Kreeka filosoofiast pärinev "paradoks"): et läbida teatud vahemaa, tuleb esmalt läbida pool sellest ($1/2$), siis järelejäänud teest veel pool ($1/4$), siis järelejäänud teest veel pool ($1/8$), jne. Selleks kulub aega kokku $$\sum_{n=1}^\infty 2^{-n} = \frac{1}{2} + \frac{1}{4} + \frac{1}{8}+ \ldots$$ On ilmne, et selle rea summa peab tulema 1:

In [121]:
var('n')
summation(2**-n, (n, 1, oo))
Out[121]:
$$1$$

Üldise geomeetrilise rea summa $1+q+q^2+\ldots$ pole enam nii kergesti läbinähtav. SymPy nuputab ise välja, et rida koondub vaid $|q|<1$ korral:

In [122]:
var('n,q')
S = Sum(q**n, (n, 0, oo))
Eq(S, S.doit())
Out[122]:
$$\sum_{n=0}^{\infty} q^{n} = \begin{cases} \frac{1}{- q + 1} & \text{for}\: \left|{q}\right| < 1 \\\sum_{n=0}^{\infty} q^{n} & \text{otherwise} \end{cases}$$
Harjutus. Arvuta summa: $$\sum_{k=0}^\infty \frac{(-1)^k}{2k+1}.$$

Võrrandite lahendamine¶

Algebraliste võrrandite (või nende süsteemide) lahendamiseks sobib käsk solve (lineaarsete süsteemide jaoks on ka linsolve). Kui võrrand on antud kujul $f(x)=0$, siis piisab kui käsu solve esimese argumendina anda vaid avaldis $f(x)$, st solve leiab sisuliselt avaldise nullkohad. Kui lahendeid on mitu, siis need tagastatakse järjendina. Ruutvõrrandi lahendamine:

In [123]:
var('a,b,c,x')
solve(a * x**2 + b * x + c, x)
Out[123]:
$$\left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]$$

Diferentsiaalvõrrandite (või nende nende süsteemide) lahendamiseks on analoogne käsk dsolve. Võtame näitena pendli võnkumise, kus taastav jõud (raskus- või elastsusjõud) on võrdeline pendli kõrvalekaldega $x$ tasakaaluasendist. Newtoni II seaduse alusel $m\ddot x=-kx$, millest $\ddot x+\omega^2x=0$, kus $\omega^2=k/m$. Selle lahendamiseks peame kõigepealt käsuga Function defineerima muutuja x, mis esindab koordinaati $x$ kui funktsiooni ajast. Sümboli $\omega$ defineerime kohe reaalarvulise suurusena, et vältida kompleks-kujul lahendit.

In [132]:
x = Function('x') # koordinaat (ajast sõltuv)
t = Symbol('t')   # aeg (sõltumatu muutuja)
w = Symbol('omega', positive=True)
lahend = dsolve(x(t).diff(t,2) + w**2 * x(t), x(t))
lahend
Out[132]:
$$x{\left (t \right )} = C_{1} \sin{\left (\omega t \right )} + C_{2} \cos{\left (\omega t \right )}$$

Et saada konkreetset liikumist kirjeldavat erilahendit, tuleb integreerimiskonstandid $C_1$ ja $C_2$ paika panna algtingimustega. Kuna võrrand oli 2. järku, siis algtingimusi (nagu ka integreerimiskonstante) on 2 tükki, näiteks: alghetkel ($t=0$) pendli kõrvalekalle on maksimaalne (amplituud $A$) ja algkiirus null. Funktsioon dsolve tagastab võrduse (sympy.Equality), mille vasaku ja parema poole saab eraldi kätte atribuutidega lhs (left hand side) ja rhs (right hand side). Kombineeritdes viimast algtingimustega, saab koostada lineaarvõrrandisüsteemi integreerimiskonstantide määramiseks:

In [135]:
var('A,C1,C2')
avaldis = lahend.rhs
konst = solve( (avaldis.subs(t,0) - A, avaldis.diff(t).subs(t,0)), (C1,C2) )
konst
Out[135]:
$$\left \{ C_{1} : 0, \quad C_{2} : A\right \}$$

Viimaks olemegi kätte saanud diferentsiaalvõrrandi erilahendi, mis tuleb ootuspärane:

In [136]:
lahend.subs(konst)
Out[136]:
$$x{\left (t \right )} = A \cos{\left (\omega t \right )}$$

Lineaaralgebra¶

SymPy-s esindab maatrikseid klass Matrix. Maatriksi elemendid võib ette anda Pythoni järjendiga, mille iga element on omakorda järjend, mis kirjeldab ühe rea maatriksis. Alternatiivselt võib kirjeldada maatriksi mõõdud ja anda funktsiooni, mis genereerib elemendid rea- ja veeruindeksite kaudu.

In [125]:
from IPython.display import Math

# abifunktsioon, mida läheb siin ja edaspidi vaja
def kuva(sumbol, objekt):
    display(Math(sumbol + '=' + latex(objekt)))

var('a:d')
A = Matrix([[a, b], [c, d]])
B = Matrix(2, 3, lambda i, j: 3*i+j+1)

kuva('A', A)
kuva('B', B)
$$A=\left[\begin{matrix}a & b\\c & d\end{matrix}\right]$$
$$B=\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\end{matrix}\right]$$

Hästi lihtsaid maatrikseid saab funktsioonidega zeros (kõik elemendid nullid), ones (ühed) ja eye (ühikmaatriks). Indeksite kaudu saab maatriksi elemente nii lugeda kui ka kirjutada. Näiteks veidi kohmakam viis maatriksi B koostamiseks:

In [128]:
B = zeros(2, 3)
for i in range(0,2):
    for j in range(0,3):
        B[i,j] = 3*i+j+1

kuva('B', B)
$$B=\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\end{matrix}\right]$$

Olemas on kõik põhioperatsioonid maatriksitega (liitmine, korrutamine, determinant, transponeerimine, pöördmaatriks, jne):

In [130]:
C = B * B.T
kuva('|A|', A.det())
kuva('C=B\cdot B^T', C)
kuva('C^{-1}', C.inv())
kuva('C\cdot C^{-1}', C*C.inv())
$$|A|=a d - b c$$
$$C=B\cdot B^T=\left[\begin{matrix}14 & 32\\32 & 77\end{matrix}\right]$$
$$C^{-1}=\left[\begin{matrix}\frac{77}{54} & - \frac{16}{27}\\- \frac{16}{27} & \frac{7}{27}\end{matrix}\right]$$
$$C\cdot C^{-1}=\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$

Üheveerulist maatriksit võib vaadelda kui vektorit (õigemini vektori esitust ristbaasis). Selliste vektorite jaoks omavad mõtet skalaar- ja vektorkorrutise (dot/cross product) ning pikkuse mõisted:

In [137]:
var('a1,a2,a3,b1,b2,b3')
a = Matrix([a1,a2,a3])
b = Matrix([b1,b2,b3])
kuva(r'a\cdot b', a.dot(b))
kuva(r'a\times b', a.cross(b))
kuva('|a|', a.norm())
$$a\cdot b=a_{1} b_{1} + a_{2} b_{2} + a_{3} b_{3}$$
$$a\times b=\left[\begin{matrix}a_{2} b_{3} - a_{3} b_{2}\\- a_{1} b_{3} + a_{3} b_{1}\\a_{1} b_{2} - a_{2} b_{1}\end{matrix}\right]$$
$$|a|=\sqrt{\left|{a_{1}}\right|^{2} + \left|{a_{2}}\right|^{2} + \left|{a_{3}}\right|^{2}}$$
Töölehe lõpp
In [10]:
%%html
<style>
div.harjutus { background-color: Khaki; padding: 1ex; }
code { background-color: transparent !important;
      padding: 0px 0.3ex !important; }
table.cmdlist { width: 100%; }
table.center th, table.center td { text-align: center; }
table.left th, table.left td { text-align: left; }
kbd { color: black; background-color: #E0E2EB; border: solid #000000 1px;
    display: inline-block; text-align: center; min-width: 0.9em;
    font-family: "Courier New", Courier, monospace;
    font-size: 90%; font-weight: bold; padding: 0px 2px;
    border-radius: 0.2em; white-space: pre; }
.cm-s-ipython.CodeMirror { background: WhiteSmoke; }
div.untrusted { display: none; }
</style>