8.1 Fra ODE til kode: Euler–Cromer-løkka

$$t_0 = 0, \quad t_1 = dt, \quad t_2 = 2\,dt, \quad \dots, \quad t_i = i \cdot dt$$

Tallet $i$ kalles indeksen, og det er den som teller hvilket tidspunkt vi står i. I koden svarer $t_i$ til t[i].

Å erstatte kontinuerlig tid med et slikt rutenett kalles å diskretisere. Ordet er ikke mer mystisk enn «å dele opp i skiver» — og det er hele grunnideen bak numerisk løsning av en differensiallikning.

Tidsavstanden mellom to nabotidspunkt i rutenettet, målt i sekunder.

Steglengden er den ene knappen du kan skru på i en numerisk beregning. Liten $dt$ gir et tettere rutenett og et mer nøyaktig svar, men flere steg og lengre regnetid. Stor $dt$ går fort, men svaret kan bli ubrukelig.

En brukbar tommelfingerregel for svingninger og kast er at $dt$ skal være et par tusen ganger mindre enn den karakteristiske tiden i problemet — svingeperioden, falltiden eller tidskonstanten. Har systemet en periode på ett sekund, er $dt = 0{,}001\ \text{s}$ et trygt valg.

Sammenhengen mellom antall punkt og sluttiden er $t_{\text{slutt}} = (N-1)\,dt$: velger du to av de tre størrelsene, er den tredje låst.

Tabeller som opprettes med full lengde før løkka starter, slik at hvert tidspunkt har sin ferdige plass å skrive svaret sitt inn i.

t = np.zeros(N)
x = np.zeros(N)
v = np.zeros(N)

np.zeros(N) lager en tabell med $N$ nuller. Etter løkka inneholder x posisjonen ved hvert av de $N$ tidspunktene, og du kan lese ut, plotte eller lete i hele forløpet.

Hvorfor gjøre det på forhånd? Fordi løkka da bare skal fylle inn, ikke utvide. Det er raskere, og — viktigere for eksamen — det tvinger fram indeksskrivemåten x[i], som er nettopp det sensor ser etter.

$$\frac{d^2x}{dt^2} = \frac{1}{m}\sum F \qquad \longrightarrow \qquad \texttt{a = ...}$$

Alt på høyresiden må være noe koden allerede kjenner: konstantene ($m$, $k$, $g$, $D$) og tilstanden i indeks $i$ (x[i] og v[i]).

Dette er de tre poengene. En akselerasjonslinje med riktig fortegn, riktige krefter og riktig masse-deling er selve fysikken i deloppgaven — resten av løkka er mal.

$$v_{i+1} = v_i + a_i\,dt, \qquad x_{i+1} = x_i + v_{i+1}\,dt$$

I kode:

v[i+1] = v[i] + a*dt
x[i+1] = x[i] + v[i+1]*dt

Metoden er oppkalt etter Alar Cromer, som beskrev den i 1981 som en «forbløffende enkel» rettelse av standardmetoden for mekanikk-simuleringer.

Dette er metoden dette emnet bygger på, og den sensor forventer. Grunnen er at den holder energien i sjakk over lange tider — se løkke 4.

$$x_{i+1} = x_i + v_i\,dt, \qquad v_{i+1} = v_i + a_i\,dt$$

I kode:

x[i+1] = x[i] + v[i]*dt
v[i+1] = v[i] + a*dt

Metoden er den mest opplagte oversettelsen av differensiallikningen, og den brukes mye i andre sammenhenger. For svingninger og baner er den systematisk dårlig: den tilfører energi til systemet steg for steg, slik at amplituden vokser og banene spiraler utover.

Legg merke til at forskjellen fra Euler–Cromer er én indeks. Bytter du rekkefølgen på de to linjene uten å endre indeksene, har du byttet metode.

Kravet om at hver størrelse i løkka skal skrives med den indeksen som sier hvilket tidspunkt den hører til: v[i], v[i+1], x[i], x[i+1].

Skriver du v = v + a*dt uten indeks, forsvinner nettopp den opplysningen løkka handler om. På eksamen koster det typisk 2 poeng, selv når meningen er åpenbar for en leser — og i en kode som faktisk skal kjøre, ødelegger det lagringen av forløpet, siden bare den siste verdien overlever.

Unntaket er akselerasjonen. a regnes ut på nytt i hvert eneste gjennomløp og brukes umiddelbart, så den trenger ikke lagres. Vil du likevel ha den som array, skriver du a[i] — det er også riktig, bare mer arbeid for hånd.

At den beregnede mekaniske energien i et system endrer seg systematisk gjennom simuleringen, selv om den fysisk skal være bevart.

