# Programovanie

Letná škola FKS 2018

Maťo Gažo, Fero Dráček
(& vykradnuté materiály od Mateja Badina, Feriho Hermana, Kuba, Peťa, Jarných škôl FX a kade-tade po internete)

V tomto kurze si ukážeme základy programovania a naučíme sa programovať matematiku a fyziku.
Takéto vedomosti sú skvelé a budete vďaka nim: 
* vedieť efektívnejšie robiť domáce úlohy
* kvalitnejšie riešiť seminárové a olympiádové príklady
* lepšie rozumieť svetu (IT je dnes na trhu najrýchlejšie rozvíjajúcim sa odvetvím)

Počítač je blbý a treba mu všetko povedať a vysvetliť. Komunikovať sa s ním dá na viacerých úrovniach, my budeme používať Python. Python (názov odvodený z Monty Python's Flying Circus) je všeobecný programovací jazyk, ktorým sa dajú vytvárať webové stránky ako aj robiť seriózne vedecké výpočty. To znamená, že naučiť sa ho nie je na škodu a možno vás raz bude živiť.

Rozhranie, v ktorom píšeme kód, sa volá Jupyter Notebook. Je to prostredie navrhnuté tak, aby sa dalo programovať doslova v prehliadači a aby sa kód dal kúskovať. Pre zbehnutie kúskov programu stačí stlačiť Shift+Enter. 



# Dátové typy a operátory

### Čísla

 podľa očakávaní, vracia trojku


```python
3
```




    3




```python
2+3 # scitanie
```




    5




```python
6-2 # odcitanie
```




    4




```python
10*2 # nasobenie
```




    20




```python
35/5 # delenie
```




    7.0




```python
5//3 # celociselne delenie TODO je toto treba?
```




    1




```python
7%3 # modulo
```




    1




```python
2**3 # umocnovanie
```




    8




```python
4 * (2 + 3) # poradie dodrzane
```




    20



### Logické výrazy


```python
1 == 1 # logicka rovnost
```




    True




```python
2 != 3 # logicka nerovnost
```




    True




```python
1 < 10
```




    True




```python
1 > 10
```




    False




```python
2 <= 2
```




    True



# Premenné

Toto je premenná.

Po stlačení Shift+Enter program v okienku zbehne a premenná sa uloží do pamäte (RAMky, všetko sa deje na RAMke).


```python
a = 2
```

Teraz s ňou možno pracovať ako s bežným číslom.


```python
2 * a
```




    4




```python
a + a
```




    4




```python
a + a*a
```




    6



Možno ju aj umocniť.


```python
a**3
```




    8



Pridajme druhú premennú.


```python
b = 5
```

Nasledovné výpočty dopadnú podľa očakávaní.


```python
a + b
```




    7




```python
a * b
```




    10




```python
b**a
```




    25



Reálne čísla môžeme zobrazovať aj vo vedeckej forme: $2.3\times 10^{-3}$.


```python
d = 2.3e-3
```

### Priklad [0]

FERO DAJ SEM NIECO

# ---------------------------------------------------------------------------

Teraz sa môžeme posunúť k niečomu viac zmysluplnejšiemu. Môžeme začať počítaču zadávať úlohy, ktoré by sme už sami nezvládli !

# Funkcie

Spravme si jednoduchú funkciu, ktorá za nás sčíta dve čísla, aby sme sa s tým už nemuseli trápiť my:


```python
def scitaj(a, b):
    print("Číslo a je {} a číslo b je {}".format(a, b))
    return a + b
```


```python
scitaj(10, 12) # vypise vetu a vrati sucet
```

    Číslo a je 10 a číslo b je 12





    22



Funkcia funguje na celých aj reálnych číslach.

Naša sčítacia funkcia má __štyri podstatné veci__:
1. `def`: toto slovo definuje funkciu.
2. dvojbodka na konci prvého riadku, odtiaľ začína definícia.
3. Odsadenie kódu vnútri funkcie o štyri medzery.
4. Samotný kód. V ňom sa môže diať čokoľvek, Python ho postupne prechádza.
5. `return`: kľúčová vec. Za toto slovo sa píše, čo je output funkcie.

### Úloha 1
Napíšte funkciu `priemer`, ktorá zoberie dve čísla (výšky dvoch chlapcov) a vypočíta ich priemernú výšku.

Ak máš úlohu hotovú, prihlás sa vedúcemu.


```python
# Tvoje riesenie:
def priemer(prvy, druhy):
    return ((prvy+druhy)/2)

priemer(90,20)

```




    55.0



# Poďme na fyziku

V tomto momente môžeme začať používať Python ako sofistikovanejšiu kalkulačku a počítať ňoz základné fyzikálne problémy. 
Predstavme si napríklad, že potrebujeme zistiť, koľko mólov atómov je v dvoch litroch vody.


```python
rho = 1000.0                # hustota
V = 2.0 * 1e-3              # treba premeniť na metre kubické
m = rho * V                 # hmotnosť vody
Mm = (16 + 1 + 1) * 1e-3    # kg/mol

n = m / Mm
print(n)                    # vypiseme
```

    111.1111111111111


A koľko molekúl je v jednom litri vody?


```python
NA = 6.022e23               # Avogadrova konštanta
V = 1e-3
m = rho * V
N = m / Mm * NA             # zamyslite sa nad poradím násobenia a delenia
print(N)
```

    3.345555555555555e+25


Vcelku dosť...

## Úloha 2
Spočítajte objem, ktorý v priemere zaberá jedna molekula ľubovoľnej kvapaliny. Vyjadrite ho v nanometroch kubických.

Urobte to tak, že napíšete funkciu, ktorá bude ako vstup brať:
* objem kvapaliny
* hustotu kvapaliny
* molárnu hmotnosť kvapaliny
a ako výstup to dá objem jednej molekuly kvapaliny. 

Spočítajte to pre metanol, etanol a benzén a výsledky potom ukážte vedúcemu. Neváhajte používať Google.


```python
# Tvoje riesenie:

def ObjemMolekuly(objem, hustota, molarna_hmotnost):
    pocet_molekul = (objem*hustota/molarna_hmostnost)*6.022e23
    return (objem/pocet_molekul)
```

# Zoznamy

Zatiaľ sme sa zoznámili s číslami (celé, reálne), stringami a trochu aj logickými hodnotami.
Zo všetkých týchto prvkov vieme vytvárať množiny, v informatickom jazyku `zoznamy`.

Na úvod sa teda pozrieme, ako s vytvára zoznam (po anglicky `list`). Takúto vec všeobecne nazývame dátová štruktúra.


```python
li = [] # prazdny list
```


```python
v = [4, 2, 3] # list s cislami
```


```python
v
```




    [4, 2, 3]




```python
v[0] # indexovat zaciname nulou!
```




    4




```python
v[1]
```




    2




```python
type(v)
```




    list




```python
w = [5, 'ahoj', True]
```


```python
type(w[0])
```




    int



Čo sa sa stane, ak zoznamy sčítame? Spoja sa.


```python
v + w
```




    [4, 2, 3, 5, 'ahoj', True]



Môžeme ich násobiť?


```python
v * v
```

Smola, nemôžeme. Ale všimnime si, aká užitočná je chybová hláška. Jasne nám hovorí, že nemožno násobiť `list`y.

So zoznamami môžeme robiť rôzne iné užitočné veci. Napríklad ich sčítať.


```python
sum(v)
```




    9



Alebo zistiť dĺžku:


```python
len(v)
```




    3



Alebo ich utriediť:


```python
sorted(v)
```




    [2, 3, 4]



Alebo na koniec pridať nový prvok:


```python
v.append(10)
v
```




    [4, 2, 3, 10]



Alebo odobrať:


```python
v.pop()
v
```




    [4, 2, 3]



### Interval
v zozname sa dá hľadať pomocou intervalov, ktoré sú $\langle x,y)$, čiže uzavretý - otvorený.


