Tietokone LUMA-opetuksen inspiroijana

import datetime as dt
print('Viimeksi muokattu: ', dt.date.today())

Viimeksi muokattu:  2024-09-05

Olen käynyt kouluni 60- ja 70-luvuilla ja opettajan pedagiset opinnot suoritin kymmenisen vuotta sitten. En siis tiedä, miten tietokoneita nykyään käytetään LUMA-opetuksessa.

Tarkoitukseni oli esitellä, miten tietokoneen avulla voi ratkaista käytännön esimerkkejä niin, että se innostaisi LUMA-aineiden opiskeluun.

Samalla tulisi esiteltyä, mihin kaikkeen myöhemmässä elämässä hyödylliseen ja mielenkiintoiseen LUMA-aineista erityisen kiinnostunut lukiolainen voisi perehtyä.

Innostuin kuitenkin heti ensimmäisestä esimerkistä liikaa, joten esitykseni jäi kapeaksi — eli siis perusteelliseksi ja syvälliseksi.

Pedagogisesti väärä tapa

Vuosia sitten minua jäi vaivaamaan Helsingin Sanomissa juttu, jossa kritisoitiin tietokoneen soveltamista matematiikan opetuksessa sillä perusteella, että kokeisiin vastaaminen LaTeX-ohjelman avulla on kömpelöä. Enää en löydä juttua, enkä tiedä, onko väite LaTeXin käytöstä.Mutta sen tiedän, että LaTeX on kynään ja paperiin verrattuna kömpelö eikä sen käyttäminen opeta matematiikasta mitään. LaTeX on tutkijoille sopiva matemaattisten tekstien julkaisuun kehitetty ladontaohjelma.

Lukijoiden kommentit toivat esiin hyviä näkökohtia matematiikasta ja tietokoneista, mutta myös sen, että harva tuntee tietokoneistetun matematiikan mahdollisuuksia ja haasteita eikä kukaan miettinyt, mitä matematiikasta oliisi oikeasti hyödyllistä oppia koulussa.

Tavoitteet ajalta ennen tietokoneita

Niin kauan kuin oppimistavoitteet on määritelty ja harjoitustehtävät laadittu kynä-ja-paperi 'aikakaudella', tietokoneesta on vähänlaisesti hyötyä.

Tietokoneen avulla voi monin tavoin havainnollistaa opittavia asioita ja tutkia sellaisia mielenkiintoisia "elävän elämän" ilmiöitä, joiden vaatima laskenta käsityönä puuduttaa istumalihakset ja aiheutta jännetupen tulehduksen.

Aloituskynnys tietokoneen hyödyntämisessä on kuitenkin melkoinen niin opiskelijalle kuin opettajallekin. 'Digiloikkaa' yritettiin ylhäältä tulleella käskyllä 'Loikatkaa'. Rahaa kului ja laitteita ostettiin, mutta lisääntyikö opiskelijoiden osaaminen?

Ohjelmoinnin periaatteiden tunteminen on hyödyllinen kansalaistaito. Ohjelmoinnin opetus kannattaa suunnata niin, että se tukee matematiikan ja fysiikan opetusta. Matematiikan ja fysiikan opetus rikastuu ja opiskelijat saavat harjoitusta ohjelmoinnin soveltamisessa.

Työkaluja

Tunnen rajoitetusti opetuksen tueksi sopivia työkaluja:

• Kaupallisia työkaluja en ole käyttänyt, koska avoimen lähdekoodin työkalut ovat tarjonneet kaiken tarvitsemani, vaikka olen tehnyt vaativampia juttuja kuin koulussa tarvitsee.
• Työni ja mielenkiintoni — kokemukseni — kattaa vain kapean alueen matematiikasta jafysiikasta. Minkä, se ilmenee myöhemmin tässä jutussa.

Macsymaa olen joskus käyttänyt ja sitä kai kouluissakin käytetään. Sopii symboliseen ja numeeriseen laskentaan ja kuvaajien piirtelyyn. Ihan hyvä työkalu. Notaatioltaan ja filosofialtaan "matematiikkalähtöinen". Tätä juttu arten installoin wxMaximan, jossa on graafinen käyttöliittymä: Menuista löytyy matemaattiset perusoperaatiot, kuten derivointi, integrointi, sieventäminen, yhtälöiden ratkaiseminen. Lausekkeita voi myös kirjoittaa samaan tyyliin kuin ohjelmointikielissä.

Maxima kokeilu

Käytän nykyään Jyputer-lab nimistä ympäristöä, koska se yhdistää symbolisen ja numeerisen laskennan Python ohjelmointikielen mahdollisuuksiin. Matemaatikon mielestä lopputulos saattaa vaikuttaa huonolta kompromissilta, mutta ohjelmoinnin opetuksen kannalta kompromissi on hyvä. Eikö ohjelmoinnin perusteiden ymmärtämistä pidetä nykyään olennaisena kansalaistaitona?

Käytännön haaste

Minkään työkalun käytön opettelu ei ole helppoa ja vie aikaa. Sama pätee ohjelmointiin. Aloittelijalta tietoteknisten työkalujen 'oikut' vievät helposti huomion itse asiasta.

Onko järkeä uhrata aikaa tietotekniikan käytön harjoitteluun, jos oppimisvaatimukset ovat ajalta ennen tietokoneita? Jos vaatimukset ja tehtävät ovat ajalta kyna-ja-paperi, opetukseen riittää kynä ja paperi. Kynä ja paperi on oivallinen työkalu perustehtävien suorittamiseen.

Minun mielestäni ohjelmointi kuuluu opetusohjelmaan. Se ei tarkoita, että koulun pitäisi kasvattaa koodareita vaan, että jokaisella olisi ymmärrys siitä, mitä ohjelmoimalla voi tehdä. Eikö matematiikan kouluopetusta perustella matemaattisen ajattelun edistämiselle. Ohjelmointi edistää sitä lisää. Eikä ymmärrystä mihinkään synny kuin tekemällä jotain konkreettista. Ratkomalla matikan tehtäviä, ohjelmoimalla.

Koulutuntien puitteissa ei ole mahdollista opiskella ohjelmointia niin, että ensin luetaan kirja ja sitten ohjelmoidaan. Itse haen esimerkin, jota muokkaamalla saan aikaiseksi, mitä haluan. Senkin oppiminen, miten löytää kuhunkin tarpeeseen sopivia esimerkkejä ja ohjeita vie aikaa.