Energidrift er en metodefeil, ikke fysikk. Symptomene er lette å kjenne igjen:

- en oscillator der amplituden vokser jevnt periode for periode;
- en planet- eller kometbane som spiraler utover i stedet for å lukke seg;
- en pendel som svinger høyere og høyere uten at noe tilfører den energi.

Forward Euler har energidrift oppover for svingninger og baner. Euler–Cromer har det ikke: energien svinger rundt riktig verdi uten å vandre av gårde. Det er derfor metoden er standard i mekanikk-simuleringer, og det er nettopp dette du skal kunne si i de to setningene som begrunner metodevalget.

Et spesialtilfelle av det samme problemet der du kjenner svaret fra en formel, og som du kjører koden på før du stoler på den.

For en kode som skal håndtere en fjær med luftmotstand, er kontrollcaset å sette motstanden til null: da må koden gi $x(t) = A\cos(\omega_0 t)$ med $\omega_0 = \sqrt{k/m}$. For en kode som skal håndtere et kast med drag, er kontrollcaset $D = 0$: da må rekkevidden bli $(v_0^2/g)\sin 2\theta$.

Regelen er enkel: slå av det leddet som gjør problemet vanskelig, og sjekk at koden gjenskaper det du kan regne ut for hånd. Gjør den ikke det, er feilen i infrastrukturen — og da hjelper det ikke å finpusse på det vanskelige leddet.

Å kjøre den samme koden på nytt med halvert steglengde og se om svaret endrer seg vesentlig.

Endrer svaret seg lite, er $dt$ liten nok. Endrer det seg mye, er den ikke det — og da forteller ingen enkeltkjøring deg hva som er riktig.

Testen skiller også metodefeil fra modellfeil: en feil som skyldes for grov diskretisering, krymper når $dt$ krymper. En feil som ligger i fysikken — et tapt kraftledd, et galt fortegn, en feil initialbetingelse — står helt uendret.

Dette er den mest undervurderte kontrollen som finnes, og den er verdt en setning i en besvarelse selv når du ikke har en maskin å kjøre på: «steglengden bør kontrolleres ved å halvere $dt$ og se at svaret er stabilt».

En metode der feilen i sluttsvaret er tilnærmet proporsjonal med steglengden: halverer du $dt$, halveres feilen.

Både Euler–Cromer og Forward Euler er førsteordens i denne forstand. For fjærpendelen i eksempel 1, målt mot den analytiske løsningen ved $t = 10{,}0\ \text{s}$:

De to tallene som beskriver starttilstanden, skrevet inn i indeks 0 av tabellene:

x[0] = ...      # hvor legemet er ved t = 0
v[0] = ...      # hvor fort det gaar ved t = 0

De er verdt 1 poeng i en typisk kodedeloppgave, og de er det eneste stedet i koden der opplysningene fra oppgaveteksten om starttilstanden dukker opp.

Den vanligste feilen er å sette begge til null av vane. «Slippes fra ro i utslaget $A$» gir x[0] = A, v[0] = 0; «passerer likevekt med farten $v_0$» gir x[0] = 0, v[0] = v0. Les setningen to ganger.

$$i = 0, 1, \dots, N-2 \quad \Longrightarrow \quad i+1 = 1, 2, \dots, N-1$$

Med range(N) ville siste gjennomløp prøvd å skrive til x[N], som ikke finnes, og programmet stopper. På papiret straffes det ikke, men det er et unødvendig signal om at løkka ikke er tenkt igjennom.

Huskeregelen: antall steg er alltid én mindre enn antall punkt. Med $N$ punkt går det $N-1$ mellomrom.

En metode der den nye tilstanden regnes ut direkte fra den gamle, uten at man må løse en likning underveis.

Både Euler–Cromer og Forward Euler er eksplisitte: høyresiden inneholder bare størrelser du allerede kjenner, og oppdateringen er ren innsetting.

Alternativet er implisitte metoder, der den nye tilstanden står på begge sider og må løses ut. De er mer stabile for enkelte problemer, men langt tyngre å skrive for hånd, og de er ikke pensum i dette emnet.

$$\frac{dv}{dt} = \frac{1}{m}\sum F \quad \longrightarrow \quad \texttt{v[i+1] = v[i] + a*dt}$$

$$\frac{dx}{dt} = v \quad \longrightarrow \quad \texttt{x[i+1] = x[i] + v[i+1]*dt}$$

Å se koblingen mellom de to formene er selve broa fra kap. 7.1 til dette kapitlet: hver førsteordens likning i tilstandsformen blir én oppdateringslinje.

