7.3 ODE-systemer: pakke ut u og plotte komponentene
Den tunge sammensatte oppgaven: skrive høyresiden for et system, sette opp løseren og plotte komponentene — vitenskapelig innpakning som roterer.
Dette er kapitlet for den tyngste enkeltoppgaven på eksamen. Sjanger J — et ODE-system løst med det utleverte hierarkiet og plottet — står i hvert av de åtte gjennomgåtte settene (100 %) og er verdt 9–10 poeng, altså rundt en åttendedel av hele den avsluttende eksamenen.
Det som roterer, er innpakningen: epidemimodeller (SEIR, SIRD, SEID), rovdyr og byttedyr, Lorenz-likningene for vær, planetbaner, nervecellemodeller. Øvelsen er den samme hver gang, og den har tre trinn: skriv høyresiden, løs, pakk ut komponentene. Klarer du det på én modell, klarer du det på alle.
Oppgaven er nesten alltid flerdelt: (a) skriv høyresiden, (b) løs med hierarkiet, (c) plott komponentene med forklarende legende. Deloppgavene kan løses uavhengig — får du ikke (a) til, kan du fortsatt skrive (b) og (c) og få uttelling for dem.
Prioritet: dette må sitte — høyeste prioritet i hele boka sammen med kodesporing, klasser, rekker og fillesing.
Løkke 1 — hva et system er, og hvordan høyresiden skrives (~16 min) · Løkke 2 — å pakke ut komponentene (~14 min) · Løkke 3 — parametre og andreordens likninger (~16 min) · Løkke 4 — plottet (~14 min)
Forkunnskaper — sist du var her
1. De tre kallene fra kap. 7.2. De er nøyaktig de samme for et system; det eneste som endres, er at startverdien blir en liste:
import numpy as np
from ODESolver import RungeKutta4
def f(t, u):
return -0.5*u + np.sin(t)
solver = RungeKutta4(f)
solver.set_initial_condition(2)
t, u = solver.solve((0, 10), 100)
print("u.shape =", u.shape)Utskrift:
u.shape = (101,)2. Utpakking av flere verdier på én linje (kap. 1.2). Står det like mange navn på venstre side som verdier på høyre, fordeles de i rekkefølge:
verdier = [3, 7, 11]
a, b, c = verdier
print(c, b, a)Utskrift:
11 7 33. Klasse med __call__ (kap. 5.2). Et objekt med __call__ kan brukes som om det var en funksjon. Det er slik du lagrer parametre og likevel har noe løseren kan kalle:
class Lineaer:
def __init__(self, a, b):
self.a = a
self.b = b
def __call__(self, x):
return self.a*x + self.b
g = Lineaer(2, 5)
print(g(0), g(3))Utskrift:
5 114. Plotting med label og legend (kap. 4.3): flere kurver i samme figur, hver med sitt navn i forklaringsboksen.
Matematikken bak systemer — hvordan oppfører seg, og hva egenverdiene forteller — står i Differensiallikningssystem x′=Ax og diskret dynamikk. Du trenger den ikke for å kode, men den forklarer hvorfor løsningene svinger eller dør ut.
Løkke 1 — hva et system er, og hvordan høyresiden skrives
I en skog lever gauper og hare. Antall hare vokser av seg selv, men spises av gauper; antall gauper vokser bare hvis det finnes hare. Ingen av de to bestandene kan beskrives alene — de endrer seg sammen.
Med for hare og for gauper blir modellen to likninger som deler ukjente:
Det er et system av differensiallikninger. Det er fremdeles bare ett initialverdiproblem, men nå er tilstanden en vektor , og høyresiden returnerer en vektor av deriverte.
Og her er den gode nyheten: løseren du fikk utlevert, takler dette uten en eneste endring. Alt som skjer, er at u blir en liste av tall i stedet for ett tall.
Høyresiden får hele vektoren inn og returnerer alle de deriverte, i samme rekkefølge.
Rekkefølgen er ditt valg, men den må være den samme på alle tre stedene: når du pakker ut i f, når du returnerer de deriverte, og når du oppgir startverdiene. Velger du , må u0 være [x0, y0] og returverdien [dx, dy].
Ordet «vektor» betyr her bare «flere tall i rekkefølge» — i Python en liste eller et NumPy-array.
Tre trinn, alltid de samme. Dette er malen du skal kunne skrive fra bunnen.
1. Pakk ut, regn ut, returner liste.
def f(t, u):
x, y = u
dx = 0.8*x - 0.04*x*y
dy = 0.02*x*y - 0.5*y
return [dx, dy]Utskrift: ingen (bare høyresiden; hele programmet står i Eksempel 1).
2. Løs, med u0 som liste.
solver = RungeKutta4(f)
solver.set_initial_condition([40, 9])
t, u = solver.solve((0, 40), 800)Utskrift: ingen (fragment; krever f og importen over).
3. Pakk ut komponentene som kolonner.
x, y = u[:, 0], u[:, 1]Utskrift: ingen (én tilordning).
Tre detaljer avgjør poengene: u pakkes ut med like mange navn som likninger; f returnerer en liste, ikke et enkelt tall; og komponentene hentes som kolonner, u[:, k], ikke rader.
Løs systemet , med på (måneder) med RungeKutta4 og 800 steg. Hvilken form har løsningen, og hvor store blir bestandene?
import numpy as np
from ODESolver import RungeKutta4
a, b, c, d = 0.8, 0.04, 0.02, 0.5
def f(t, u):
x, y = u
dx = a*x - b*x*y
dy = c*x*y - d*y
return [dx, dy]
solver = RungeKutta4(f)
solver.set_initial_condition([40, 9])
t, u = solver.solve((0, 40), 800)
print("u.shape =", u.shape)
byttedyr = u[:, 0]
rovdyr = u[:, 1]
print(f"start: byttedyr = {byttedyr[0]:.1f} rovdyr = {rovdyr[0]:.1f}")
print(f"t = 40: byttedyr = {byttedyr[-1]:.2f} rovdyr = {rovdyr[-1]:.2f}")
print(f"største byttedyrbestand: {byttedyr.max():.2f}")
print(f"største rovdyrbestand: {rovdyr.max():.2f}")Utskrift:
u.shape = (801, 2)
start: byttedyr = 40.0 rovdyr = 9.0
t = 40: byttedyr = 11.86 rovdyr = 10.02
største byttedyrbestand: 60.13
største rovdyrbestand: 40.91u.shape er (801, 2) — 801 tidspunkter nedover, 2 komponenter bortover. Én rad per tidspunkt, én kolonne per ukjent. Det er hele forskjellen fra kap. 7.2, der u var 1-dimensjonal.
Legg merke til a, b, c, d = 0.8, 0.04, 0.02, 0.5 øverst: parametrene står ett sted, og høyresiden leser dem. Det er den enkleste måten, og den er helt akseptabel på eksamen. Skal parametrene kunne varieres uten å redigere f, bruker du klasse-varianten i Løkke 3.
Bestandene svinger: haren tar seg opp, gaupene følger etter, haren kollapser, gaupene sulter, og syklusen gjentar seg. Det er den klassiske oppførselen til denne modellen, og at du kan si det, er en del av svaret på en plotteoppgave.
(Innsteg — formen på løsningen.) Du løser et system med tre ukjente funksjoner på med 100 steg.
a) Hvilken form (shape) har u?
b) Hvor mange tall er det i u til sammen?
c) Hvordan henter du ut hele kurven for den andre komponenten?
Løkke 2 — å pakke ut komponentene
Her mistes flest poeng i denne sjangeren, og feilen er alltid den samme: å blande rader og kolonner.
Tenk på u som en tabell. Nedover går tiden, bortover går komponentene. Da er en hel kurve en kolonne, og en hel rad er «alle komponentene i ett øyeblikk».
Kolonne k av løsningsmatrisen: hele tidsforløpet for komponent nummer k. Kolon betyr «alle rader».
S, I, R = u[:, 0], u[:, 1], u[:, 2]Utskrift: ingen (én tilordning; forutsetter at u er løst).
Resultatet er 1-dimensjonalt og like langt som t — altså klart til å plottes mot t.
Til sammenligning: u[n, :] (eller bare u[n]) er én rad: alle komponentene ved tidspunkt n. Det er en vektor med like mange tall som systemet har likninger, og den bruker du nesten aldri etter at løsningen er ferdig.
Regelen: kolon foran kommaet gir en kurve, kolon bak kommaet gir et øyeblikk.
Bruk løsningen fra Eksempel 1 og skriv ut formen på u[:, 0] og på u[0, :], samt rad nummer 5.
import numpy as np
from ODESolver import RungeKutta4
def f(t, u):
x, y = u
return [0.8*x - 0.04*x*y, 0.02*x*y - 0.5*y]
solver = RungeKutta4(f)
solver.set_initial_condition([40, 9])
t, u = solver.solve((0, 40), 800)
print("u[:, 0] har form", u[:, 0].shape, "— hele byttedyrkurven")
print("u[0, :] har form", u[0, :].shape, "— begge bestandene ved t = 0:", u[0, :])
print("u[5] er en rad:", np.round(u[5], 3))Utskrift:
u[:, 0] har form (801,) — hele byttedyrkurven
u[0, :] har form (2,) — begge bestandene ved t = 0: [40. 9.]
u[5] er en rad: [44.481 9.81 ]Legg merke til at u[5] betyr u[5, :] — når du bare gir én indeks til et 2D-array, er det raden du får.
Den vanligste konsekvensen av feilen: du skriver plt.plot(t, u[0]) og får en feilmelding om at lengdene ikke passer, fordi t har 801 tall og u[0] bare 2. Ser du en slik lengdefeil i et plott, er det nesten alltid rad-i-stedet-for-kolonne.
u. u[0] er raden ved første tidspunkt, ikke den første komponentens kurve. Kurven er u[:, 0]. Symptomet er en lengdefeil i plottet, eller en kurve med bare to punkter.Returnere et tall i stedet for en liste. Glemmer du hakeparentesene og skriver return dx, dy, får du en tuppel — det virker faktisk, fordi modulen pakker resultatet i et array. Men skriver du return dx + dy, har du summert de deriverte, og løsningen blir tull uten feilmelding. Én linje, én komponent, i riktig rekkefølge.
Rekkefølgen glipper mellom de tre stedene. Pakker du ut S, I, R = u men oppgir u0 = [0, 10, 4990], har du byttet om S og R. Skriv rekkefølgen én gang og hold deg til den — mange skriver den som en kommentar over f, og det er en legitim bruk av kommentarer.
#1 — av-én på antall punkter. solve((0, T), N) gir N+1 rader. Det gjelder like fullt for systemer.
#14 — API-drift. Høyresiden er f(t, u), tid først. I eldre sett var det f(u, t). Feilen gir ingen feilmelding når t ikke brukes i høyresiden, som er nettopp tilfellet i de fleste systemene her. Les vedlegget.
Rekursjon. Systemene her er koblede, og det kan friste å tenke «komponent 1 kaller komponent 2». Ikke gjør det: rekursjon er utenfor pensum i dette emnet og markeres som uønsket i de offisielle løsningsforslagene. Alle komponentene regnes i samme kall til f, og framdriften skjer i løserens løkke.
der . Bruk , , og .
a) Skriv høyresiden f(t, u).
b) Løs på med RungeKutta4 og 1500 steg.
c) Skriv ut når smittetoppen kommer, hvor mange som dør, og hvor mange som aldri blir smittet.
Løkke 3 — parametre og andreordens likninger
— naturlig pausepunkt —
To utvidelser dekker resten av det oppgavene spør om. Den første er hvordan du lagrer parametre uten globale variabler. Den andre er hvordan en likning med dobbelderivert blir et system — det er slik planetbaner, fjærer og nerveceller kommer inn.
En spesialmetode (en metode med to understreker på hver side av navnet, som Python kaller automatisk) som gjør en instans kallbar: du kan sette parenteser bak objektet, og løseren merker ingen forskjell fra en vanlig funksjon.
Mønsteret er alltid det samme: parametrene inn i __init__, modellen i __call__.
class Fjaer:
def __init__(self, k, c):
self.k = k
self.c = c
def __call__(self, t, u):
x, v = u
return [v, -self.k*x - self.c*v]Utskrift: ingen (bare klassen; den brukes i Eksempel 3).
Deretter solver = RungeKutta4(Fjaer(4, 0.5)) — eller lagre instansen i en variabel først. Fordelen: du kan løse samme modell med flere parametersett uten å redigere koden, og hver parameter har et navn. Husk self foran hvert attributt inne i __call__; å glemme det er feilkode #3.
Første komponent i høyresiden er alltid v — det er definisjonen av farten. Andre komponent er akselerasjonen.
def f(t, u):
x, v = u
return [v, -x]Utskrift: ingen (høyresiden for ; brukes i Eksempel 3).
Dette er oppskriften bak alle mekanikk-innpakningene. En planetbane i planet har to romretninger og to farter, altså fire komponenter: .
En masse i en fjær med demping oppfyller . Løs med , , , , på med 2000 steg. Skriv høyresiden som en klasse.
import numpy as np
from ODESolver import RungeKutta4
class Fjaer:
def __init__(self, k, c):
self.k = k
self.c = c
def __call__(self, t, u):
x, v = u
return [v, -self.k*x - self.c*v]
modell = Fjaer(4, 0.5)
solver = RungeKutta4(modell)
solver.set_initial_condition([2, 0])
t, u = solver.solve((0, 20), 2000)
x = u[:, 0]
print(f"start: x = {x[0]:.3f}")
print(f"største utslag etter t = 5: {np.max(np.abs(x[t > 5])):.4f}")
print(f"x ved t = 20: {x[-1]:.5f}")Utskrift:
start: x = 2.000
største utslag etter t = 5: 0.5326
x ved t = 20: -0.00390Utslaget faller fra 2 til 0,53 i løpet av fem tidsenheter og er nesten borte ved — dempingen virker, som forventet.
To ting å legge merke til. x[t > 5] er boolsk indeksering (kap. 4.1): den plukker ut de elementene i x der tilsvarende t er større enn 5. Og modell sendes til løseren uten parenteser, akkurat som en funksjon — det er __call__ som gjør at det virker.
Skulle du løse for flere dempinger, er det nå trivielt: for c in [0, 0.5, 2]: solver = RungeKutta4(Fjaer(4, c)). Det er hele gevinsten med klasse-varianten.
Den udempede svingningen med , har den eksakte løsningen .
a) Skriv om til et system og skriv høyresiden f(t, u).
b) Løs på med RungeKutta4 og 1000 steg, og mål avviket fra for posisjonen og fra for farten.
c) Regn ut i første og siste punkt. Hva bør den være?
Løkke 4 — plottet
Nesten alle J-oppgaver har en siste deloppgave: «plott komponentene». Den er verdt et par poeng, den er rask, og den mistes ofte fordi legenden glemmes. Mønsteret er kort nok å pugge.
Flere kurver i samme figur, hver med sitt navn i forklaringsboksen:
plt.plot(t, S, label='mottakelige')
plt.plot(t, I, label='smittsomme')
plt.xlabel('dager')
plt.ylabel('antall personer')
plt.legend()
plt.show()Fire faste ledd: én plot-linje per kurve med label=, aksenavn, legend() for å vise navnene, og show() til slutt. Uten legend() blir label= aldri synlig — det er den vanligste tapte poengen i plottedelen.
Blokken over er merket text og ikke python fordi den er et fragment: den mangler både importen og løsningen den plotter. Det fullstendige programmet står i Eksempel 4.
Plott hare- og gaupebestanden fra Eksempel 1 i samme figur, med navn på aksene og en legende. Beskriv figuren i ord.
import matplotlib.pyplot as plt
from ODESolver import RungeKutta4
def f(t, u):
x, y = u
return [0.8*x - 0.04*x*y, 0.02*x*y - 0.5*y]
solver = RungeKutta4(f)
solver.set_initial_condition([40, 9])
t, u = solver.solve((0, 40), 800)
plt.plot(t, u[:, 0], label='byttedyr')
plt.plot(t, u[:, 1], label='rovdyr')
plt.xlabel('tid (måneder)')
plt.ylabel('bestand')
plt.legend()
plt.show()Plottekode får ingen «Utskrift»-blokk — den lager et vindu, ikke tekst. Til gjengjeld skal du kunne beskrive figuren, og det er nettopp det en eksamensbesvarelse med «kommenter figuren» ber om:
Figuren i ord. Vannrett akse er tid i måneder fra 0 til 40, loddrett akse er bestandsstørrelse fra rundt 0 til rundt 60. To kurver svinger opp og ned med samme periode på omtrent 12 måneder, men i utakt: byttedyrkurven (som starter i 40 og topper seg rundt 60) når toppen sin før rovdyrkurven (som starter i 9 og topper seg rundt 41). Etter hver byttedyrtopp faller byttedyrbestanden bratt mens rovdyrbestanden fortsetter å stige en stund, og deretter faller også den. Svingningene gjentar seg uten å dø ut. Legenden nede til høyre skiller de to kurvene.
Skriv plt.show() selv om du ikke kan kjøre koden på eksamen — den hører til mønsteret. Vil du plotte den ene bestanden mot den andre i stedet for mot tiden (en fasekurve), er det plt.plot(u[:, 0], u[:, 1]); da får du en lukket løkke som viser syklusen.
(Midtveisnivå, sjanger A — kodesporing, altså «hva skriver dette programmet ut?».)
import numpy as np
from ODESolver import ForwardEuler
def f(t, u):
a, b = u
return [b, -a]
solver = ForwardEuler(f)
solver.set_initial_condition([1, 0])
t, u = solver.solve((0, 1), 2)
print(u.shape)
print(np.round(u, 3))med . Bruk , , og .
a) Skriv høyresiden som en klasse med __init__ og __call__, slik at parametrene kan varieres.
b) Løs på med RungeKutta4 og 1800 steg, og pakk ut de fire komponentene.
c) Skriv ut smittetoppen med tidspunkt, sluttverdiene for og , og summen av alle gruppene til slutt.
d) Skriv plottekoden for de fire kurvene, og beskriv figuren i ord.
Bruk , startposisjon og startfart — som gir en sirkelbane.
a) Skriv om til et system med fire komponenter og skriv høyresiden.
b) Løs på med RungeKutta4 og 2000 steg.
c) Skriv ut avstanden fra sentrum i første og siste punkt, og posisjonen etter ett omløp. Hva bør du få?
Begrepsbank
Begrepsbanken er flashcard-/repetisjonsstoff — den gjentar det du nettopp har lest. Hopp trygt over ved førstegangslesing; tidsanslaget for kapitlet gjelder kjernestoffet.
Et 2-dimensjonalt array med form (N+1, m): rader er tidspunkter, kolonner er komponenter.
u[n] er hele tilstanden ved tidspunkt n; u[:, k] er hele kurven for komponent k. Skal du plotte, er det kolonner du trenger.
Rekkefølgen på komponentene velger du selv, men den må være den samme på tre steder: i utpakkingen S, I, R = u, i returlista [dS, dI, dR], og i u0.
Bytter du om på ett av stedene, får du et annet system uten noen feilmelding. Mange skriver rekkefølgen som en kommentar over høyresiden — det er en legitim bruk av kommentarer.
u[:, 0] er en kurve (like lang som t); u[0, :] er ett øyeblikk (like langt som antall likninger).Symptomet på forveksling: plottet klager over at lengdene ikke passer, eller kurven har bare to punkter. Kolon foran kommaet gir kurven.
Mange systemer har en størrelse som ikke skal endre seg: befolkningstotalen i en epidemimodell, energien i en svingning, avstanden fra sentrum i en sirkelbane.
Regn den ut i første og siste punkt. Holder den seg, er oppsettet nesten sikkert riktig — dette er den raskeste kontrollen du har.
En klasse med __init__ for parametrene og __call__(self, t, u) for modellen. Instansen sendes til løseren uten parenteser.
Fordelen er at du kan løse samme modell for flere parametersett uten å endre koden. Husk self. foran hvert attributt inne i __call__ (feilkode #3).
Første komponent i returlista er alltid v. Alle mekanikk-innpakningene bygger på dette, og en bane i planet gir fire komponenter.
I.argmax() gir indeksen til den største verdien, ikke verdien selv. Tidspunktet er t[I.argmax()].Å skrive ut indeksen i stedet for tidspunktet er en klassisk slurvefeil i epidemioppgaver: «toppen kom ved 610» når svaret var «ved dag 61».
label= gir en kurve navn, men navnet blir bare synlig hvis du kaller plt.legend().Er det flere kurver i samme figur, er legenden en del av oppgavesvaret — den forteller leseren hvilken kurve som er hvilken. Å glemme den er den vanligste tapte poengen i plottedelen.
Å beskrive et plott presist i tekst: hvilke akser, hvilket område, hvilke kurver, hvilken retning, hvor toppene er, og hva legenden sier.
Dette er en egen ferdighet oppgavene ber om («kommenter figuren»), og den kan trenes uten å kjøre koden.
Et plott av én komponent mot en annen i stedet for mot tiden: plt.plot(u[:, 0], u[:, 1]).
For rovdyr og byttedyr blir det en lukket løkke som viser syklusen; for en dempet svingning en spiral inn mot likevekten.
Parametre kan stå som variabler over høyresiden (beta = 0.35) eller lagres i en instans (self.beta).
Begge gir full uttelling. Den første er kortere og fin i en eksamensbesvarelse; den andre er nødvendig når oppgaven ber om å variere parametrene.
Plukker ut de elementene i x der betingelsen på t er sann. Arrayene må være like lange.
Brukes til å se på en del av løsningen — for eksempel det største utslaget etter at en innsvingingsfase er over.
NumPy skriver et 2D-array med doble hakeparenteser, én rad per linje, og kolonnene justert med mellomrom:
import numpy as np
print(np.array([[1.0, 0.0], [0.75, -1.0]]))Utskrift:
[[ 1. 0. ]
[ 0.75 -1. ]]Merk mellomrommet der minustegnet ellers ville stått — NumPy holder kolonnebredden.
I SIR-, SEIR- og SIRD-modellene forekommer hvert ledd to ganger med motsatt fortegn: det som forlater én gruppe, kommer inn i en annen.
Derfor er summen av alle de deriverte null, og totalen er konstant. Går totalen din, har du mistet eller duplisert et ledd.
Dette kapitlet er skrevet av Anthropics toppmodeller (Claude Opus og Claude Fable) og er foreløpig ikke manuelt gjennomgått — kvalitetskontrollen gjøres av uavhengige KI-agenter, og innmeldte feil rettes fortløpende. Funnet en feil? Meld fra, så retter vi den. Les mer om hvordan innholdet lages.
Skolesaga er en uavhengig læringsressurs og er ikke tilknyttet eller godkjent av Universitetet i Oslo. Dette er ikke offisielt studiemateriell. Les mer.