Esimerkiksi Jupyter-lab:in käyttäminen matematiikan opetuksessa edellyttäisi valmiita malleja selityksineen siitä, miten tietyt perusasiat tehdään.

Tätä juttua laatiessa aika ajoin tuntui, että tällainen on liian vaativaa opittavaksi koulun tuntimäärien puitteissa. Toisaalta tuntuu, että opiskelijoille on ehdottomasti annettava mahdollisuus oppia näin mielenkiintoisia taitoja, joista monelle on jatkossakin hyötyä. Pitää löytää oikea tapa opettaa.

Miten välttää juuttuminen (ohjelmoinnin) perusteiden opiskeluun

Ennen kuin tietokoneesta on hyötyä ja hupia opiskelussa pitää sekä opettajan että opiskelijoiden nähdä vaivaa ohjelmoinnin opettelusta. Ohjelmoinnin opiskelu on osin kiehtovaa, osin puuduttavan tylsää.

Ohjelmointia voi kuitenkin hyödyntää LUMA-opetuksessa ilman, että ensin opiskellaan vuosi ohjelmointia. Sama pätee moniin sovelletun matematiikan menetelmiin Pitää vain kääntää ajatus opetuksesta ylösalaisin.

Kun aikoinaan TKK:lla kysyin, mitä järkeä tällaisia matemaattisia kieputuksia on opetella, minulle vastattiin: "Ymmärrät valmistuttuasi." Piti osin paikkansa, mutta motivaatio olisi ollut parempi, jos käymällä käsiksi käytännön ongelmaan olisi nähnyt, miksi nuo 'kieputukset' pitää osata ja niiden perusteet ymmärtää.

Onko matemaattisia menetelmiä, joita ei voi kokeilla käytännössä ennen kuin ymmärtää perusteet, eikä perusteita voi opettaa ennen kuin on opetettu perusteiden perusteet. Kun perusteiden perusteita on paljon, ei koskaan päästä opitun soveltamiseen asti. Oppi jää pölyttymään omaan lokeroonsa. "Matematiikasta ei ole hyötyä käytännössä!"

Entä jos ajateltaisiin, että tietokoneella voivat kirjoittaa vain ne, jotka ovat perehtyneet keskusyksikön toimintaan?

Matemaattisia menetelmiä voi kokeilla, vaikkei niiden matemaattista todistusta tuntisikaan. Kokeilu saattaa inspiroida perehtymään matemaattiseen taustaan.

Vaikkapa lentokoneita ei tietenkään pidä ruveta suunnittelemaan menetelmillä, joiden toimintaperiaatetta ja rajoituksia ei ymmärrä, mutta videopeleissä voi luonnonlakeja soveltaa vähemmällä varovaisuudella.

---

Ohjelmia laadittaessa tarvitaan matemaattista ymmärrystä ja oivallusta — matemaattista ajattelua, mutta ohjelmoinnissa tarvitaan myös paljon 'nippelitietoa' merkintätavoista. Sympyllä voi ratkaista matemaattisia tehtäviä kutsumalla sitä varten laadittuja funktioita, mutta loogisesti päättelemällä ei voi selvittää funktioiden nimiä saati niiden vaatimia parametrejä. Yritys-erehdys-menetelmällä olen oppinut, mistä tietoa kannattaa hakea. Ohjelmallisesti on kätevä piirtää kuvaajia, mutta vasta, kun on kaivanut esille millaisilla komennoilla piirtäminen tehdään.

Ainakin alkuvaiheessa nippelitieto on syytä tarjoilla valmiina. Niitä, jotka innostuvat ohjelmointia harrastamaan, voi myöhemmin opastaa selvittelemään ohjelmoinnin kiemuroita.

LUMA-opetusta varten on hyvä laatia esimerkkejä ratkaisuista ja/tai ohjelmarunkoja, joihin opiskelijan tarvitsee lisää vain kyseisen tehtävän LUMA-ajattelua vaativat osat.

---

LUMA-opetuksen tueksi sopivia työkaluja on monenlaisia.

Esittelen python-ohjelmointikielen tarjoamia mahdollisuuksia, koska harrastan simulaatiota ja sellaistan pelien ohjelmointia, jotka perustuvat lukion fysiikkaan ja matematiikkaan täydennettynä muutamalla käytännössä välttämättömällä lisäviisaudella.

Ohjelmointi on matemaattisten algoritmien toteuttamista tietokoneella. Se edellyttää ja kehittää matemaattista ajattelua, mutta puhtaan matemaattinen merkintätapa ei ohjelmoinnissa riitä. On pedagogisesti tärkeää huomata merkintätapojen ero.

---

Tätäkin esimerkkiä laatiessani työ keskeytyi monta kertaa, kun ohjelma ei toiminut niin kuin halusin. Loogisia ajatusvirheitään on toisinaan vaikea havaita.

Aloittelija tekee paljon virheitä, joiden ratkaiseminen omin avuin on niin työlästä, että se vie helposti innostuksen ohjelmoinnin opiskeluun. Opettajatkin voi joutua koville mystisiä ongelmia selvitellessä.

'Ilmiö' tarkasteltavaksi — vapaasti lentävä kappale

Tietokoneella voi ratkaista ongelmia, joita ei kynällä ja paperilla kykene ratkaisemaan. Tietokoneella voi ratkaista mielenkiintoisia 'elävän elämän' ongelmia. Näyttää, että matematiikasta on hyötyä.

Seuraavassa ratkaisen kappaleen liikeyhtälön tietokoneen avulla. Ennen kuin ryhtyy moiseen, opiskelijolla pitää tietysti olla alustava ymmärrys mekaniikasta ja differentiaaliyhtälöistä. Trigonometriaakin tarvitaan. Yksinkertaisin esimerkki kannattaa ratkaista ensin kynällä ja paperilla. Tietokoneella ei pidä temppuja, joista ei mitään ymmärrä.

70-luvun lukio-opetuksesta muistan, että jotkut fysiikan ilmiöt selitettiin hyvin kökösti, koska "matematiikassa ei ole vielä käsitelty tarvittavia menetelmiä". Mielenkiintoisia käytännön esimerkkejä ei voitu ratkaista, koska matematiikassa ei ...