Kode skrevet for å bli lest av et menneske, ikke kjørt av en maskin: riktig struktur og riktige uttrykk, men uten krav om at hvert komma sitter.

På eksamen er dette formen som forventes. Syntaksfeil straffes ikke, alle programmeringsspråk godtas, og du skriver ikke et helt program — bare initialbetingelser og løkke.

Det som likevel må være presist, er indeksene og rekkefølgen: v[i+1] før x[i+1], og den oppdaterte farten i posisjonslinja. Der er pseudokode like streng som ekte kode.

Feilen som oppstår i ett enkelt tidssteg, fordi vi behandler akselerasjonen som konstant gjennom steget mens den i virkeligheten endrer seg.

Feilen kommer av at vi kutter Taylor-rekka for posisjonen etter førsteordensleddet: det neste leddet er $\tfrac12 a\,dt^2$, og det er størrelsesordenen på feilen per steg.

Feilen per steg går altså som $dt^2$, mens den samlede feilen etter $t/dt$ steg går som $dt$. Det er nettopp derfor metodene kalles førsteordens, selv om enkeltsteget er andreordens nøyaktig.

$$t_{\text{slutt}} = (N-1)\,dt$$

Velger du hvor lenge simuleringen skal vare og hvor fint du vil regne, er antall punkt bestemt. Vil du dekke fem svingeperioder à $1{,}26\ \text{s}$ med $dt = 0{,}0010\ \text{s}$, trenger du $N = 6284$.

Den praktiske konsekvensen er at finere oppløsning koster proporsjonalt med regnetid: halverer du $dt$, dobles antall steg for samme sluttid.

Et krav om at løkka skal avsluttes når noe fysisk inntreffer, i stedet for etter et fast antall steg.

For et kast er det naturlige kriteriet at legemet treffer bakken:

if y[i+1] < 0:
    break

Alternativet er å velge $N$ stort nok og lete opp nedslaget i etterkant. Begge deler godtas.

Vær oppmerksom på at et stoppkriterium alltid treffer litt for sent: løkka oppdager først at $y$ er blitt negativ, altså i steget etter nedslaget. Trenger du nedslagspunktet nøyaktig, interpolerer du mellom de to siste punktene.

Regelen om at alle tall som mates inn i koden, skal være i SI-enheter — meter, kilogram, sekund — og at svaret dermed også kommer i SI.

Koden regner uten å vite noe om enheter. Blander du $\text{g}$ og $\text{kg}$, eller $\text{cm}$ og $\text{m}$, kommer det ingen feilmelding: du får et pent tall som er galt med en faktor tusen.

Det praktiske grepet er å skrive enheten som kommentar bak hver konstant når du setter opp koden. Det tar to sekunder og fanger den ene feilen ingen kontroll ellers oppdager.

Skillet mellom størrelser som må lagres for hvert tidspunkt, og størrelser som bare brukes og kastes.

Lagres (egen array): posisjon, fart, tid — alt du vil kunne plotte, lete i eller lese av til slutt.

Kastes (vanlig variabel): akselerasjonen, farten regnet ut av komponentene, avstanden til sentrum. De regnes på nytt i hvert gjennomløp og trengs ikke etterpå.

Denne inndelingen forklarer hvorfor a skrives uten indeks mens v alltid har indeks — og det er verdt å si i en besvarelse, slik at den manglende indeksen på a ikke leses som slurv.

Kravet om at du skriver hvorfor du valgte Euler–Cromer, ikke bare at du gjorde det.

Sensorveiledningene i H2023, V2024 og V2025 trekker dette fram som et skille mot toppkarakter. En fullgod begrunnelse har to ledd: hva Euler–Cromer gjør (oppdaterer farten først og bruker den oppdaterte farten videre, slik at feilen kansellerer over en periode), og hva alternativet gjør galt (Forward Euler bruker den gamle farten og tilfører energi systematisk).

To setninger. De koster et halvt minutt og er blant de billigste poengene i hele settet.

Å regne ut én egenskap ved bevegelsen på en helt annen måte enn koden gjør, og sammenlikne.

De tre beste kontrollene i mekanikk er:

Energi. Er systemet konservativt, skal $E = \tfrac12 mv^2 + U$ være konstant. Amplituder og maksfarter kan regnes ut av energibevaring uten å løse likningen.

Grensetilfelle. Slå av det vanskelige leddet og sammenlikn med den kjente løsningen.

Størrelsesorden. Er tallet fysisk rimelig? En pendel med sekundperiode, en ball som når tjue meter, en terminalfart på noen titalls meter per sekund.

Ingen av dem krever en maskin, og alle tre kan skrives inn i en håndskrevet besvarelse som en kontrollsetning.