```python
li = [2, 5, 7, 8, 10, 11, 14, 18, 20, 25]
len(li)
```




    10




```python
li[1:5] # indexovanie na zaciatku nulou + polouzavrety interval
```




    [5, 7, 8, 10]




```python
li[:3] # prve tri prvky
```




    [2, 5, 7]




```python
li[6:] # prvky zacinajuce na 6tom mieste
```




    [14, 18, 20, 25]




```python
li[2:9:2] # zaciatok:koniec:krok
```




    [7, 10, 14, 20]




```python
del li[2]
li
```




    [2, 5, 8, 10, 11, 14, 18, 20, 25]



Zoznam možno zadefinovať aj cez rozsah:


```python
range(10)
type(range(10))
```




    range




```python
list(range(10))
```




    [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]




```python
list(range(3, 9))
```




    [3, 4, 5, 6, 7, 8]



## Úloha 3
Spočítajte:
* súčet všetkých čísel od 1 do 1000.

Vytvorte zoznam `letnaskola`, ktorý bude obsahovať vašich 5 obľúbených celých čísel. 
* Pridajte na koniec zoznamu číslo 100
* Vymažte druhé číslo zo zoznamu.
* Prepíšte prvé číslo v zozname tak, aby sa rovnalo poslednému v zozname.
* Vypočítajte súčet prvého čísla, posledného čísla a dĺžky zoznamu.