Eikö opetuksessa voitaisi ottaa tarkasteltavaksi vaikkapa tykillä ampumisen kaltainen 'ilmiö'. Tarkastellaan kuulaa, joka singotaan taivaalle tietyllä alkunopeudella tietyssä kulmassa. Jatkossa voi laskea, miten jousi sinkoaa kuulan taivaalle. Tai ruuti.

Mitä tapahtuu? Aikansa lennettyään kuula putoaa maahan. Miksi? ...

Newtonin laki, differentiaaliyhtälö, x-y koordinaatisto, komponentteihin jako,

Entä ilmanvastus? Se vaikuttaa paljon, mutta sitä ei voinut lukiossa tarkastella, koska yhtälö ei ratkea analyyttisesti eikä numeerinen matematiikka kuuluu tavoitteisiin.

Eli mikään mielenkiintoinen ja käytännössä hyödyllinen ei kuuluu oppimäärään.

---

On tärkeää tehdä ero matematiikan ja 'insinöörilaskennon' välillä. Käytännön ongelmia ratkaistaan insinöörilaskennolla — derivoidaan, integroidaan, ratkotaan kolmioiden kulmia, ratkaistaan differentiaaliyhtälöitä jne.

Matematiikkaa taas ... No, sen selittämiseen, mitä matematiikka on, tarvitaan matemaatikkoa. Matemaatikot ovat kehittäneet insinöörilaskennan menetelmät ja todistaneet, että menetelmät ovat päteviä — matemaatikkojen määrittelemin rajoituksin.

Joitain opiskelijoita kiinnostaa matematiikka, toisia insinöörilaskennan soveltaminen käytännön esimerkkeihin.

Matemaattisen tekstin tuottaminen tietokoneella

Voin kirjoittaa tietokoneelle

$$ \underline{F}(t) = m \underline{a}(t) \iff \begin{bmatrix} F_x(t) \\ F_y(t) \end{bmatrix} = \begin{bmatrix} a_x(t) \\ a_y(t) \end{bmatrix} = \begin{bmatrix} \dot{v}_x(t) \\ \dot{v_y}(t) \end{bmatrix} $$ $$ \begin{bmatrix} v_x(t) \\ v_y(t) \end{bmatrix} = \begin{bmatrix} \dot{x}(t) \\ \dot{y}(t) \end{bmatrix} $$ $\underline{F}(t)$: Kappaleeseen vaikuttava voima $m$: Kappaleen massa $\underline{a}(t)$: Kappaleen kiihtyvyys $s_x, s_y$: Muuttujan $s$ x ja y-komponentit $\dot{s}(t)$: Muuttujan $s$ derivaatta ajan suhteen

ja saan näytölle tällaista

$$ \underline{F}(t) = m \underline{a}(t) \iff \begin{bmatrix} F_x(t) \\ F_y(t) \end{bmatrix} = \begin{bmatrix} a_x(t) \\ a_y(t) \end{bmatrix} = \begin{bmatrix} \dot{v}_x(t) \\ \dot{v_y}(t) \end{bmatrix} $$

$$ \begin{bmatrix} v_x(t) \\ v_y(t) \end{bmatrix} = \begin{bmatrix} \dot{x}(t) \\ \dot{y}(t) \end{bmatrix} $$ $\underline{F}(t)$: Kappaleeseen vaikuttava voima

$m$: Kappaleen massa

$\underline{a}(t)$: Kappaleen kiihtyvyys

$s_x, s_y$: Muuttujan $s$ x ja y-komponentit

$\dot{s}(t)$: Muuttujan $s$ derivaatta ajan suhteen

Entä sitten?

Onko vaivan arvoista? Joskus on:

• Tieteelliseen raporttiin en voi liittää kännykkäkuvaa käsin tuhertelemistani yhtälöistä.

• Jos joudun korjailemaan yhtälöitä, minun ei tarvitse kummata tai kirjoittaa kaikkea uudelleen.

Ajattelun tukena kynällä ja paperilla 'tuhertelu' on kuitenkin yleensä kätevämpää.

Mutta yhtälöt pitää myös ratkaista ja usein ne ovat valtavan työläitä ratkaista kynällä ja paperilla ja jos teen virheen, koko homma pitää aloittaa alusta.

Kynällä ja paperilla voi ratkaista vain niin pelkistettyjä ongelmia, että käytännön esimerkeistä pelkistyy käytäntö pois.

Monimutkaisten matemaattisten ongelmien ratkaisussa tietokone on välttämätön.

Ylläolevasta merkinnästä ei kuitenkaan ole apua ongelman ratkaisussa. Kannattaako mokomaa merkintätapaa opetella?

Symbolisen matematiikan työkalut osaavat näyttää matemaattiset lausekkeet tyylikkäästi, mutta sen lisäksi niillä voi ratkaista matemaattisia ongelmia.

Esimerkiksi Mathematica, Macsyma ja python-ohjelmointikielen lisäosa sympy, jota esittelen seuraavassa.

Ohjelmointia

Matemaattisten ongelmien ratkaisemiseen tietokoneella tarvitaan muutakin kuin yhtälöiden kirjoittaminen matematiikan sääntöjen mukaan. Työkalusta riippuen yhtälöitäkään ei — valitettavasti — voi kirjoittaa matematiikasta tutulla tavalla.

Opetuksessa monet ohjelmassa tarvittavat määrittelyt kannattaa antaa valmiina. Opiskelijoille voi antaa eräänlaisen työkirjan, johon täytetään matemaattista ymmärrystä vaativat kohdat. Joskus on järkevintä antaa esimerkki, joka muokataan käsillä olevaan ongelmaan sopivaksi.

Esimerkiksi allaolevasta riittää aluksi todeta "selitän myöhemmin". Tällainen lienee vastoin perinteistä ajattelua (matematiikan) opetuksesta. "Elämä on"

import math
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

Kappaleen liikeyhtälö xy-koordinaatistossa

Ratkaistaan xy-koordinaatistossa annetulla alkunopeudella lentoon lähetetyn kuulan lentorata.

Määritellään symbolisia muuttujia ja funktioita.

Sympyllä et voi kirjoittaa yhtälöä $ax² + bx + c = 0$ ja ratkaista sitä kertomatta, että $a, b, c$ ja $x$ ovat symbolisia muuttujia.

Funktiota $f(x) = ax² + bx + c$ et voi derivoida $x$:n suhteen ellet kerro, että $f$ on funktio.

