Funktioiden piirtäminen ja differentiaaliyhtälöiden numeerinen ratkaiseminen pythonilla

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 selaimeen, 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ää.

Funktio

funktio $f(x) = x*(x-1)*(x+1)$ pythonilla koodattuna

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)] 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

Monissa videopeleissä kappaleiden törmäysten seuraukset lasketaan ratkaisemalla differentiaaliyhtlöitä; savupilvien kehittyminen, meren aaltoilu saatetaan animoida ratkaisemalla osittaisdifferentiaaliyhtälöitä. kasvihuonekaasujen vaikutuksia maapallon ilmastoon selvitetään osittaisdifferentiaaliyhtälöiden avulla.

Vain todella yksinkertaisissa tapauksissa differentiaaliyhtälö voidaan ratkaista tarkasti matematiikan sääntöjen avulla. Yleensä kaikissa käytännön tapauksissa differentiaaliyhtälö joudutaan ratkaisemaan numeerisesti tietokoneella.

Numeerinen ratkaisu ei ole koskaan tarkasti oikea, eivätkä ratkaistavat yhtälötkään aina kuvaa tarkasti todellisuutta. Todellisia ongelmia ratkoessa pitää aina arvioida, kuvaavatko yhtälöt todellisuuden ilmiöitä riittävän hyvin ja onko ratkaisumenetelmä riittävän tarkka.

Vesisäiliö

Ämpäri
Kuva1: Vuotava vesisäiliö

Kun kuvan 1 vesisäiliön venttiilin asentoa muuttaa, säiliön pinnankorkeus ja ulosvirtaus säiliöstä muuttuvat. Miksi?

Säädä venttiiliä ja odota, että pinnankorkeus tasaantuu. Laita ulosvirtauksen alle mitta-astia ja katso kellosta, kauanko sen täyttyminen kestää. Kirjoita muistiin pinnankorkeus ja mitta-astian täyttymisaika. Toista tämä niin, että saat mittausarvot ainakin melkein tyhjälle, melkein täydelle ja puolillaan olevalle säiliölle. Aina parempi, jos ehdit tehdä mittauksia muillakin pinnankorkeuksilla.

Piirrä koordinaatistoon ulosvirtauksen $q_{out}$ nopeus (litraa sekunnissa) pinnakorkeuden $h$ (metrejä) funktiona.

Fysiikasta tiedämme, että virtaus putkessa on verrannollinen paine-eron neliöjuureen. Tarkistetaan, päteekö.

Piirrä samaan koordinaatistoon kuvaaja $ q_{out} = k \sqrt(h)$ Valitse k niin, että kuvaajat asettuvat mahdollisimman lähelle toisiaan.

Mistä säiliön pinnankorkeuden vaihtelu johtuu? Mitä tapahtuu, jos reikä säiliön pohjassa tukitaan ja venttiili avataan niin, että $q_{in}$ on esimerkiksi 1 l/s? Säiliön vesimäärä kasvaa litran sekunnissa. Jos säiliö on lieriön muotoinen ja sen pohjan pinta-ala on A, tilavuus on $V = Ah$ ja vastaavasti tilavuuden muutosnopeus on A kertaa pinnankorkeuden muutosnopeus. Jos pohjan pinta-ala olisi neliömetrin ja vettä tulisi säiliöön kuutiometri sekunissa, pinta nousisi metrin sekunnissa.

Pinnankorkeuden muutosvauhti, kun ulosvirtaus on tukittu, voidaan esittää differentiaaliyhtälöllä $$\frac{d}{dt}h(t) = q_{in}(t)$$

Ulosvirtauksen $q_{out}$ vaikutuksen pinnankorkeuteen voi järkeillä samalla periaatteella kuin sisäänvirtauksenkin. Kuten yllä todettiin $q_{out}(t) = k\sqrt(h)$ Kun ulostuloputki ei ole tukossa, saadaan differentiaaliyhtälö

$$\frac{d}{dt}h(t) = q_{in}(t) - k\sqrt(h(t))$$

Mitä voimme tuosta päätellä? Emme paljon mitään. Sen ehkä huomaa, että pinnankorkeuden muutosnopeus on nolla, kun virtaukset ovat yhtäsuuret. Mukava olisi nähdä, mitä tapahtuu ajan kuluessa, kun ventiiliä käännellään.

Numeerinen ratkaisu

Allaoleva ei ole matemaattisesti kunnollinen selitys, mutta auttaa toivottavasti mieltämään, mitä numeerinen integrointi tässä tapauksessa on. Matematiikan oppikirjoista kannattaa perehtyä derivaattaan ja differentiaaliyhtälöihin.

Jos pinnankorkeuden muutosnopeus $\frac{d}{dt}h(t)$ olisi vakio $\frac{d}{dt}h$, silloin pätisi $$h(t+\Delta t) = h(t) + \frac{d}{dt}h \Delta t $$

(Mieti, millä perusteella ylläoleva on totta.)

Koska ulosvirtaus riippuu pinnankorkeudesta, se muuttuu koko ajan pinnankorkeuden muuttuessa, eikä ylläoleva tarkkaan ottaen pidä paikkaansa. Jos kuitenkin $\Delta t$ on tarpeeksi pieni verrattuna pinnakorkeuden muutosnopeuteen, virhe on pieni. Tähän oletukseen perustuu yksinkertaisin tapa ratkaista differentiaaliyhtälöitä numeerisesti.

Eulerin algoritmi

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 heitettyä 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ä.

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)

Kirjoittamillani rutiineilla voi siis ratkaista usean muuttujan differentiaaliyhtälöryhmiä.

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.

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 nopeudelle 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ä.

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.

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

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.

Lasketaan lentoradat kolmella eri lähtökulmalla

Muuta ohjelmoinnista