```python
# Tvoje riesenie:
zoznam = list(range(1,1001))
print(sum(zoznam))

letnaskola = [1,1995,12,6,42]
print(letnaskola)
letnaskola.append(100)
print(letnaskola)
del letnaskola[1]
print(letnaskola)
letnaskola[0] = letnaskola[len(letnaskola)-1]
print(letnaskola)

print(letnaskola[0]+letnaskola[len(letnaskola)-1], len(letnaskola))
```

    500500
    [1, 1995, 12, 6, 42]
    [1, 1995, 12, 6, 42, 100]
    [1, 12, 6, 42, 100]
    [100, 12, 6, 42, 100]
    200 5


# For cyklus

Indexy zoznamu môžeme postupne prechádzať. For cyklus je tzv. `iterátor`, ktorý iteruje cez zoznam.


```python
for i in li:
    print(i)
```

    2
    5
    8
    10
    11
    14
    18
    20
    25



```python
for i in li:
    print(i**2)
```

    4
    25
    64
    100
    121
    196
    324
    400
    625


Ako úspešne vytvoriť For cyklus? Podobne, ako pri funkciách:

* `for`: toto slovo je na začiatku.
* `i`: iterovana velicina
* `in`: pred zoznamom, cez ktorý prechádzame (iterujeme).
* dvojbodka na konci prvého riadku.
* kod, ktory sa cykli sa odsadzuje o štyri medzery.

Za pomoci for cyklu môžeme takisto sčítať čísla. Napr. čísla od 0 do 100:


```python
sum = 0

for i in range(101):   # uvedomme si, preco tam je 101 a nie 100
    sum = sum + i      # skratene sum += i
    
print(sum)
```

    5050


## Úloha 4
Spočítajte súčet druhých mocnín všetkých nepárnych čísel od 1 do 100 s využitím for cyklu.


```python
# Tvoje riesenie:
sum = 0
for i in range(101):
    if (i%2 == 1): sum = sum + i**2
print(sum)

```

    166650


# Podmienky

Pochopíme ich na príklade. Zmeňte `a` a zistite, čo to spraví.


```python
a = 5

if a == 3:
    print("cislo a je rovne trom.")
elif a == 5:
    print("cislo a je rovne piatim")
else:
    print("cislo a nie je rovne trom ani piatim.")
```

    cislo a je rovne piatim


Za pomoci podmienky teraz môžeme z for cyklu vypísať napr. len párne čísla. Párne číslo identifikujeme ako také, ktoré po delení dvomi dáva zvyšok nula.

Pre zvyšok po delení sa používa percento:


```python
for i in range(10):
    if i % 2 == 0:
        print(i)
```

    0
    2
    4
    6
    8


Cyklus mozeme zastavit, ak sa porusi nejaka podmienka


```python
for i in range(20):
    print(i)
    if i>10:
        print('Koniec.')
        break
```

    0
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    Koniec.


## Úloha 5
Spočítajte:
* súčet všetkých čísel od 1 do 1000 deliteľné jedenástimi.
* súčet tretích mocnín čísel od 1 do 1000 deliteľných dvanástimi.


```python
# Tvoje riesenie:
sum = 0

for i in range(1001):
    if (i%11 == 0): sum = sum + i

print(sum)

sum = 0

for i in range(1001):
    if (i%12 == 0): sum = sum + i**3
        
print(sum)
```

    45045
    20998994688


## Úloha 6

Teraz, keď už vieme, ako sa zisťuje deliteľnosť, môžeme tiež zistiť, či je zadané číslo prvočíslom. Vymyslite algoritmus, ktorý overí prvočíselnú vlastnosť nejakého čísla.

Výstup by mal byť nasledovný:
```Python
>>> prvocislo(10)
nie
>>> prvocislo(13)
ano
```


