Tilbake
7.3

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.

60 min
6 oppgaver
ODE-systemerpakke ut uplotte komponentene
Din fremgang i kapitlet
0 / 6 oppgaver

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 3

3. 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 11

4. Plotting med label og legend (kap. 4.3): flere kurver i samme figur, hver med sitt navn i forklaringsboksen.

Matematikken bak systemer — hvordan x=Axx' = Ax 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 xx for hare og yy for gauper blir modellen to likninger som deler ukjente:

x=axbxy,y=cxydyx' = ax - bxy, \qquad y' = cxy - dy

Det er et system av differensiallikninger. Det er fremdeles bare ett initialverdiproblem, men nå er tilstanden en vektor u=[x,y]u = [x, y], 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.

Tilstandsvektoren uu i et system
For et system med mm likninger samles de ukjente i én vektor:

u=[u0,u1,,um1],u=f(t,u)u = [u_0, u_1, \ldots, u_{m-1}], \qquad u' = f(t, u)

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 [x,y][x, y], 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.

📜System-ODE-malen

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.

✏️Eksempel 1: hare og gaupe (rovdyr og byttedyr)

Løs systemet x=0,8x0,04xyx' = 0{,}8x - 0{,}04xy, y=0,02xy0,5yy' = 0{,}02xy - 0{,}5y med u0=[40,9]u_0 = [40, 9][0,40][0, 40] (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.91

u.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.

📝Oppgave 1

(Innsteg — formen på løsningen.) Du løser et system med tre ukjente funksjoner på [0,20][0, 20] 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».

`u[:, k]` — komponent nummer kk

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.

✏️Eksempel 2: kolonne eller rad?

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.

📝Oppgave 2
Eksamensnivå, sjanger J
SIRD-modellen deler en befolkning i mottakelige SS, smittsomme II, friske RR og døde DD:

S=βSI/N,I=βSI/NγIμI,R=γI,D=μIS' = -\beta S I / N, \quad I' = \beta S I / N - \gamma I - \mu I, \quad R' = \gamma I, \quad D' = \mu I

der N=S+I+R+DN = S + I + R + D. Bruk β=0,35\beta = 0{,}35, γ=0,1\gamma = 0{,}1, μ=0,01\mu = 0{,}01 og u0=[4990,10,0,0]u_0 = [4990, 10, 0, 0].

a) Skriv høyresiden f(t, u).
b) Løs på [0,150][0, 150] 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.

Høyresiden som klasse med `__call__(self, t, u)`

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.

Andreordens likning skrevet om til et system
En likning med dobbelderivert, x=g(t,x,x)x'' = g(t, x, x'), gjøres om til to førsteordens likninger ved å innføre farten som egen ukjent:

u=[x,v],x=v,v=g(t,x,v)u = [x, v], \qquad x' = v, \qquad v' = g(t, x, v)

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 x=xx'' = -x; brukes i Eksempel 3).

Dette er oppskriften bak alle mekanikk-innpakningene. En planetbane i planet har to romretninger og to farter, altså fire komponenter: u=[x,y,vx,vy]u = [x, y, v_x, v_y].

✏️Eksempel 3: dempet fjær, med parametre i en klasse

En masse i en fjær med demping oppfyller x=kxcxx'' = -kx - cx'. Løs med k=4k = 4, c=0,5c = 0{,}5, x(0)=2x(0) = 2, x(0)=0x'(0) = 0, på [0,20][0, 20] 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.00390

Utslaget faller fra 2 til 0,53 i løpet av fem tidsenheter og er nesten borte ved t=20t = 20 — 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.

📝Oppgave 3
Eksamensnivå, sjanger J

Den udempede svingningen x=xx'' = -x med x(0)=1x(0) = 1, x(0)=0x'(0) = 0 har den eksakte løsningen x(t)=costx(t) = \cos t.

a) Skriv om til et system og skriv høyresiden f(t, u).
b) Løs på [0,10][0, 10] med RungeKutta4 og 1000 steg, og mål avviket fra cost\cos t for posisjonen og fra sint-\sin t for farten.
c) Regn ut x2+v2x^2 + v^2 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.

`plt.plot(t, kurve, label='navn')` og `plt.legend()`

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.

✏️Eksempel 4: plotte begge bestandene

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.

📝Oppgave 4

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

📝Oppgave 5
Eksamensnivå, sjanger J
SEIR-modellen legger til en gruppe EE av smittede som ennå ikke er smittsomme:

S=βSI/N,E=βSI/NσE,I=σEγI,R=γIS' = -\beta S I / N, \quad E' = \beta S I / N - \sigma E, \quad I' = \sigma E - \gamma I, \quad R' = \gamma I

med N=S+E+I+RN = S + E + I + R. Bruk β=0,4\beta = 0{,}4, σ=0,2\sigma = 0{,}2, γ=0,1\gamma = 0{,}1 og u0=[9990,10,0,0]u_0 = [9990, 10, 0, 0].

a) Skriv høyresiden som en klasse med __init__ og __call__, slik at parametrene kan varieres.
b) Løs på [0,180][0, 180] med RungeKutta4 og 1800 steg, og pakk ut de fire komponentene.
c) Skriv ut smittetoppen med tidspunkt, sluttverdiene for SS og RR, og summen av alle gruppene til slutt.
d) Skriv plottekoden for de fire kurvene, og beskriv figuren i ord.

📝Oppgave 6
Eksamensnivå, sjanger J
En satellitt i bane rundt en planet i planet oppfyller

x=GMxr3,y=GMyr3,r=x2+y2x'' = -\frac{GMx}{r^3}, \qquad y'' = -\frac{GMy}{r^3}, \qquad r = \sqrt{x^2 + y^2}

Bruk GM=1GM = 1, startposisjon (1,0)(1, 0) og startfart (0,1)(0, 1) — som gir en sirkelbane.

a) Skriv om til et system med fire komponenter og skriv høyresiden.
b) Løs på [0,2π][0, 2\pi] 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å?

API- og konstruksjonsliste

Begrepsbank

Begrepsbanken er flashcard-/repetisjonsstoff — den gjentar det du nettopp har lest. Hopp trygt over ved førstegangslesing; tidsanslaget for kapitlet gjelder kjernestoffet.

Løsningsmatrisen `u` for et system

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ølgeregelen

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.

Kolonne mot rad
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.

Bevart størrelse som fornuftssjekk

Mange systemer har en størrelse som ikke skal endre seg: befolkningstotalen i en epidemimodell, energien x2+v2x^2 + v^2 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.

Kallbar høyreside med parametre

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

Reduksjon av andreordens likning
x=g(t,x,x)x'' = g(t, x, x') blir et system ved å innføre farten: u=[x,v]u = [x, v], x=vx' = v, v=gv' = g.

Første komponent i returlista er alltid v. Alle mekanikk-innpakningene bygger på dette, og en bane i planet gir fire komponenter.

`argmax` og tidspunktet for en topp
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».

Legenden
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.

Figur-i-ord

Å 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.

Fasekurve

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.

Global parameter mot lagret parameter

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.

Boolsk indeksering `x[t > 5]`

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.

Utskriftsformen for et 2D-array

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.

Sum-kansellering i en kompartmentmodell

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.

Repetisjon

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.