Miksi? Pitkään miettimällä ja dokumentaatiota lukemalla se selviää, mutta sympy voi käyttää, vaikkei tätä(kään) miettisi.

t, m, g, v_x0, v_y0 = sp.symbols('t, m, g, v_x0, v_y0')
y = sp.Function('y')
x = sp.Function('x')
v_x = sp.Function('v_x')
v_y = sp.Function('v_y')
a_x = sp.Function('a_x')
a_y = sp.Function('a_y')
F_x = sp.Function('F_x')
F_y = sp.Function('F_y')

$F = ma$ on toisen kertaluvun differentiaaliyhtälö. Muutan sen ensimmäisen kertaluvun yhtälöryhmäksi syistä, jotka selitän myöhemmin, kun otamme huomioon ilmanvastuksen.

Alla käytettävä merkintä poikkeaa 'pahasti' matemaattisesta merkinnästä:

A = B tarkoittaa, että B:n arvo sijoitetaan muuttujalle A.

sp.diff(s(x),x) on sympy-kirjaston funktio, joka muodostaa derivaatan $\frac{ds(x)}{dx}$ lausekkeen

sp.Eq(A, B) muodostaa yhtälön $A=B$ lausekkeen

[a, b] on lista, johon voidaan koota yhteen kuuluvat muuttujat, lausekkeet tai yhtälöt

Yhtälöt kirjoitetaan sympy-kirjaston vaatimassa muodossa, mutta funktio display() muotoilee ne ihmiselle helpommin ymmärrettäviksi, joten lausekkeita kirjoittaessaan saa keskittyä tarkoitukseen ja tietokone näyttää, mitä tuli kirjoitettua.

Liikeyhtälöt yleisessä muodossa

# eq_a = sp.Eq([(a_x(t), sp.diff(v_x(t),t)), (a_y(t), sp.diff(v_y(t),t))])
a_x = sp.diff(v_x(t),t)
a_y = sp.diff(v_y(t),t)
eqs_g = [
    sp.Eq(v_x(t), sp.diff(x(t),t)), 
    sp.Eq(v_y(t), sp.diff(y(t),t)),
    sp.Eq(F_x(t), m*a_x),
    sp.Eq(F_y(t), m*a_y)
]

for eq in eqs_g:
    display(eq)

$\displaystyle v_{x}{\left(t \right)} = \frac{d}{d t} x{\left(t \right)}$

$\displaystyle v_{y}{\left(t \right)} = \frac{d}{d t} y{\left(t \right)}$

$\displaystyle F_{x}{\left(t \right)} = m \frac{d}{d t} v_{x}{\left(t \right)}$

$\displaystyle F_{y}{\left(t \right)} = m \frac{d}{d t} v_{y}{\left(t \right)}$

Ratkaistavat muuttujat

funcs = [x(t), y(t), v_x(t), v_y(t)]

Sympy-kirjaston dokumentaatio on mielestäni vaikeasti navigoitavaa ja vaikeaselkoista, mutta matemaatikoillekin kai riittävän seikkaperäinen.

Opiskelijoille on vähintäänkin annettava valmis linkki dokumentaation oikeaan kohtaan. Aluksi lienee tarpeen antaa esimerkki tai parametrejä vaille valmis koodinpätkä.

Lentoradan ratkaiseminen

Vaakasuorassa suunnassa kuulaan ei vaikuta voimia. y-suunnassa Maa vetää kuulaa puoleensa voimalla $-mg$

Sijoitetaan yleiseen liikeyhtälöön voimille arvot $F_x(t) = 0, F_y(t) = -mg$

Muuttujan korvaaminen, subs

# Alla muodostan listan seuraavalla python-maisella tavalla

xx = [i for i in range(10)]
yy = [n**2 for n in xx]
display(xx, yy)

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