```python
# Tvoje riesenie:
from math import *

def prvocislo(cislo):
    vysledok = "ano"
    for i in range(2,int(floor(sqrt(cislo)))+1):
        if (cislo%i == 0): 
            vysledok = "nie"
    return vysledok
            
print(prvocislo(10))
print(prvocislo(13))
```

    nie
    ano


## Úloha 7
Predstavme si, že chceme sčítať nekonečný počet čísel:
$$ \sum_{n=1}^\infty \frac{1}{n^2}.$$
Analytický výsledok takejto sumy je $\pi^2/6$. Koľko členov potrebujeme sčítať, aby presnosť s analytický výsledkov bola určená na tri desatinné miesta?


```python
# Tvoje riesenie:
from math import *
print('Presne ',round(pi**2/6,3))

sum = 0
for i in range(1,2005):
    sum = sum + 1.0/i**2

print(i,' ',sum)
```

    Presne  1.645
    2004   1.6444351893330087


# Konečne poďme na niečo zaujímavé!

V predchádzajúcej časti sme sa zoznámili so základnou syntaxou Pythonu, zoznámili sme sa s Jupyterom a naučili sme sa niečo málo o tom akoby sme mohli sčítať nejaké rady. Nastal však čas pustiť sa do niečo zaujímavejšieho.

V nasledujúcej časti sa pozrieme ako sa dajú pomocou numerickej matematiky a nejakých tých šikovných matematických vzťahov dajú vypočítať konštanty, s ktorými ste sa už stretli.

## Hľadanie hodnoty zlatého rezu $\varphi$

Jednoduché cvičenie na oboznámenie sa s tzv. selfkonzistentným problémom a for cyklom


```python
x = 1;

for i in range (0,20):
    x = 1+1/x
    print (x)
```

## Hľadanie hodnoty Eulerovho čísla $e$

Hoci je to málo známe, Eulerove čislo $e$ sa dá nájsť ako odpoveď na nasledujúcu úlohu (rozdeľuj a ponásob):
Na aké veľké časti treba rozdeliť hocaké číslo tak, aby súčin týchto častí bol maximálny.


```python
import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np
```


```python
Num0 = 25

delitele = np.arange(1, 20 , 1)

casti = []
suciny = []

for i in range (0, 19):
        
    casti.append(Num0/delitele[i])
    
    suciny.append(casti[i]**delitele[i])
```


```python
plt.plot(casti, suciny, "ro-")
plt.show()
```


```python
MinDel = delitele[np.argmax(suciny)-1]
MaxDel = delitele[np.argmax(suciny)+1]
```


```python
delitele =  np.arange(MinDel, MaxDel, (MaxDel - MinDel)/10)

casti = []
suciny = []

for i in range (0, 9):
        
    casti.append(Num0/delitele[i])
    
    suciny.append(casti[i]**delitele[i])
```


```python
plt.plot(casti, suciny, "ro-")
plt.show()
```


```python
casti[np.argmax(suciny)]
```

# Ako sa numericky derivuje, integruje a riešia difky?

## Derivovanie

Na to, aby sme numericky zderivovali nejakú funkciu si musíme najprv uvedomiť, čo tá derivácia vlastne intuitívne robí.
Predstavme si obrázok a v ňom nakreslenú dotyčnicu a malý trojuholníček. Buď sme sa učili alebo sme sa práve teraz dozvedeli, že derivácia funkcie je smernica jej dotyčnice v danom bode. Odtiaľ by už mohlo byť (po nejakých tých obrázkoch) jasné, že to bude nejako takto: 

$$ \frac{df}{dx} = \frac{f(x+h) - f(x-h)}{2h} $$

Zvolíme si malé $h$ a derivujeme.


```python
def func(x):
    return cos(x)

def deriv(f, x, h=0.01):
    return (f(x+h)-f(x-h))/(2*h)

x = np.linspace(0, 2*pi, 101)
#print(x)
y = [func(i) for i in x]
dydx = [deriv(func, i) for i in x]

plt.plot(x, y, label="$f(x)$")
plt.plot(x, dydx, label="$\\frac{df}{dx}$")
plt.xlim([0, 2*pi])
plt.xticks([0, pi, 2*pi])
plt.legend(loc="best")
plt.show()
```

### Išlo by to však aj s lepšou presnosťou?

Odkiaľ spadol magický vzorček vyššie? Tí čo už na podobné akcie chodia dlhšie sa už možno stretli s pojmom Taylorov rozvoj, intuitívne v princípe každá slušná a poslušná funkcia, s ktorou sa stretneme sa dá rozviť v istom okolí bodu $x_0$ do súčtu mocninových funkcií nasledovne:

$$ f(x) \approx f(x_0) + \left.\frac{df}{dx}\right|_{x=x_0}(x-x_0) + \frac{1}{2!}{\left(\left.\frac{df}{dx}\right|_{x=x_0}\right)}^2{(x-x_0)}^2 + \frac{1}{3!}{\left(\left.\frac{df}{dx}\right|_{x=x_0}\right)}^3{(x-x_0)}^3 + ...$$

Ak takto rozvinieme funkcie nielen v okolí bodu $x_0$, ale aj $x_0 + h$ a $x_0 - h$, a následne rozvejo vhodne medzi sebou odčítame, tak sme schopný nájsť hodnotu derivácie funkcie $\frac{df}{dx}$ v bode $x_0$, teda 
$$
\left.\frac{df}{dx}\right|_{x=x_0}
$$

Skúste rozviť funkciu aj v bodoch $x_0+2h$ a $x_0-2h$ a získať vzorček, ktorý nám umožní vypočítať deriváciu s lepšou presnosťou.

## Integrovanie

Dve základné metódy ako sa dá niečo numericky integrovať sú:
* Kvadratúra
* Monte Carlo

Dva spôsoby:
* Kvadratúra (dobrá v 1D, ale so zvyšovaním rozmerov presnosť klesá)
* Monte Carlo (presnosť vždy $O(N^{-1/2})$ (N je počet bodov, ktoré použijeme na výpočet inegr), použiť pri viac ako troch rozmeroch)


Pomôžu nám knižnice.
Teraz si ukážeme kvadratúru.


```python
#from scipy.integrate import quad
import scipy.integrate as sp

def func2(x):
    return exp(-x**2)

sp.quad(func2, -20, 20)   # vysledok je 1/4
```


```python
sqrt(pi)
```

## Hľadanie hodnoty Ludolfovho čísla $\pi$
Pomocou Monte Carlo metódy integrovania sa naučíme ako napríklad vypočítať $\pi$.



```python
import random as rnd

NOP = 50000

CoordXList = [];
CoordYList = [];

for j in range (NOP):
    CoordXList.append(rnd.random())
    CoordYList.append(rnd.random())
```


```python
CircPhi = np.arange(0,np.pi/2,0.01)
```


```python
plt.figure(figsize=(7,7))

plt.plot(
    CoordXList,
    CoordYList,
    color = "red",
    linestyle= "none",
    marker = ","
)
plt.plot(np.cos(CircPhi),np.sin(CircPhi))
#plt.axis([0, 1, 0, 1])
#plt.axes().set_aspect('equal', 'datalim')
plt.show()
```
NumIn = 0

for j in range (NOP):
    
    #if (CoordXList[j] - 0.5)*(CoordXList[j] - 0.5) + (CoordYList[j] - 0.5)*(CoordYList[j] - 0.5) < 0.25:
    if CoordXList[j]*CoordXList[j] + CoordYList[j]*CoordYList[j] <= 1:
        
        NumIn = NumIn + 1;   NumIn/NOP*4
## Diferenciálne rovnice

Príklad neriešiteľnej difky:
$$ y'(x) = \sqrt{1+xy} $$
Príklad Eulerovej (najjednoduchšej) metódy.


```python
N = 101
x = np.linspace(0, 1, N)
dx = x[1] - x[0]
y_eu = np.zeros(N)
y_eu[0] = 1      # pociatocna podmienka

def func(x, y):
    return sqrt(1.0 + x*y)
#    return sin(x)

for i in range(1, N):
    y_eu[i] = y_eu[i-1] + func(x[i-1], y_eu[i-1])*dx
    
plt.plot(x, y_eu)
plt.show()
```

## Difky za pomoci knižníc

Použijeme funkciu `ode` z knižnice `scipy`. Znova riešime
$$y'(x) = \sqrt{1+xy} .$$


```python
from scipy.integrate import odeint

N = 101

def func(y, x):
    return sqrt(1 + x*y)

x = np.linspace(0, 1, N)
y0 = 1.0
y = odeint(func, y0, x) ### Pozri Google Scipy.integrate.odeint

plt.plot(x, (y_eu-y.T[0])/y.T[0])
plt.title("Rozdiel medzi odeint a Eulerom")
plt.show()
```

