Yritän kertoa lyhyesti, mitä ovat funktiot ja differentiaaliyhtälöt ja miten niitä voi käsitellä pythonilla. Ennen tai jälkeen tähän perehtymisen kannattaa lukea jostain matikan kirjasta matemaattisen asiallinen esittely funktioista.
Selitän ohjelmakoodia melko seikkaperäisesti, mutta tuskin niin perusteellisesti, että tätä juttua voisi kunnolla ymmärtää, ellei tunne ohjelmoinnin perusteita. Oudoille python-komennoille löytyy selitys googlaamalla, joten kaikkea minulta selittämättä jäänyttä ei tarvitse tietää ennakolta. Googlailu on muutenkin olennainen osa ohjelmointia, joten on hyvä harjoitella, millaisilla hauilla parhaiten löytää tarvitsemansa.
On todella vaikeaa ja suuritöistä selittää kirjoittamalla kaikki oleellinen. On paljon helpompaa auttaa pientä ryhmää ratkomaan ongelmia ja selittää asiat juuri heidän tarvitsemallaan tavalla. Isossa luokassa taas ei tiedä, kuka tarvitsee minkäkinlaista selitystä. Oppijan kannalta tylsintä lienee kuunnella opettajan pitkiä yksinpuheluita, joten ehkä sittenkin parasta on kirjoittaa niin perusteelliset ohjeet, että opiskelijat kykenevät työskentelemään osittain omin päin jolloin opettajalle jää aikaa pienten ryhmien opastamiseen. Samalla opettaja oppii, miten parantaa ohjeitaan.
Tämä tiedosto jupyter-notebook -muodossa (Ei lataudu selailemaan, mutta voit ladata sen koneellesi)
Tämä tiedosto python-koodina Voit ladata koneellesi ja suorittaa, jos sinulla on python installoituna.
Valmiit aliohjelmakirjastot¶
Ohjelman aluksi kerron, mitä valmiita aliohjelmakirjastoja haluan käyttää.
pylab-kirjastossa on kuvaajien piirtelyyn sopivia funktioita.import pylab as plt-komennon jälkeen voin kutsua esimerkiksipylab-kirjaston funktiotaplotseuraavasti:plt.plot(x,y)%matplotlib inlinekomento kertoo, että haluan kuvaajat näkyviin tähän samaan dokumenttiin.math-kirjastosta voin kutsua luettelemiani funktiota suoraan; esimerkiksi:y = sin(pi)numpy-kirjastossa on erilaisia numeerisen laskennan funktiota jasympy-kirjastossa symbolisen laskennan funktiota.
import matplotlib.pyplot as plt
# %matplotlib inline
from math import sin, cos, sqrt, pi, exp
import numpy as np
import sympy as sp
# from IPython.display import display
sp.init_printing()
Funktio¶
funktio $f(x) = x*(x-1)*(x+1)$ pythonilla koodattuna
def f(x):
return x * (x - 1) * (x + 1)
# kokeillaan funktiota
# print-käskylle voi antaa monta argumenttia pilkulla erotettuna,
# alla merkkijono ja funktion f antama luku
print("f(1) = ", f(1))
print("\n for-silmukka") # \n tarkoittaa rivinvaihtoa
for i in range(-4, 5): # Huom: range(-4,5) menee -4:stä 4:ään
x = 0.5 * i
print(i, ": f(", x, ") = ", f(x))
f(1) = 0 for-silmukka -4 : f( -2.0 ) = -6.0 -3 : f( -1.5 ) = -1.875 -2 : f( -1.0 ) = 0.0 -1 : f( -0.5 ) = 0.375 0 : f( 0.0 ) = -0.0 1 : f( 0.5 ) = -0.375 2 : f( 1.0 ) = 0.0 3 : f( 1.5 ) = 1.875 4 : f( 2.0 ) = 6.0
Funktioista saa paremman käsityksen piirtelemällä niitä. Funktioiden piirtämiseen käytän pylab-kirjastoa. Esimerkiksi plot-funktio ottaa argumenteiksi listan x- ja y-koordinaatteja.
xx = [i / 10.0 for i in range(-15, 16)]
print("xx = ", xx)
# f(x):n arvot kullakin x:n arvolla, joka kuuluu listaan xx
yy = [f(x) for x in xx]
print()
print("yy = ", yy)
print("\nTulostetaan y neljän desimaalin tarkkuudella")
# Tuntuu, että tähän olisi siistimpikin tapa, mutta en nyt keksi
yys = ["{0:.4g}".format(y) for y in yy]
print("yy = ", yys)
plt.plot(xx, yy)
plt.grid(True)
plt.tick_params(axis='both', colors='green')
plt.show()
xx = [-1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5] yy = [-1.875, -1.3439999999999996, -0.897, -0.5279999999999999, -0.23100000000000026, 0.0, 0.17099999999999996, 0.288, 0.35700000000000004, 0.384, 0.375, 0.33599999999999997, 0.27299999999999996, 0.192, 0.09900000000000002, -0.0, -0.09900000000000002, -0.19200000000000003, -0.273, -0.33599999999999997, -0.375, -0.384, -0.35700000000000004, -0.288, -0.17099999999999996, 0.0, 0.23100000000000023, 0.5279999999999999, 0.8970000000000001, 1.3439999999999996, 1.875] Tulostetaan y neljän desimaalin tarkkuudella yy = ['-1.875', '-1.344', '-0.897', '-0.528', '-0.231', '0', '0.171', '0.288', '0.357', '0.384', '0.375', '0.336', '0.273', '0.192', '0.099', '-0', '-0.099', '-0.192', '-0.273', '-0.336', '-0.375', '-0.384', '-0.357', '-0.288', '-0.171', '0', '0.231', '0.528', '0.897', '1.344', '1.875']
xx = [i/10.0 for i in range(-15,16)] ja yy = [f(x) for x in xx] ovat kätevä pythonmainen tapa luoda lista. Oma tapani on nimetä esimerkiksi useita x:n arvoja sisältävä lista xx:ksi.
pylab on hyvin monipuolinen kirjasto. Käytän siitä vain muutamia perusfunktiota yksinkertaisella tavalla. Käyttämäni komennot olen copy-pastennut googlaamalla löytämistäni esimerkeistä.
Kuvaa voi "täydentää" funktiokutsuilla kuten plt.grid(True), kunnes kutsuu funktiota plt.show(), joka piirtää kuvan.
Differentiaaliyhtälöt¶
Newtonin toinen laki¶
Fysiikan kirjoissa kerrotaan Newtonin laeista. Newtonin toisen lain $F = ma \Leftrightarrow a = F/m$ mukaan kappaleen kiihtyvyys a on sitä suurempi, mitä suurempi voima F siihen kohdistuu ja mitä pienempi on sen massa m.
(Tässä kohdin aloin miettiä, miksi liian kevyttä kiveä ei saa lentämään niin pitkälle kuin sopivan painoista. Ehkä se johtuu siitä, että kiven lisäksi lihasten pitää saada kiihdytettyä käsikin vauhtiin. Tilannetta voisi tarkastella laskemalla, miten tietyn mittainen ja massainen vipuvarsi kiihtyy, kun sitä väännetään tietyllä momentilla. Samalla tulisi tutuksi pyörimisliikkeen ilmiöt. Ei kuitenkaan ryhdytä siihen nyt ettei homma karkaa kokonaan käsistä.)
Kiihtyvyys on nopeuden muutosnopeus eli nopeuden derivaatta ajan suhteen ja ja nopeus on paikan derivaatta ajan suhteen. Newtonin toinen laki on siis differentiaaliyhtälö
$$F = m\frac{d²}{dt²}x(t) = m\frac{d}{dt}v(t) \Leftrightarrow \frac{d}{dt}v(t) = F/m$$
Mörssäri¶
Sovelletaan Newtonin lakia vanhan ajan mörssärin kuulaan, jonka massa on $m$. Mörssärin putki on $\alpha$ asteen kulmassa maahan nähden ja kuulan lähtee sen putkesta nopeudella $v_0$. Unohdamme toistaiseksi ilmanvastuksen, jolloin kuulaan vaikuttaa y-suuntainen voima $-mg$ ja y-suuntainen kiihtyvyys on $a_y(t) = -mg/m$
Paikan ja nopeuden derivaatoista saamme neljän differentiaaliyhtälön ryhmän, jonka kirjoitamme matriisimuotoon:
$$\begin{bmatrix} \frac{d}{dt}v_x(t) \\ \frac{d}{dt}v_y(t) \\ \frac{d}{dt}x(t)\\ \frac{d}{dt}y(t)\\ \end{bmatrix} = \begin{bmatrix} 0\\ -g\\ v_x(t)\\ v_y(t)\\ \end{bmatrix} $$
Kuulan tullessa ulos piipusta $x = 0, y = 0, v_x = v_0*cos(\alpha), v_y = v_0*sin(\alpha)$
Ensimmäinen yhtälö kertoo, että vaakasuora nopeus ei muutu laukaisun jälkeen, koska mikään voima ei vaikuta kuulaan x-suunnassa. y-suuntainen kiihtyvyys on -g, koska kuulaan vaikuttaa voima $-mg$
Numeerinen integrointi¶
Nyt olisi hauska nähdä, miten kuula lentää.
$\frac{d}{dt}x(t) = v(t)\\$ tarkoittaa, että $\frac{x(t+\Delta)-x(t)}{\Delta t} \text{ lähestyy } v(t)\text{:tä} $ kun $\Delta t$ lähestyy nollaa eli on äärettömän pieni. Matemaattisesti tämä esitetään seuraavasti:
$$\lim_{\Delta t \to 0} \frac{x(t+\Delta)-x(t)}{\Delta t} = v(t)$$
$\frac{x(t+\Delta)-x(t)}{\Delta t} = v(t) \Leftrightarrow x(t+\Delta t) = x(t) + v(t)\Delta t$ ei pidä tarkalleen paikkaansa, mutta ei ole ihan huonokaan arvio, jos $\Delta t$ on tarpeeksi lyhyt. Tähän arvioon perustuu differentiaaliyhtälöiden numeerinen ratkaiseminen yksinkertaisimmalla mahdollisella tavalla, Eulerin menetelmällä.
Unohdetaan mörssärin kuula hetkeksi ja otetaan testiesimerkiksi differentiaaliyhtälö $\frac{d}{dt}y(t) = y(t)$. Ylläolevan mukaisesti saamme sille approksimaation: $y(t + \Delta t) = y(t) + y(t) \Delta(t) $
Ratkaistaan ensin yhtälö analyyttisesti ja piirretään kuva helpottamaan approksimaation ymmärtämistä.
t = sp.symbols('t')
y = sp.Function('y')(t)
g = sp.dsolve(sp.diff(y, t) - y)
display(g)
dt = 1.0 # aika-askel
ta = 2.0 # aikavälin alkupiste
tb = ta + dt # aikavälin loppupiste
# aikapisteet, joissa piirretään funktion kuvaaja arvo
tt = [i * 0.1 for i in range(15, 40)]
# funktion arvot ym. aikapisteissä, piirtämistä varten, c_1 = 1
yy = [exp(t) for t in tt]
# funktion arvo aikavälin alkupisteessä
ya = exp(ta) # y = dy/dt = exp(t)
# funktion derivaatta aikavälin alkupisteessä
# Eulerin menetelmällä saatu funktion arvo aikavälin loppupisteessä
dya = ya # y = dy/dt = exp(t)
yb0 = ya + dya * dt
# derivaatan arvo ym. pisteessä
dyb = yb0
# approksimaatio funktion arvolle käyttäen ym. derivaattaa
yb1 = ya + dyb * dt
C1 = yb0 / exp(tb) #
yy_e = [C1 * exp(t) for t in tt]
yb = ya + (dya + dyb) / 2.0
plt.plot(tt, yy, 'k-', label='analyyttinen ratkaisu')
plt.plot([ta, tb], [ya, yb0], 'b-', label='euler')
plt.plot([ta, tb], [ya, yb1], 'g-')
plt.plot([ta, tb], [ya, yb], 'r-', label='pred-corr')
plt.plot(tt, yy_e, 'g:') # apukuvio
plt.plot([tb - 0.5 * dt, tb + 0.5 * dt],
[yb0 - 0.5 * dt * dyb, yb0 + 0.5 * dt * dyb], 'g--')
plt.ylim([0, 35])
plt.grid(True)
plt.legend(loc='upper left')
plt.tick_params(axis='both', colors='green')
plt.show()
Eulerin menetelmällä edetään aikavälin alkupisteestä aika-askeleen verran funktion tangentin suuntaan. Lopputulos olisi tietysti oikea vain, jos funktion derivaatta olisi vakio.
Kuvankin perusteella tuntuisi hyvältä käyttää keskiarvoa funktion derivaattojen arvoista aikavälin päätepisteissä: $\frac{1}{2}(\frac{d}{dt}y(t+\Delta t) + \frac{d}{dt}y(t))$. Harmi kyllä funktion derivaatan arvo $\frac{d}{dt}y(t+\Delta t)$ riippuu funktion arvosta $y(t+\Delta t)$, mitä olemme juuri laskemassa. Eulerin menetelmällä saamme kuitenkin yleensä oikeansuuntaisen arvion y:lle ja sen derivaatalle hetkellä $t+ \Delta t$ joten voimme käyttää sitä approksimaationa derivaatan arvosta. Tätä integroimismenetelmää kutsutaan predictor-corrector menetelmäksi. Kuvaan on piiretty katkoviivalla Eulerin menetelmän mukainen arvio derivaatasta aikavälin päätepisteessä.
Integrointirutiineja¶
Normaalisti käytetään matematiikkakirjastoissa olevia valmiita rutiineja, mutta koodasin tähän uteliaisuudesta ja harjoituksen vuoksi kolme niistä. Ne ovat seuraavan mallin mukaisia, missä integ on menetelmän nimi, joko euler, pred_corr tai rk4.
y(t+dt) = integ(fdydt,yy,uu,dt)
fdydt(yy,uu): käyttäjän kirjoittama funktio, joka palauttaa tilamuuttujien yy derivaatat.yy: järjestelmän tilamuuttujieny(t)lista (=taulukko, vektori, ...)uu = u(t): ulkoiset ohjaukset ja parametritdt: integrointiaskel- funktio
integpalauttaa listanyy1, joka ony(t+dt)
Kirjoittamillani rutiineilla voi siis ratkaista usean muuttujan differentiaaliyhtälöryhmiä.
# enumerate palauttaa vuorollaan kunkin listan yy alkioista y = yy[i]
# ja sitä vastaavan indeksin i
def euler(fdydt, yy, uu, dt):
dy = fdydt(yy, uu)
yy1 = [y + dy[i] * dt for i, y in enumerate(yy)]
return yy1
def pred_corr(fdydt, yy, uu, dt):
dy_a = fdydt(yy, uu)
yy_pred = [y + dy_a[i] * dt for i, y in enumerate(yy)]
dy_b = fdydt(yy_pred, uu)
yy_cor = [y + (dy_a[i] + dy_b[i]) * dt / 2.0 for i, y in enumerate(yy)]
return yy_cor
Allaolevan Runge-Kutta 4 algoritmiin en ole koskaan kunnolla perehtynyt. Runge-Kutta algoritmeja on useita ja ymmärtääkseni ylläoleva predictor-corrector on yksinkertaisin Runge-Kutta algoritmi.
# Runge-Kutta
# algoritmi: http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
def rk4(fdydt, yy, uu, dt):
kk1 = fdydt(yy, uu)
yk1 = [yy[i] + k1 * dt / 2.0 for i, k1 in enumerate(kk1)]
kk2 = fdydt(yk1, uu)
yk2 = [yy[i] + k2 * dt / 2.0 for i, k2 in enumerate(kk2)]
kk3 = fdydt(yk2, uu)
yk3 = [yy[i] + k3 * dt for i, k3 in enumerate(kk3)]
kk4 = fdydt(yk3, uu)
yy1 = [
y + dt / 6.0 * (kk1[i] + 2.0 * kk2[i] + 2.0 * kk3[i] + kk4[i])
for i, y in enumerate(yy)
]
return yy1
Mörssärin kuulan lento¶
Testataan integrointirutiineja mörssärin kuulan lentoon. Halusin nähdä, miten ilmanvastus vaikuttaa kuulan lentoon, joten lisäsin yhtälöihin ilmanvastuksen.
Ilmanvastus on yleensä verrannollinen nopeuden neliöön: $F_{\varrho} = c_{\varrho} v^2$. Koska ilmanvastuksen aiheuttama voima vastustaa kappaleen liikettä sen suunta on nopeudella vastakkainen.
Pallonmuotoiselle kappaleelle löytynee wikipediasta järkevä arvo $c_{\varrho}$:lle, mutta hain sille kokeilemalla arvon, joka vaikutti pikkuisen lentorataan.
Kirjoitin funktion, joka palauttaa ilmanvastuksen aiheuttaman kiihtyvyyden (=jarrutuksen) x- ja y-komponentit. Jätän trigonometrian harjoitukseksi miettiä, miten koodi toimii.
Aluksi asetin kaikki derivaatat nolliksi kuulan osuttua maahan eli kun y < 0, mutta pian innostuin kokeilemaan, miltä näyttäisi, jos kuula ei heti pysähtyisi, vaan vastus kasvaisi kuulan porautuessa syvemmälle maahan.
if (y < 0):
Fv = Fv*(1 - 10*y)
Vastus kasvaa kymmenkertaiseksi kuulan porauduttua metrin syvyyteen. Melko lailla epärealistinen malli kuulan tunkeutumiselle maahan. Vastuksen pitäisi olla ehkä ainakin tuhatkertainen. Numeerisen integroinnin kannalta suuret äkkinäiset muutokset järjestelmän parametreissa ovat kuitenkin hankalia. Ilmalennon voi laskea paljon pidemmällä integrointiaskeleella kuin maahan törmäämisen ja siihen porautumisen. Integrointiaskeleen sopivan pituuden pohdiskelun jätän kuitenkin tekemättä.
def ilmanvastus(xx, roo):
[vx, vy, x, y] = xx
v2 = vx**2 + vy**2
v = sqrt(v2)
Fv = -roo * v2
if (y < 0):
Fv = Fv * (1.0 - 4.0 * y)
a = Fv / m
v = sqrt(v2)
ax = a * vx / v
ay = a * vy / v
return ax, ay
Mörssärin kuula tuskin pomppaa tennispallon lailla uudestaan ilmaan maahan osuttuaan, mutta halusin kuitenkin kokeilla, miten pomppisen voisi mallintaa niin, että se näyttäisi "luonnolliselta". Lisäsin jousivoiman tapaisen voiman, joka työntää kuula ylöspäin sitä voimakkaammin, mitä syvemmällä maassa se on. Mallinnin siis maaperän jonkinlaiseksi kumimatoksi. Kuula olisi alun perin kannattanut ajatella kokoonpainuvaksi osittain kimmoisaksi palloksi niin pomppiminen olisi ollut luontevampi mallintaa.
Täysin kimmoisan törmäyksen voisi mallintaa niin, että kun korkeus on pienempi kuin nolla ja pystysuora nopeus alaspäin, käännetään pystysuora nopeus osoittamaan ylöspäin.
if (y < 0 and vy < 0):
vy = -vy
Halusin kuitenkin törmäyksen täysin kimmoisen ja täysin kimmottoman törmäyksen väliltä. Nyt vastus vie osan liike-energiasta eikä törmäys siksi ole täysin kimmoinen. Kuula kuitenkin kimpoaa maasta, joten törmäys ei ole täysin kimmotonkaan.
def pomppu(xx):
[vx, vy, x, y] = xx
if (y < 0):
ay = -2.5 * g * y
else:
ay = 0
return ay
Tuulesta temmaten arvioin, että kuulaa ylöspäin heittävä voima on viisinkertainen kuulan painoon verrattuna, kun kuula on metrin syvyydessä. Mallinin voima suoraan verrannolliseksi siirtymään, kuten tavallisella jousellakin.
integrointirutiinien vaatima järjestelmän tilamuuttujien aikaderivaatat palauttava funktio
def fdxdt(xx, uu):
[vx, vy, x, y] = xx
[roo] = uu
ax_iv, ay_iv = ilmanvastus(xx, roo)
ay_pmp = pomppu(xx)
dvx = ax_iv
dvy = (ay_iv + ay_pmp - g)
dx = vx
dy = vy
return [dvx, dvy, dx, dy]
g = 9.81 # maan vetovoiman kiihtyvyys
m = 1.0 # kuulan massa
# Mörssärin putken kulma maahan nähden, 45 astetta, joka muunnetaan radiaaneiksi
alpha = 45.0 * (2.0 * pi / 360.0)
# lennon lähtöpisteen x-koordinaatti
x0 = 0.0
# lennon lähtöpisteen y-koordinaatti.
#Pieni positiivinen, koska y0 = 0 on maahan törmäämisen ehto
y0 = 0.1
v0 = 50 # lähtönopeus
vx0 = cos(alpha) * v0
vy0 = sin(alpha) * v0
T0 = 12.0 # Simuloinnnin loppuaika
dt = 0.01 # integrointiaskel
dim = int(T0 / dt) # aikapisteiden lukumäärä
print("aikapisteiden lukumäärä: ", dim)
tt = [i * dt for i in range(0, dim)] # aikapisteet
# Alustetaan taulukko, jossa on tilavektori kullakin ajanhetkellä.
# Riittää ilmoittaa, että siinä on dim-alkiota
xxa = [0 for i in range(0, dim)]
xxa[0] = [vx0, vy0, x0, y0] # tilamuuttujien arvot alkuhetkellä t = 0
xxb = [0 for i in range(0, dim)]
xxb[0] = [vx0, vy0, x0, y0]
xxc = [0 for i in range(0, dim)]
xxc[0] = [vx0, vy0, x0, y0]
roo = 0.02 # ilmavastuksen kerroin
# simuloidaan kolmella eri integrointirutiinilla
for i in range(0, dim - 1):
xxa[i + 1] = euler(fdxdt, xxa[i], [roo], dt)
xxb[i + 1] = pred_corr(fdxdt, xxb[i], [roo], dt)
xxc[i + 1] = rk4(fdxdt, xxc[i], [roo], dt)
# Piirtämistä varten listoista pitää tehdä numpy-taulukoita
zza = np.array(xxa)
zzb = np.array(xxb)
zzc = np.array(xxc)
ax, sxarr = plt.subplots(3, 2, figsize=(10, 8))
plt.suptitle("korkeus paikan funktiona", color='green')
sxarr[0,0].plot(zza[:, 2], zza[:, 3], 'r-', label='euler')
sxarr[1,0].plot(zzb[:, 2], zzb[:, 3], 'b-', label='pred-corr')
sxarr[2,0].plot(zzc[:, 2], zzc[:, 3], 'g-', label='rk4')
# "osasuurennos" edellisestä
sxarr[0,1].plot(zza[:, 2], zza[:, 3], 'r-', label='euler')
sxarr[1,1].plot(zzb[:, 2], zzb[:, 3], 'b-', label='pred-corr')
sxarr[2,1].plot(zzc[:, 2], zzc[:, 3], 'g-', label='rk4')
for i in range(3):
for j in range(2):
sxarr[i,j].tick_params(axis='both', colors='green')
sxarr[i,j].legend(loc='upper left')
sxarr[i,j].grid(True)
sxarr[i,j].axis('equal')
for i in range(3):
sxarr[i, 1].set(xlim=(60, 80), ylim=(-5, 8))
plt.show()
aikapisteiden lukumäärä: 1200
Lentoradan huippukohta on 40 metrin kohdalla. Ellei ilmanvastusta olisi otettu huomioon, kuula lentäisi yli 80 metriä ja lentorata olisi symmetrinen — parabeli
Aikapisteitä on 1200, mikä on paljon, mutta silti eulerin menetelmä "eksyy" pahasti kahdesta tarkemmasta. Tässä tehtävässä eksyminen johtuu maahan osumisen aiheuttamista suurista kiihtyvyyksistä ilmalentoon verrattuna. Tämä näkyy parhaiten vaaka- ja pystynopeuden nopeista muutoksista.
On olemassa algoritmeja, jotka arvioivat integroinnin virhettä ja sovittavat integrointiaskeleen sen mukaan.
otsikot = [
"vaakanopeus ajan funtiona", "pystynopeus ajan funktiona",
"etaisyys tykista ajan funktiona", "korkeus ajan funktiona"
]
for j in range(0, 4):
plt.plot(tt, zza[:, j], 'r-')
plt.plot(tt, zzb[:, j], 'b-')
plt.plot(tt, zzc[:, j], 'g-')
plt.grid(True)
plt.title(otsikot[j], color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
for j in range(0, 2):
plt.plot(tt, zza[:, j], 'r-')
plt.plot(tt, zzb[:, j], 'b-')
plt.plot(tt, zzc[:, j], 'g-')
plt.xlim([3.5, 8.0])
plt.ylim([-20, 12])
plt.grid(True)
plt.title(otsikot[j], color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
vv = [sqrt(zza[t, 0]**2 + zza[t, 1]**2) for t in range(0, dim)]
plt.plot(zza[:, 2], vv, 'r-')
vv = [sqrt(zzb[t, 0]**2 + zzb[t, 1]**2) for t in range(0, dim)]
plt.plot(zzb[:, 2], vv, 'b-')
vv = [sqrt(zzc[t, 0]**2 + zzc[t, 1]**2) for t in range(0, dim)]
plt.plot(zzc[:, 2], vv, 'g-')
plt.grid(True)
plt.title("nopeus paikan funktiona", color='green')
plt.tick_params(axis='both', colors='green')
# plt.legend(loc='upper left')
plt.show()
vv = [sqrt(zza[t, 0]**2 + zza[t, 1]**2) for t in range(0, dim)]
plt.plot(zza[:, 2], vv, 'r-')
vv = [sqrt(zzb[t, 0]**2 + zzb[t, 1]**2) for t in range(0, dim)]
plt.plot(zzb[:, 2], vv, 'b-')
vv = [sqrt(zzc[t, 0]**2 + zzc[t, 1]**2) for t in range(0, dim)]
plt.plot(zzc[:, 2], vv, 'g-')
plt.xlim([60, 75])
plt.ylim([0, 20])
plt.grid(True)
plt.title("nopeus paikan funktiona", color='green')
plt.tick_params(axis='both', colors='green')
#plt.legend(loc='upper left')
plt.show()
Lasketaan lentoradat kolmella eri lähtökulmalla
def lento(roo, korotus, v0):
alpha = korotus * (2.0 * pi / 360.0)
vx0 = cos(alpha) * v0
vy0 = sin(alpha) * v0
xx = [0 for i in range(0, dim)]
xx[0] = [vx0, vy0, x0, y0]
for i in range(0, dim - 1):
xx[i + 1] = rk4(fdxdt, xx[i], [roo], dt)
return np.array(xx)
g = 9.81
m = 1.0
x0 = 0.0
y0 = 0.1
T0 = 10.0
dt = 0.01
dim = int(T0 / dt)
tt = [i * dt for i in range(0, dim)]
zza = lento(0.02, 35, 50)
zzb = lento(0.02, 45, 50)
zzc = lento(0.02, 55, 50)
plt.plot(zza[:, 2], zza[:, 3], 'r-')
plt.plot(zzb[:, 2], zzb[:, 3], 'b-')
plt.plot(zzc[:, 2], zzc[:, 3], 'g-')
plt.grid(True)
plt.title("korkeus paikan funktiona", color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
plt.plot(zza[:, 2], zza[:, 3], 'r-')
plt.plot(zzb[:, 2], zzb[:, 3], 'b-')
plt.plot(zzc[:, 2], zzc[:, 3], 'g-')
plt.xlim((55, 85))
plt.ylim((-5, 10))
plt.grid(True)
plt.title("korkeus paikan funktiona", color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
vv = [sqrt(zza[t, 0]**2 + zza[t, 1]**2) for t in range(0, dim)]
plt.plot(zza[:, 2], vv, 'r-')
vv = [sqrt(zzb[t, 0]**2 + zzb[t, 1]**2) for t in range(0, dim)]
plt.plot(zzb[:, 2], vv, 'b-')
vv = [sqrt(zzc[t, 0]**2 + zzc[t, 1]**2) for t in range(0, dim)]
plt.plot(zzc[:, 2], vv, 'g-')
plt.grid(True)
plt.title("nopeus paikan funktiona", color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
vv = [sqrt(zza[t, 0]**2 + zza[t, 1]**2) for t in range(0, dim)]
plt.plot(zza[:, 2], vv, 'r-')
vv = [sqrt(zzb[t, 0]**2 + zzb[t, 1]**2) for t in range(0, dim)]
plt.plot(zzb[:, 2], vv, 'b-')
vv = [sqrt(zzc[t, 0]**2 + zzc[t, 1]**2) for t in range(0, dim)]
plt.plot(zzc[:, 2], vv, 'g-')
plt.xlim((55, 80))
plt.ylim((0, 20))
plt.grid(True)
plt.title("nopeus paikan funktiona", color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
otsikot = [
"vaakanopeus ajan funtiona", "pystynopeus ajan funktiona",
"etaisyys tykista ajan funktiona", "korkeus ajan funktiona"
]
for j in range(0, 4):
plt.plot(tt, zza[:, j], 'r-')
plt.plot(tt, zzb[:, j], 'b-')
plt.plot(tt, zzc[:, j], 'g-')
plt.grid(True)
plt.title(otsikot[j], color='green')
plt.tick_params(axis='both', colors='green')
plt.show()