[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

F_subs = [(F_x(t), 0), (F_y(t), -m*g)]

eqs = [eq.subs(F_subs) for eq in eqs_g]

for eq in eqs:
    display(eq)

$\displaystyle v_{x}{\left(t \right)} = \frac{d}{d t} x{\left(t \right)}$

$\displaystyle v_{y}{\left(t \right)} = \frac{d}{d t} y{\left(t \right)}$

$\displaystyle 0 = m \frac{d}{d t} v_{x}{\left(t \right)}$

$\displaystyle - g m = m \frac{d}{d t} v_{y}{\left(t \right)}$

Ratkaistaan liikerata

Differentiaaliyhtlöiden ratkaiseminen

liikerata = sp.dsolve(
    eqs,   # ratkaistava differentiaaliyhtälöryhmä 
    funcs, # ratkaistavat muuttujat
    # reunaehdot: kuula lähtee origosta nopeudella v = (v_x0, v_y0)
    ics={x(0): 0, y(0): 0, v_x(0): v_x0, v_y(0): v_y0})

for fun in liikerata:
    display(fun)

$\displaystyle x{\left(t \right)} = t v_{x0}$

$\displaystyle y{\left(t \right)} = - \frac{g t^{2}}{2} + t v_{y0}$

$\displaystyle v_{x}{\left(t \right)} = v_{x0}$

$\displaystyle v_{y}{\left(t \right)} = - g t + v_{y0}$

Lentoaika saadaan ratkaisemalla, milloin kuula on maan tasalla eli milloin $y(t) = 0$

$y(t)$:n lauseke on liikerata[1] eli listan liikerata toinen alkio

Yhtälön ratkaiseminen: Solve Equations

T_f = sp.solve(liikerata[1].rhs,t)
display(T_f)

[0, 2*v_y0/g]

T_f = T_f[1]
display(T_f)

$\displaystyle \frac{2 v_{y0}}{g}$

Piiretään lentorata eri alkuarvoilla

piirtäminen yleensä

subplots-funktio

plot-funktio

def vektorin_komponentit(pituus, suunta_asteina):
    kulma_radiaaneina = suunta_asteina*math.pi/180.0
    return (pituus*math.cos(kulma_radiaaneina),
            pituus*math.sin(kulma_radiaaneina))

fig, ax = plt.subplots()
# Kun aspect ratio on equal, lentorata pysyy luonnollisen muotoisena
ax.set_aspect('equal', 'box')
N = 20  # lasketaan lentorata N:ssä pisteessä
V = 20.0 # alkunopeus

# Piirretään kuvaaja kolmen eri lähtökulman arvolla
for suunta_kulma in (30.0, 45.0, 60.0):
    (vx_alku, vy_alku) = vektorin_komponentit(V, suunta_kulma)
    sijoitukset = [(v_x0, vx_alku), (v_y0, vy_alku), (g, 9.81), (m, 1.0)]

    # Sijoitetaan liikeradan funktioiden x(t) ja y(t) lausekkeille numeeriset arvot
    xn = liikerata[0].rhs.subs(sijoitukset)
    yn = liikerata[1].rhs.subs(sijoitukset)

    # Lasketaan N aikapisteessä
    dt = (T_f/N).subs(sijoitukset)

    # sijoitetaan arvo ajalle kussakin pisteessä
    # ja tehdään lista x- ja y-arvoista 
    xx = [xn.subs(t,dt*i) for i in range(N+1)]
    yy = [yn.subs(t,dt*i) for i in range(N+1)]
    ax.plot(xx, yy, marker='o', label='kulma ' + str(suunta_kulma))
    ax.legend()

Entä ilmanvastus?

Ilmanvastuksen aiheuttama voima on vastakkainen nopeuteen nähden. Täytyy piirrellä vektoreita, että osaa jakaa ilmanvastuksen x- ja y-komponentteihin.

Minulta piirtäminen sujuu parhaiten käsin, mutta aina en piirustuksistani tunnista samanmuotoisia kolmioita, ristikulmia, yms. Eikä tuherruksiani tietenkään kehtaa liittää mihinkään artikkeliin

Jotkut ehkä tuntevat käteviä geometrian työkaluja. Pitäisi kai minunkin niihin perehtyä.

Tavalliset piirtelytyökalut käyvät hätätilassa, mutta ovat aika kömpelöitä varsinkin, jos tulee tarvetta korjata piirustusta.

Seuraavassa esittelen kömpelöimmän tavan piirtää vektoreita. Ajattelun — ratkaisun hahmottamisen — tueksi tästä ei ole. Etuna tässä parametrointi, eli voin kokeilla, miltä eri suuntiin osoittavat vektorit näyttävät. Voin piirtää kätevästi eri kuvat tilanteesta, jossa kuula nousee ja tilanteesta, jossa se laskee

## Vektoreita

x1, x2, y1, y2 = sp.symbols("x1, x2, y1, y2")
vektori1 = sp.Matrix([x1, y1])
vektori2 = sp.Matrix([x2, y2])
display(vektori1, vektori2, vektori1 + vektori2)

$\displaystyle \left[\begin{matrix}x_{1}\\y_{1}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}x_{2}\\y_{2}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}x_{1} + x_{2}\\y_{1} + y_{2}\end{matrix}\right]$

def piirra_nuoli(ax, alkupiste, siirtyma, vari='black', nimi='', apunuoli=False, paksuus=0.1):
    x0, y0 = float(alkupiste[0]), float(alkupiste[1])
    dx, dy = float(siirtyma[0]), float(siirtyma[1])
    if apunuoli:
        vari = 'blue'
    # else:
    #    vari = 'black'
    ax.arrow(x0, y0, dx, dy, color=vari, fill=not apunuoli, 
             length_includes_head=True, width=paksuus)
    if nimi != '':
        ax.text(x0+dx/2.0, y0+dy/2.0, nimi, backgroundcolor='white', fontsize='small')

def piirra_kuula(ax, alkupiste, koko):
    p0 = (alkupiste[0], alkupiste[1])
    circle1 = plt.Circle(alkupiste, koko, color='black')
    ax.add_artist(circle1)     

alkupiste = sp.Matrix([1.0, 2.0])
vektori_x = sp.Matrix([20.0, 0.0])
vektori_y = sp.Matrix([0.0, 10.0])
display(alkupiste, vektori_x, vektori_y)

$\displaystyle \left[\begin{matrix}1.0\\2.0\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}20.0\\0\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}0\\10.0\end{matrix}\right]$

fig, axs = plt.subplots(3,3, figsize=[10, 10])
row = 0
for x0 in [-20, 10, 20]:
    col = 0
    for y0 in [-10, 5, 10]:
        vektori_x = sp.Matrix([x0, 0.0])
        vektori_y = sp.Matrix([0.0, y0])
        axs[row, col].set_aspect('equal', 'box')
        axs[row, col].set_axis_off()
        piirra_nuoli(axs[row, col], alkupiste, vektori_x, apunuoli=True, nimi='vektori_x')
        piirra_nuoli(axs[row, col],  alkupiste + vektori_x, vektori_y, apunuoli=True, nimi='vektori_y')
        piirra_nuoli(axs[row, col], alkupiste, vektori_x + vektori_y, nimi='vektori')
        piirra_kuula(axs[row, col], alkupiste, 0.5)
        col = col + 1
    row = row + 1

Ilmanvastus

c_r = 0.1
kuulan_paikka = (0.0, 0.0)
V = 20.0
alpha = 30.0
(vx, vy) = vektorin_komponentit(V, alpha)
nopeus = (vx, vy)

K = c_r*math.sqrt(vx**2 + vy**2)
Frx = -K*vx
Fry = -K*vy
Fr = (Frx, Fry)

fig, ax = plt.subplots()
ax.set_aspect('equal', 'box')
ax.set_axis_off()
piirra_nuoli(ax, kuulan_paikka, (vx, 0), apunuoli=True, nimi='v_x', vari='blue')
piirra_nuoli(ax, kuulan_paikka, (0, vy), apunuoli=True, nimi='v_y', vari='blue')
piirra_nuoli(ax, kuulan_paikka, nopeus, nimi='nopeus', paksuus=0.5, vari='green')

piirra_nuoli(ax, kuulan_paikka, (Frx, 0), apunuoli=True, nimi='Fr_x', vari='green')
piirra_nuoli(ax, kuulan_paikka, (0, Fry), apunuoli=True, nimi='Fr_y', vari='green')
piirra_nuoli(ax, kuulan_paikka, Fr, nimi='Fr', paksuus=0.5, vari='black')
piirra_kuula(ax, kuulan_paikka, 1.0)
ax.text(4.0, 0.5, "$\\alpha$")

Text(4.0, 0.5, '$\\alpha$')

Ilmanvastus on yleensä verrannollinen nopeuden neliöön: $Fr = c_r v^2$. Koska ilmanvastuksen aiheuttama voima vastustaa kappaleen liikettä sen suunta on nopeudelle vastakkainen.