## Iný príklad na Difky: Exponenciálny rozpad

Ešte sa kukneme na numerické riešenie jednoduchých diferenciálnych rovníc. Metódu demonštrujeme na príklade exponenciálneho rozpadu. Znovu pomocou najjednoduchšej - Eulerovej metódy.

$$\frac{d n}{dt}=-\lambda n$$

$$\frac{n(t+dt)-n(t)}{dt}=-\lambda n(t)$$

$$n(t+dt)=n(t)(1 -\lambda dt), n(0)=N_0$$


```python
N0 = 500

#τ=λdt#
τ = 0.5
```


```python
tauka=[]
NumOfPart=[]

for i in range (15):
    
    if i == 0:
        
        n = N0
    
    else:
        
        n = n*(1 - τ)
        
    tauka.append(i*τ)
    NumOfPart.append(n)
```


```python
plt.plot(tauka, NumOfPart, "ro-")
plt.show()
```

# Obiehanie Zeme okolo Slnka

Fyziku (dúfam!) všetci poznáme.

* gravitačná sila:
$$ \mathbf F(\mathbf r) = -\frac{G m M}{r^3} \mathbf r $$

### Eulerov algoritmus (zlý)
$$\begin{align}
a(t) &= F(t)/m \\
v(t+dt) &= v(t) + a(t) dt \\
x(t+dt) &= x(t) + v(t) dt \\
\end{align}$$

### Verletov algoritmus (dobrý)
$$ x(t+dt) = 2 x(t) - x(t-dt) + a(t) dt^2 $$


```python
from numpy.linalg import norm

G = 6.67e-11
Ms = 2e30
Mz = 6e24
dt = 86400.0
N = int(365*86400.0/dt)
#print(N)

R0 = 1.5e11
r_list = np.zeros((N, 2))
r_list[0] = [R0, 0.0]     # mozno miesat listy s ndarray

v0 = 29.7e3
v_list = np.zeros((N, 2))
v_list[0] = [0.0, v0]

# sila medzi planetami
def force(A, r):
    return -A / norm(r)**3 * r

# Verletova integracia
def verlet_step(r_n, r_nm1, a, dt):  # r_nm1 -- r n minus 1
    return 2*r_n - r_nm1 + a*dt**2

# prvy krok je specialny
a = force(G*Ms, r_list[0])
r_list[1] = r_list[0] + v_list[0]*dt + a*dt**2/2


# riesenie pohybovych rovnic
for i in range(2, N):
    a = force(G*Ms, r_list[i-1])
    r_list[i] = verlet_step(r_list[i-1], r_list[i-2], a, dt)
    
    
plt.plot(r_list[:, 0], r_list[:, 1])
plt.xlim([-2e11, 2e11])
plt.ylim([-2e11, 2e11])
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)
plt.gca().set_aspect('equal', adjustable='box')
#plt.axis("equal")
plt.show()
```

## Pridajme Mesiac


```python
Mm = 7.3e22
R0m = R0 + 384e6
v0m = v0 + 1e3
rm_list = np.zeros((N, 2))
rm_list[0] = [R0m, 0.0]
vm_list = np.zeros((N, 2))
vm_list[0] = [0.0, v0m]

# prvy Verletov krok
am = force(G*Ms, rm_list[0]) + force(G*Mz, rm_list[0] - r_list[0])
rm_list[1] = rm_list[0] + vm_list[0]*dt + am*dt**2/2

# riesenie pohybovych rovnic
for i in range(2, N):
    a = force(G*Ms, r_list[i-1]) - force(G*Mm, rm_list[i-1]-r_list[i-1])
    am = force(G*Ms, rm_list[i-1]) + force(G*Mz, rm_list[i-1]-r_list[i-1])
    r_list[i] = verlet_step(r_list[i-1], r_list[i-2], a, dt)
    rm_list[i] = verlet_step(rm_list[i-1], rm_list[i-2], am, dt)
    
plt.plot(r_list[:, 0], r_list[:, 1])
plt.plot(rm_list[:, 0], rm_list[:, 1])
plt.xlabel("$x$", fontsize=20)
plt.ylabel("$y$", fontsize=20)
plt.gca().set_aspect('equal', adjustable='box')
plt.xlim([-2e11, 2e11])
plt.ylim([-2e11, 2e11])
plt.show()         # mesiac moc nevidno, ale vieme, ze tam je
```

## Úloha pre Vás: Treba pridať Mars :)
Pridajte Mars!

## Matematické kyvadlo s odporom 
Nasimulujte matematické kyvadlo s odporom $\gamma$,
$$ \ddot \theta = -\frac g l \sin\theta -\gamma \theta^2,$$
za pomoci metódy `odeint`.

Alebo pád telesa v odporovom prostredí:
$$ a = -g - kv^2.$$


```python
from scipy.integrate import odeint

def F(y, t, g, k):
    return [y[1], g -k*y[1]**2]

N = 101
k = 1.0
g = 10.0
t = np.linspace(0, 1, N)
y0 = [0.0, 0.0]
y = odeint(F, y0, t, args=(g, k))

plt.plot(t, y[:, 1])
plt.xlabel("$t$", fontsize=20)
plt.ylabel("$v(t)$", fontsize=20)
plt.show()
```

## Harmonický oscilátor pomocou metódy Leapfrog (modifikácia Verletovho algoritmu)


```python
N = 10000
t = linspace(0,100,N)
dt = t[1] - t[0]

# Funkcie
def integrate(F,x0,v0,gamma):
    x = zeros(N)
    v = zeros(N)
    E = zeros(N)    
    
    # Počiatočné podmienky
    x[0] = x0
    v[0] = v0
    
    # Integrovanie rovníc pomocou metódy Leapfrog (wiki)
    fac1 = 1.0 - 0.5*gamma*dt
    fac2 = 1.0/(1.0 + 0.5*gamma*dt)
    
    for i in range(N-1):
         v[i + 1] = fac1*fac2*v[i] - fac2*dt*x[i] + fac2*dt*F[i]
         x[i + 1] = x[i] + dt*v[i + 1]
         E[i] += 0.5*(x[i]**2 + ((v[i] + v[i+1])/2.0)**2)
    
    E[-1] = 0.5*(x[-1]**2 + v[-1]**2)
    
    # Vrátime riešenie
    return x,v,E
```


```python
# Pozrime sa na tri rôzne počiatočné podmienky
F = zeros(N)
x1,v1,E1 = integrate(F,0.0,1.0,0.0) # x0 = 0.0, v0 = 1.0, gamma = 0.0
x2,v2,E2 = integrate(F,0.0,1.0,0.05) # x0 = 0.0, v0 = 1.0, gamma = 0.01
x3,v3,E3 = integrate(F,0.0,1.0,0.4) # x0 = 0.0, v0 = 1.0, gamma = 0.5

# Nakreslime si grafy

plt.rcParams["axes.grid"] = True
plt.rcParams['font.size'] = 14
plt.rcParams['axes.labelsize'] = 18
plt.figure()
plt.subplot(211)
plt.plot(t,x1)
plt.plot(t,x2)
plt.plot(t,x3)
plt.ylabel("x(t)")

plt.subplot(212)
plt.plot(t,E1,label=r"$\gamma = 0.0$")
plt.plot(t,E2,label=r"$\gamma = 0.01$")
plt.plot(t,E3,label=r"$\gamma = 0.5$")
plt.ylim(0,0.55)
plt.ylabel("E(t)")

plt.xlabel("Čas")
plt.legend(loc="center right")

plt.tight_layout()
```

A čo ak bude oscilátor aj tlmenný?


```python
def force(f0,t,w,T):
          return f0*cos(w*t)*exp(-t**2/T**2) 

F1 = zeros(N)
F2 = zeros(N)
F3 = zeros(N)
for i in range(N-1):
    F1[i] = force(1.0,t[i] - 20.0,1.0,10.0)
    F2[i] = force(1.0,t[i] - 20.0,0.9,10.0)
    F3[i] = force(1.0,t[i] - 20.0,0.8,10.0)
```


```python
x1,v1,E1 = integrate(F1,0.0,0.0,0.0)
x2,v2,E2 = integrate(F1,0.0,0.0,0.01)
x3,v3,E3 = integrate(F1,0.0,0.0,0.1)

plt.figure()
plt.subplot(211)
plt.plot(t,x1)
plt.plot(t,x2)
plt.plot(t,x3)
plt.ylabel("x(t)")

plt.subplot(212)
plt.plot(t,E1,label=r"$\gamma = 0$")
plt.plot(t,E2,label=r"$\gamma = 0.01$")
plt.plot(t,E3,label=r"$\gamma = 0.1$")
pt.ylabel("E(t)")

plt.xlabel("Time")
plt.rcParams['legend.fontsize'] = 14.0
plt.legend(loc="upper left")

plt.show()
```