Yhdenmuotoisista kolmioista:

$$ \frac{Fr}{v} = c_r \cdot v \iff \frac{Fr_x}{v_x} = c_r \cdot v ; \frac{Fr_y}{v_y} = c_r \cdot v$$

$$ Fr_x = c_r \cdot \sqrt{v_x^2 + v_y^2} \cdot v_x ; Fr_y = c_r \sqrt{v_x^2 + v_y^2} \cdot v_y $$

Ilmanvastuksen yhtälö

c_r, alpha, m, Fr_x, Fr_y, g = sp.symbols("c_r, alpha, m, Fr_x, Fr_y, g", real=True)
# Fr_x = sp.Function('Fr_x')
# Fr_y = sp.Function('Fr_y')
v_x = sp.Function('v_x')
v_y = sp.Function('v_y')
x = sp.Function('x')
y = sp.Function('y')

def nopeus_xy(nopeus, alpha): kulma_radiaaneina = alpha*math.pi/180.0 return (nopeus*math.cos(kulma_radiaaneina), nopeus*math.sin(kulma_radiaaneina))def vektori_diff(vektori, t, n=1): return [diff(f, t, n) for f in vektori]

Fr = c_r*(v_x(t)**2 + v_y(t)**2)
eqx = sp.Eq(Fr_x/Fr, -v_x(t)/sp.sqrt(v_x(t)**2 + v_y(t)**2))
eqy = sp.Eq(Fr_y/Fr, -v_y(t)/sp.sqrt(v_x(t)**2 + v_y(t)**2))
sol = sp.solve([eqx, eqy], [Fr_x, Fr_y], dict=False)
display('Fr:', Fr, 'eqx:', eqx, 'eqy:', eqy, sol)
display(sol[Fr_x], sol[Fr_y])

'Fr:'

$\displaystyle c_{r} \left(v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}\right)$

'eqx:'

$\displaystyle \frac{Fr_{x}}{c_{r} \left(v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}\right)} = - \frac{v_{x}{\left(t \right)}}{\sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}}}$

'eqy:'

$\displaystyle \frac{Fr_{y}}{c_{r} \left(v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}\right)} = - \frac{v_{y}{\left(t \right)}}{\sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}}}$

{Fr_x: -c_r*sqrt(v_x(t)**2 + v_y(t)**2)*v_x(t),
 Fr_y: -c_r*sqrt(v_x(t)**2 + v_y(t)**2)*v_y(t)}

$\displaystyle - c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{x}{\left(t \right)}$

$\displaystyle - c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{y}{\left(t \right)}$

Liikeyhtälöt ilmanvastus huomioituna

Fr_x = sol[Fr_x]
Fr_y = sol[Fr_y]
XY = sp.Matrix([x(t), y(t)])
V = sp.Matrix([v_x(t), v_y(t)])
F = sp.Matrix([Fr_x, Fr_y - m*g])
display(XY, V, F)

$\displaystyle \left[\begin{matrix}x{\left(t \right)}\\y{\left(t \right)}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}v_{x}{\left(t \right)}\\v_{y}{\left(t \right)}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}- c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{x}{\left(t \right)}\\- c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{y}{\left(t \right)} - g m\end{matrix}\right]$

eq1 = sp.Eq(sp.diff(XY, t), V)
eq2 = sp.Eq(sp.diff(V, t), F/m)
display(eq1, eq2)

$\displaystyle \left[\begin{matrix}\frac{d}{d t} x{\left(t \right)}\\\frac{d}{d t} y{\left(t \right)}\end{matrix}\right] = \left[\begin{matrix}v_{x}{\left(t \right)}\\v_{y}{\left(t \right)}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}\frac{d}{d t} v_{x}{\left(t \right)}\\\frac{d}{d t} v_{y}{\left(t \right)}\end{matrix}\right] = \left[\begin{matrix}- \frac{c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{x}{\left(t \right)}}{m}\\\frac{- c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{y}{\left(t \right)} - g m}{m}\end{matrix}\right]$

eqm = sp.Eq(sp.Matrix.vstack(eq1.lhs, eq2.lhs), sp.Matrix.vstack(eq1.rhs, eq2.rhs))
dSdt0 = sp.Matrix.vstack(eq1.rhs, eq2.rhs)

Kuulan paikan ja nopeuden derivaatat

dSdt = sp.Matrix([sp.simplify(expr) for expr in dSdt0])
display(eqm, dSdt)

$\displaystyle \left[\begin{matrix}\frac{d}{d t} x{\left(t \right)}\\\frac{d}{d t} y{\left(t \right)}\\\frac{d}{d t} v_{x}{\left(t \right)}\\\frac{d}{d t} v_{y}{\left(t \right)}\end{matrix}\right] = \left[\begin{matrix}v_{x}{\left(t \right)}\\v_{y}{\left(t \right)}\\- \frac{c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{x}{\left(t \right)}}{m}\\\frac{- c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{y}{\left(t \right)} - g m}{m}\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}v_{x}{\left(t \right)}\\v_{y}{\left(t \right)}\\- \frac{c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{x}{\left(t \right)}}{m}\\- \frac{c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{y}{\left(t \right)}}{m} - g\end{matrix}\right]$

Tarvittaessa sympy-lausekkeista voi tuottaa python-koodia.

pycode kieltäytyy käsittelemästä muotoa x(t) olevaa lauseketta joten se on korvattava lausekkeella x[t]

Peleissä tai muuten simuloinnissa tarvitsee yleensä funktiota, joka laskee tilayhtälöiden derivaatan tietyllä ajanhetkellä eli seuraava kelpaa.

Funktiosta kannattaa tehdä yleispätevä

vs_x, vs_y = sp.symbols('vs_x, vs_y')
poista_aika = [(v_x(t), vs_x), (v_y(t), vs_y)]

koodi = [sp.pycode(dS.subs(poista_aika)) for dS in dSdt]

print('def dfdt(xx, params):')
print('    # x_s ja y_s tässä erikoistapauksessa turhia')
print('    (x_s, y_s, vs_x, vy_x) = xx')
print('    (c_r, m, g) = params')
print(f'    return {koodi}\n')