# Lineárna algebra

## Program na dnes

* Matematické operácie s vektormi a maticami. Hľadanie vlastných čísel. 
* Zoznámenia sa s knižnicou `numpy`.

## Vektory


```python
a = np.array([1, 2, 3])
print(a)
print(type(a))

b = np.array([2, 3, 4])
print(a + b)    # spravne scitanie!
```
np.dot(a, b)    # skalarny sucin
a.dot(b)

```python
np.cross(a, b)  # vektorovy sucin
```
np.outer(a, b)  # outer product (jak sa to povie po slovensky?), premyslite si
## Matice
A = np.array([[0, 1], [1, 0]])
print(A)
type(A)AA = np.matrix(A)   # premena na iny datovy typ
type(AA)B = np.array([[1, 0], [0, -1]])
print(B)np.dot(A, B)# pozor!
A*B# Vlastnosti numpy matic
len(A)   # pocet riadkovA.shape   # rozmery# uzitocne vektory/matice, konvencia z Matlabu
N = 3
np.ones(N)      # konstantny vektornp.ones((N, N))    # jednotky
np.eye(N)       # identita
np.zeros((N, N+1))      # nulova matica NxN

```python
# spajanie matic
A = np.ones((3, 3))
B = np.eye(3)
print(A, "\n")
print(B)
np.hstack((A, B))
```

### Prístup k prvkom
# Najprv si pripravime nase vektory
A = np.arange(9)
print(A)
A = A.reshape((3, 3))
A# Pri poliach by sme urobili A[1][1]. Pri vektoroch mozme pouzit aj A[1,1]. 
# Pri maticiach musime pouzit A[1,1], lebo A[1][1] vrati prekvapivy vysledok A[1, 1]A[1][1]A[-1, -1] = 10   # znovu vieme indexovat aj od zadu
AA[0, 0]A[-3, -2]
## Slicing

Ako vyberať jednotlivé stĺpce a riadky matíc?
A = np.arange(9)
A = A.reshape((3, 3))
AA[:, 0]    # Ak nechceme vyberat nejaky usek, pouzijeme dvojbodku.A[1, :]A[[0, 2], :]   # Mozme pouzit aj pole na indexovanie# zmena riadku
A[:, 1] = 4
A# pripocitanie cisla k stlpcu
for i in range(len(A)):
    A[i, 2] += 100

AA = np.arange(25).reshape((5, 5))
AA[2:, 3]
## Generovanie (pseudo)náhodných čísel
# cisla sa menie, pokial nefixneme seed!
np.random.seed(12)

N = 5
a = 2*np.random.rand(N)   # v intervale [0,1]
aA = np.random.rand(N, N)
Aa = np.random.randn(5)
amu, sigma = 1.0, 2.0
a = np.random.randn(1000)*sigma + mu
plt.hist(a, bins=20)
plt.show()a = np.random.randint(10, size=1000)   # cele cisla, pozor na size!
u = plt.hist(a)
#plt.plot(u[1][1:], u[0])
## Vlastné čísla

$$ A v = \lambda v$$
Použijeme funkciu `eig`.

Okrem nej existujú ešte funkcie na vlastné čísla symetrických alebo hermitovských matíc `eigs` a `eigh`.


```python
from numpy.linalg import eig
```
N = 20
A = np.random.rand(N, N)
vals, vecs = eig(A)
print(vals.real)
## Vlastné čísla náhodných matíc

Vygenerujte (100x100) maticu náhodných čísel, získajte jej vlastné čísla a spravte z nich histogram.

N = 100
A = np.random.rand(N, N)
vals, vecs = eig(A)
vals = np.real(vals)
#print(vals)
plt.hist(np.real(vals))
plt.show()
9da35c83665e7c910b8dde4550a4c612a5e68ca1
ipynb
Jupyter Notebook
Programko.ipynb
matoga/LetnaSkolaFKS_notebooks
26faa2d30ee942e18246fe466d9bf42f16cc1433
MIT
Programko.ipynb
matoga/LetnaSkolaFKS_notebooks
26faa2d30ee942e18246fe466d9bf42f16cc1433
MIT
Programko.ipynb
matoga/LetnaSkolaFKS_notebooks
26faa2d30ee942e18246fe466d9bf42f16cc1433
MIT
Qwen/Qwen-72B
1. YES
2. YES