def dfdt(xx, params):
    # x_s ja y_s tässä erikoistapauksessa turhia
    (x_s, y_s, vs_x, vy_x) = xx
    (c_r, m, g) = params
    return ['vs_x', 'vs_y', '-c_r*vs_x*math.sqrt(vs_x**2 + vs_y**2)/m', '-c_r*vs_y*math.sqrt(vs_x**2 + vs_y**2)/m - g']

Simulointi

Eulerin menetelmä

Kokeillaan ensin Eulerin menetelmää

$\underline x(t+dt) = \underline x(t) + \underline{\dot x}(t) dt $

dsys_dt on funktio, joka palauttaa kuulan paikan ja nopeuden derivaattojen $\underline{\dot x}(t)$ numeeriset arvot annetussa pisteessä.

Sijoitetaan liikeyhtälön dSdt lausekkeeseen hetkeä t vastaavat numeeriset arvot

display(dSdt)

$\displaystyle \left[\begin{matrix}v_{x}{\left(t \right)}\\v_{y}{\left(t \right)}\\- \frac{c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{x}{\left(t \right)}}{m}\\- \frac{c_{r} \sqrt{v_{x}^{2}{\left(t \right)} + v_{y}^{2}{\left(t \right)}} v_{y}{\left(t \right)}}{m} - g\end{matrix}\right]$

def dsys_dt(dSdt, S, params):
    [xn, yn, vn_x, vn_y] = S
    (cr, massa, g0) = params
    # tmp = dSdt.subs([(c_r, c_r), (m, m)])
    # display(tmp)
    tmp = [dS.subs([(g, g0), (v_x(t), vn_x), (v_y(t), vn_y), (c_r, cr), (m, massa)])
           for dS in dSdt]
#     print(f'\n\ndsys:\n {S}\n {tmp}\n {params}')
#    display(dSdt)
    return sp.Matrix(tmp)  # (x(t), xn), (y(t), yn), 

Alkutila

x0 = 0.0
y0 = 0.0
V0 = 20.0
alpha = 45.0
(vx0, vy0) = vektorin_komponentit(V0, alpha)
S0 = sp.Matrix([x0, y0, vx0, vy0])

fig, ax = plt.subplots()
ax.set_aspect('equal', 'box')

max_steps = 500
S = [0 for i in range(max_steps+1)]
xx = [0 for i in range(max_steps+1)]
yy = [0 for i in range(max_steps+1)]

cr_a = 0.2
m_a = 1.0
g_a = 9.81

# Kokeillaan eripituisia simulointiaskelia
for dt in [0.2, 0.1, 0.05, 0.01, 0.005]:
    S[0] = S0
    n = 0
    # jatketaan niin kauan kuin kuula on maan pinnan yläpuolelle,
    # kuitenkin enintään max_steps askelta
    while (n < max_steps) & (S[n][1] > -0.1):
        # Euler askel
        S[n+1] = S[n] + dsys_dt(dSdt, S[n], (cr_a, m_a, g_a))*dt
        n = n+1
    xx = [S[i][0] for i in range(n)]    
    yy = [S[i][1] for i in range(n)]    
    ax.plot(xx, yy, label=f'dt: {dt}; {n} askelta')  # , marker='o'
    ax.legend()

Ennustaja-korjaaja

Numeerisen simuloinnin tarkkuus riippuu integrointiaskeleen pituudesta. Ylläolevasta voisi päätellä, että mitä pidempi askel, sitä virheellisempi tulos.

Kovin lyhyt askel hidastaa simulointia, mutta pitäisi myös miettiä, aiheuttaako liian lyhyt askel virhettä.

Kokeillaan Eulerin menetelmää tarkempaa predictor-corrector algoritmia

Olen vertaillut algoritmeja tarkemmin toisaalla

def pred_corr(dSdt, fdydt, S, dt, params):
    # fdydt: funktio, joka laskee järjestelmän muuttujien derivaattojen arvot
    # S: piste, jossa derivaatat lasketaan
    # dt: integrointiaskel
    # params: lisäparametrejä

    # derivaatat askelen alkupisteessä 
    dS_a = fdydt(dSdt, S, params)
    
    # Ennuste Eulerin menetelmällä
    S_pred = S + dS_a*dt

    # derivaatat ennustetussa askelen loppupisteessä 
    dS_b = fdydt(dSdt, S_pred, params)

    # Lasketaan arvot loppupisteessä kuten Eulerin menetelmässä,
    # mutta käytetään derivaattojen keskiarvoa
    S_cor = S + (dS_a + dS_b) * dt/2.0
    return S_cor

def simu(dSdt, S0, dt, params, tunniste, debug=False):
    max_steps = 500
    S = [0 for i in range(max_steps+1)]
    xx = [0 for i in range(max_steps+1)]
    yy = [0 for i in range(max_steps+1)]

    N = 0
    S[0] = S0
    if debug:
        print(f'debug:, {N}, {S[N][1]}, {max_steps}')
        display(S[N])
        
    while (S[N][1] > -0.01) and (N < max_steps):
        S[N+1] = pred_corr(dSdt, dsys_dt, S[N], dt, params)  # S[N] + dsys_dt(dSdt, S[N], (cr, massa))*dt
        N = N+1
        
    xx = [S[i][0] for i in range(N)]    
    yy = [S[i][1] for i in range(N)]    
    ax.plot(xx, yy, label=tunniste + f'; {N} askelta')  # , marker='o'
    ax.legend()

max_steps = 500 S = [0 for i in range(max_steps+1)] xx = [0 for i in range(max_steps+1)] yy = [0 for i in range(max_steps+1)]

cr_a = 0.2
m_a = 1.0
x_0 = 0
y_0 = 0.0
V0 = 20.0
alpha = 45.0
(vx_0, vy_0) = vektorin_komponentit(V0, alpha)

fig, ax = plt.subplots()
ax.set_aspect('equal', 'box')

for dt in [0.3, 0.2, 0.1, 0.01, 0.005]:
    S0 = sp.Matrix([x_0, y_0, vx_0, vy_0])
    params = (cr_a, m_a, g_a)
    simu(dSdt, S0, dt, params, f'dt: {dt}')

Noin kymmenen kertaa pidemmällä aika-askeleella näytään päästävän samaan tarkkuuteen kuin Eulerin menetelmällä.

Numeerisen matematiikan asiantuntijat eivät tietenkään tyydy tällaiseen silmämääräiseen arviointiin.

Itse käytän yleensä rk4-menetelmää

Ilmanvastuksen vaikutus lentorataan

fig, ax = plt.subplots()
ax.set_aspect('equal', 'box')
dt = 0.02

for cr_a in [0.0, 0.01, 0.1, 0.2]:
    S0 = sp.Matrix([x_0, y_0, vx_0, vy_0])
    params = (cr_a, m_a, g_a)
    simu(dSdt, S0, dt, params, f'c_r: {cr_a}')

Kuulan massan vaikutus lentorataan

fig, ax = plt.subplots()
ax.set_aspect('equal', 'box')
dt = 0.02
cr_a = 0.2

for massa in [0.1, 1.0, 10.0]:
    S0 = sp.Matrix([x_0, y_0, vx_0, vy_0])
    params = (cr_a, massa, g_a)
    simu(dSdt, S0, dt, params, f'massa: {massa}')

Kun lähtönopeus on sama, isomman massan liike-energia on suurempi, joten ilmanvastuksen tekemä 'työ' hidastaa sitä vähemmän, mikäli ilmanvastus on sama raskaalla ja kevyellä kuulalla.

Optimaalinen lähtökulma

fig, ax = plt.subplots()
ax.set_aspect('equal', 'box')
dt = 0.02
cr_a = 0.2
massa = 1.0

for alpha in [30.0, 35.0, 45.0, 60.0]:
    (vx_0, vy_0) = vektorin_komponentit(V0, alpha)
    S0 = sp.Matrix([x_0, y_0, vx_0, vy_0])
    params = (cr_a, m_a, g_a)
    simu(dSdt, S0, dt, params, f'$\\alpha$: {alpha}')

Jos ilmanvastusta ei ole,

• kuula lentää pisimmälle, kun se ammutaan 45° kulmassa
• kuula lentää yhtä kauas 45°±x° lähtökulmilla

Näin ei selvästikään ole, kun ilmanvastus otetaan huomioon.

Jätän harjoistustehtäväksi laskea, millä lähtökulmalla kuula lentää pisimmälle.

Mistä alkuvauhti

Jousi

Tästäkin pitäisi piirtää kuva.

$K$: jousivakio, $L$: jousen pituus, $s(t)$: kappaleen paikka niin, että $L - s(t)$ on jousen venymä

K, massa, L, T_f, v_f = sp.symbols("K massa L T_f v_f", positive=True)
s = sp.Function("s")
v_s = sp.Function("v_s")

# massa = 1.0
# K = 100.0
# L = 2.0

eq_s = [sp.Eq(K*(L - s(t)), massa*sp.diff(v_s(t), t)), 
        sp.Eq(v_s(t), sp.diff(s(t), t))]
for eq in eq_s:
    display(eq)

$\displaystyle K \left(L - s{\left(t \right)}\right) = massa \frac{d}{d t} v_{s}{\left(t \right)}$

$\displaystyle v_{s}{\left(t \right)} = \frac{d}{d t} s{\left(t \right)}$

laukaisu = sp.dsolve(eq_s, ics={s(0):0, v_s(0):0})
for x in laukaisu:
    x = sp.simplify(x)
    display(x)

$\displaystyle v_{s}{\left(t \right)} = \frac{\sqrt{K} L \sin{\left(\frac{\sqrt{K} t}{\sqrt{massa}} \right)}}{\sqrt{massa}}$

$\displaystyle s{\left(t \right)} = L \left(1 - \cos{\left(\frac{\sqrt{K} t}{\sqrt{massa}} \right)}\right)$

eqs_f = [eq.subs([(t, T_f), (s(T_f), L), (v_s(T_f), v_f)]) for eq in laukaisu]
for eq in eqs_f:
    eq = sp.simplify(eq)
    display(eq)

$\displaystyle v_{f} = \frac{\sqrt{K} L \sin{\left(\frac{\sqrt{K} T_{f}}{\sqrt{massa}} \right)}}{\sqrt{massa}}$

$\displaystyle L = L \left(1 - \cos{\left(\frac{\sqrt{K} T_{f}}{\sqrt{massa}} \right)}\right)$

sol = sp.solve(eqs_f, [T_f, v_f], dict=True)
display(sol)

[{T_f: pi*sqrt(massa)/(2*sqrt(K)), v_f: sqrt(K)*L/sqrt(massa)}]

code = sp.pycode(sol[0][v_f])
print(code)

math.sqrt(K)*L/math.sqrt(massa)

def vf_jousi(K, L, massa):
    return math.sqrt(K)*L/math.sqrt(massa)

fig, ax = plt.subplots(figsize=[10, 10])
ax.set_aspect('equal', 'box')
dt = 0.01
cr = 0.2
massa = 1.0
g0 = 9.81
alpha = 45
L = 2.0
K = 100.0
x_0 = 0.0
y_0 = 0.0

for massa in [0.2, 1.0, 2.0, 5.0, 10.0]:
    V0 = vf_jousi(K, L, massa)
    (vx_0, vy_0) = vektorin_komponentit(V0, alpha)
    S0 = sp.Matrix([x_0, y_0, vx_0, vy_0])
    params = (cr, massa, g0)
    simu(dSdt, S0, dt, params, f'massa: {massa}; v0: {V0:.1f} ')

Jätän harjoitustehtäväksi arvioida, miksi tulos on mikä ja onko se järkevä.

Lisäharjoituksia

Runge-Kutta algoritmi

Vertaile allaolevaa aiemmin käytettyihin integrointialgoritmeihin

#lasketaan y(t+dt), kun dy(t)/dt = f(y(t),u(t))
# argumentit:
# dxdt:  funktio f, palauttaa tilamuuttujien derivaatat
# yy: y(t), järjestelmän tilamuuttujat
# uu: u(t), ulkoiset ohjaukset ja parametrit
# dt: integrointiaskel

# funktio palauttaa listan yy1, joka on y(t+dt)

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

Satelliitin ratakorjaukset

• Kirjoita polaarikoordinaatistossa yhtälöt satelliitin ratakorjauksille;
• simuloi, mitä tapahtuu, kun ratanopeutta kiihdytetään käyttämällä aika ΔT moottoria;
• piirrä simulaation tulos polar plottina
• Tee peliohjelma, jolla voi harjoitella satelliitin telakoitumista avaruusasemaan. Tässä saattaa olla hyötyä näistä esimerkeistä