7.2 Det utleverte ODESolver-hierarkiet: bruke ForwardEuler og RungeKutta4

Bruke den vedlagte `ODESolver`-basisklassen med subklassene `ForwardEuler` og `RungeKutta4` på en skalar ODE — instansiering, initialbetingelse og `solve`.

55 min

7 oppgaver

Det utleverte ODESolver-hierarkietbruke ForwardEulerRungeKutta4

Din fremgang i kapitlet

0 / 7 oppgaver

Kapitlets plass i kurset

Bygger på

7.1ODE-løser som funksjon: forward Euler 5.3Arv og `super()`

Brukes videre i

7.3ODE-systemer: pakke ut u og plotte komponentene

Dette kapitlet er inngangsbilletten til den tyngste oppgaven på hele eksamen. Sjanger J — «løs et ODE-system med det utleverte løser-hierarkiet» — står i hvert av de åtte gjennomgåtte settene (100 %) og er verdt 9–10 poeng. Selve systemdelen kommer i kap. 7.3; her lærer du bruksmønsteret på en skalar likning, altså en likning med bare én ukjent funksjon.

Det du får utlevert, er et vedlegg: en PDF med en ferdigskrevet Python-modul som inneholder klassen ODESolver og metodene ForwardEuler og RungeKutta4. Du skal ikke skrive den. Du skal kunne lese den og bruke den, og det er tre kall.

Den gode nyheten: bruksmønsteret har vært identisk hvert år. Det eneste som har endret seg, er navnene og rekkefølgen på argumentene — det som i boka kalles API-drift, og som har eget nummer i feilkatalogen (#14). Derfor er regelen: les vedlegget først, så skriv koden.

Prioritet: dette må sitte. Del 7 er ett av de fem temaene med 100 % frekvens.

Løkke 1 — de tre kallene (~14 min) · Løkke 2 — hva som står i vedlegget (~16 min) · Løkke 3 — hvilken metode skal du velge (~13 min) · Løkke 4 — API-driften og fellene (~12 min)

Forkunnskaper — sist du var her

1. Løser-løkka fra kap. 7.1. Hele hierarkiet i dette kapitlet er denne løkka pakket inn i en klasse:

import numpy as np

def forward_euler(f, T, N, u0):
    t = np.linspace(0, T, N+1)
    u = np.zeros_like(t)
    u[0] = u0
    dt = t[1] - t[0]
    for n in range(N):
        u[n+1] = u[n] + dt*f(t[n], u[n])
    return t, u

def f(t, u):
    return -2*u

t, u = forward_euler(f, 3, 30, 1)
print(f"{len(t)} punkter, u(3) = {u[-1]:.6f}")

Utskrift:

31 punkter, u(3) = 0.001238

2. Klasser med __init__ og metoder (kap. 5.1). En instans er et objekt laget av en klasse; self er instansen selv, og attributtene lagres på den:

class Teller:
    def __init__(self, start):
        self.verdi = start

    def opp(self):
        self.verdi = self.verdi + 1
        return self.verdi

t = Teller(10)
print(t.opp(), t.opp(), t.verdi)

Utskrift:

11 12 12

3. Arv (kap. 5.3). En subklasse arver alt fra basisklassen og kan legge til eller overskrive enkeltmetoder:

class Grunnform:
    def hilsen(self):
        return "jeg er " + self.navn()

    def navn(self):
        return "ukjent"

class Spesial(Grunnform):
    def navn(self):
        return "spesial"

print(Spesial().hilsen())

Utskrift:

jeg er spesial

Det siste er nøkkelen til hele hierarkiet: ODESolver er grunnformen, og ForwardEuler og RungeKutta4 overskriver én metode hver.

Løkke 1 — de tre kallene

Du får utlevert en modul som allerede kan løse ODE-er. Da er jobben din redusert til fire ting: skriv høyresiden, velg metode, oppgi startverdien, og løs.

solver = RungeKutta4(f)              # 1. velg metode, send høyresiden
solver.set_initial_condition(2)      # 2. oppgi startverdien
t, u = solver.solve((0, 10), 100)    # 3. løs på [0, 10] med 100 steg

Utskrift: ingen (fragmentet mangler f og importen; det fullstendige programmet står i Eksempel 1).

Det er alt. Legg merke til at solve tar to argumenter: en tuppel med start- og sluttid, og antall steg. Og at den returnerer to arrays, i rekkefølgen tid først.

`solver = RungeKutta4(f)`

Lager en instans — et objekt — av løserklassen, og gir den høyresiden. Ingen regning skjer her; løseren bare husker f . Vil du bytte metode, bytter du bare klassenavnet : solver = ForwardEuler(f) Høyresiden f sendes uten parenteser, fordi løseren skal kalle den selv i hvert steg. Det er samme prinsipp som når du sendte f til forward_euler i kap. 7.1 — forskjellen er bare at den nå lagres på objektet som self.f .

`solver.set_initial_condition(u0)`

Oppgir startverdien. Må kalles før solve — glemmer du den, vet ikke løseren hvor den skal begynne, og programmet krasjer. For en skalar likning er u0 ett tall: solver.set_initial_condition(2) For et system er u0 en liste med én startverdi per komponent ([9990, 10, 0, 0]) — se kap. 7.3 . Det er nettopp i denne metoden modulen finner ut om problemet er skalart eller et system, ved å se på hva slags verdi du ga den.

`t, u = solver.solve(t_span, N)`

Gjør regnearbeidet og returnerer to arrays: tidspunktene og løsningen. - t_span er en tuppel (t0, T) med start- og sluttid — legg merke til dobbeltparentesen i solve((0, 10), 100) : den ytre hører til metodekallet, den indre er tuppelen. - N er antall steg. Du får N+1 punkter tilbake, som alltid. t, u = solver.solve((0, 10), 100) print(len(t)) Utskrift: 101 Dette er post-2023-signaturen, og den boka bruker overalt. Den eldre varianten tok i stedet et ferdig array av tidspunkter — se advarselen om API-drift nedenfor.

✏️Eksempel 1: en skalar likning med RungeKutta4

En sparekonto tappes jevnt samtidig som den får renter i takt med årstidene: $u' = -0{,}5u + \sin t$ , $u(0) = 2$ . Løs på $[0, 10]$ med RungeKutta4 og 100 steg. Hvor mange punkter får du, og hvilken form har løsningen?

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("antall punkter:", len(t))
print("u.shape =", u.shape)
print(f"u(0)  = {u[0]:.6f}")
print(f"u(10) = {u[-1]:.6f}")

Utskrift:

antall punkter: 101
u.shape = (101,)
u(0)  = 2.000000
u(10) = 0.472515

101 punkter — N+1, som alltid. For en skalar likning er u 1-dimensjonal: ett tall per tidspunkt. u.shape er (101,), ikke (101, 1). Det er verdt å merke seg, for i kap. 7.3 blir u 2-dimensjonal så snart du løser et system, og da endres måten du henter ut verdier.

Merk from ODESolver import RungeKutta4 på linje 2. På eksamen får du modulen som vedlegg, og du kan anta at den er tilgjengelig for import. Skriv importlinja likevel — den koster ingenting og viser at du vet hvor klassen kommer fra. Filnavnet kan variere fra år til år; står det et annet navn i vedlegget, bruker du det.

📝Oppgave 1

(Innsteg — bruksmønsteret.) Under står et program der én linje er fjernet. Programmet skal løse $u' = -0{,}5u$ , $u(0) = 10$ , på $[0, 4]$ med 8 steg.

from ODESolver import ForwardEuler

def f(t, u):
    return -0.5*u

solver = ForwardEuler(f)

t, u = solver.solve((0, 4), 8)
print("t =", t)
print("u =", u.round(4))

a) Hvilken linje mangler, og hvor skal den stå?
b) Hvor mange tall står det i u?

Løkke 2 — hva som står i vedlegget

Du trenger ikke kunne skrive modulen, men du må kunne lese den, av to grunner: du må sjekke argumentrekkefølgen, og i kap. 7.4 skal du utvide den med en egen metode. Da må du vite hvilke attributter du har tilgang til.

Strukturen er alltid den samme. Basisklassen gjør alt det kjedelige — lager arrays, setter startverdien, kjører løkka — og lar én metode stå åpen: advance, som gjør ett steg. Hver metode er en subklasse som fyller inn nettopp den.

`ForwardEuler` og `RungeKutta4`

u_{n+1} = u_n + \Delta t f(t_n, u_n)

`advance()`

Metoden som gjør ett tidssteg og returnerer den neste verdien av u . Basisklassen implementerer den ikke — den kaster NotImplementedError — og hver subklasse fyller den inn. Det er hele arvemønsteret: basisklassen kjenner løkka, subklassen kjenner formelen. Derfor er svaret på «skriv en ny ODE-metode» alltid: skriv bare advance . Se kap. 7.4 .

Attributtene `self.f`, `self.u`, `self.t`, `self.n`, `self.dt`

Det basisklassen stiller til rådighet inne i advance : Attributt Hva det er self.f høyresiden, kallbar som self.f(t, u) self.u hele løsningsarrayet så langt self.t hele tidsarrayet self.n indeksen for steget som skal gjøres nå self.dt steglengden Med dem er forlengs Euler-steget én linje: return self.u[self.n] + self.dt*self.f(self.t[self.n], self.u[self.n]) . Navnene kan variere mellom årganger. Trenger du steglengden og er usikker, bruk den robuste formen dt = self.t[self.n+1] - self.t[self.n] — den virker uansett om modulen har et self.dt eller ikke.

Attributt	Hva det er
`self.f`	høyresiden, kallbar som `self.f(t, u)`
`self.u`	hele løsningsarrayet så langt
`self.t`	hele tidsarrayet
`self.n`	indeksen for steget som skal gjøres nå
`self.dt`	steglengden

Den utleverte ODESolver-modulen

✏️Eksempel 2: ForwardEuler mot RungeKutta4

Løs $u' = -2u$ , $u(0) = 1$ , på $[0, 3]$ med begge metodene og $N = 30, 60, 120, 240$ . Mål det største avviket fra den eksakte løsningen $e^{-2t}$ , og finn ordenen til hver metode.

import numpy as np
from ODESolver import ForwardEuler, RungeKutta4

def f(t, u):
    return -2*u

print("     N   ForwardEuler   RungeKutta4")
for N in [30, 60, 120, 240]:
    rad = [N]
    for metode in [ForwardEuler, RungeKutta4]:
        solver = metode(f)
        solver.set_initial_condition(1)
        t, u = solver.solve((0, 3), N)
        rad.append(np.max(np.abs(u - np.exp(-2*t))))
    print(f"{rad[0]:6d}   {rad[1]:12.4e}  {rad[2]:12.4e}")

Utskrift:

     N   ForwardEuler   RungeKutta4
    30     4.0199e-02    5.7970e-06
    60     1.9201e-02    3.3324e-07
   120     9.3935e-03    1.9976e-08
   240     4.6470e-03    1.2227e-09

ForwardEuler-kolonnen halveres nedover (forhold 2,09 · 2,04 · 2,02): orden 1. RungeKutta4-kolonnen deles på omtrent 16 (17,4 · 16,7 · 16,3): orden 4.

Se på tallene på samme linje: med 30 steg er RungeKutta4 nesten 7 000 ganger nøyaktigere enn ForwardEuler, for fire funksjonskall i stedet for ett. Det er derfor RungeKutta4 er standardvalget når oppgaven ikke sier noe annet.

Merk også at kolonnen for ForwardEuler er de samme tallene som den egenskrevne forward_euler ga i kap. 7.1. Det er en god kontroll på at du har forstått at hierarkiet gjør nøyaktig det samme — bare pakket inn i klasser.

Typiske feil ved bruk av hierarkiet

#14 — API-forveksling. Den store. Skriver du høyresiden som f(u, t) mot en modul som kaller f(t, u), blir tid og tilstand byttet om. Det kompilerer og gir tall — bare helt gale tall. Samme problem med solve: den nye tar (t_span, N), den gamle et array av tidspunkter. Les vedlegget først.

Glemme set_initial_condition. Da mangler startverdien, og programmet stopper. Med bokas implementasjon får du siste linje:

AttributeError: 'ForwardEuler' object has no attribute 'neq'

Den eksakte meldingen avhenger av hvordan modulen er skrevet, men mønsteret er alltid det samme: en AttributeError fra solve fordi noe som skulle vært satt, ikke er satt. Ser du den, er det nesten alltid den glemte startverdien.

Sende f med parenteser. RungeKutta4(f(t, u)) sender ett tall inn i stedet for funksjonen. Skriv RungeKutta4(f).

Glemme dobbeltparentesen. solver.solve(0, 10, 100) er tre argumenter og feiler; det skal være solver.solve((0, 10), 100) — tuppelen er ett argument.

Rekursjon. Det er fristende å tenke at advance kan kalle seg selv for å «gå videre». Den skal ikke det: rekursjon er utenfor pensum i dette emnet og markeres som uønsket i de offisielle løsningsforslagene. Basisklassens løkke gjør framdriften; advance gjør bare ett steg.

📝Oppgave 2

Eksamensnivå, sjanger J

En fiskebestand følger logistisk vekst, $u' = 0{,}3u\left(1 - \dfrac{u}{100}\right)$ , med $u(0) = 5$ (målt i tusen tonn).

a) Løs likningen på $[0, 50]$ med RungeKutta4 og 200 steg.
b) Skriv ut bestanden ved $t = 50$ .
c) Finn det første tidspunktet der bestanden passerer 50.

📝Oppgave 3

Eksamensnivå, sjanger J

Bruk samme fiskebestand som i oppgave 2, men på $[0, 20]$ . Den eksakte løsningen er ikke oppgitt, så bruk RungeKutta4 med 4000 steg som referanse.

a) Skriv en hjelpefunksjon som løser problemet med en gitt metode og et gitt antall steg og returnerer sluttverdien.
b) Skriv ut avviket fra referansen for ForwardEuler og RungeKutta4 med $N = 5, 10, 20, 40$ .
c) Hvilken metode trenger færrest steg for tre riktige desimaler?

Løkke 3 — å kontrollere at svaret er riktig

— naturlig pausepunkt —

En numerisk løsning er alltid feil. Spørsmålet er om den er nok riktig, og en besvarelse som viser at du kan sjekke det, skiller seg positivt ut. Du har tre kontroller, og alle tre er raske:

1. Sammenlign med en eksakt løsning der du har en.
2. Forfin steget og se om svaret slutter å endre seg.
3. Sjekk fornuften: er høyresiden null i startpunktet, skal løsningen stå stille; en bestand skal ikke bli negativ; en temperatur skal ikke gå under romtemperaturen.

Likningen fra Eksempel 1 kan faktisk løses eksakt, og det gir den beste kontrollen.

Runge–Kutta 4 (RK4)

O(\Delta t^4)

✏️Eksempel 3: kontroll mot den eksakte løsningen

Likningen $u' = -0{,}5u + \sin t$ med $u(0) = 2$ har den eksakte løsningen $u(t) = 2{,}8e^{-0{,}5t} + 0{,}4\sin t - 0{,}8\cos t$ . Kontroller den numeriske løsningen fra Eksempel 1.

import numpy as np
from ODESolver import RungeKutta4

def f(t, u):
    return -0.5*u + np.sin(t)

def eksakt(t):
    return 2.8*np.exp(-0.5*t) + 0.4*np.sin(t) - 0.8*np.cos(t)

solver = RungeKutta4(f)
solver.set_initial_condition(2)
t, u = solver.solve((0, 10), 100)
print(f"maks avvik fra den eksakte løsningen: {np.max(np.abs(u - eksakt(t))):.2e}")

Utskrift:

maks avvik fra den eksakte løsningen: 7.43e-08

Sju siffer riktig med 100 steg. Merk hvordan eksakt(t) kalles med hele tidsarrayet og gir hele kurven tilbake — det er vektorisering (kap. 4.2), og det er grunnen til at kontrollen blir én linje.

Hvis du vil se hvor den eksakte løsningen kommer fra, er det den lineære førsteordens likningen fra Differensialligninger: separable og lineære første ordens. På eksamen får du den oppgitt hvis du skal bruke den.

📝Oppgave 4

(Midtveisnivå, sjanger A — kodesporing, altså «hva skriver dette programmet ut?».) Modulen ODESolver er den fra dette kapitlet.

from ODESolver import RungeKutta4

def f(t, u):
    return -u

solver = RungeKutta4(f)
solver.set_initial_condition(1)
t, u = solver.solve((0, 2), 4)
print(len(t), u.shape, t[1] - t[0])

📝Oppgave 5

Eksamensnivå, sjanger…

Et program fra et eldre oppgavesett er skrevet mot den GAMLE modulen:

def f(u, t):
    return -0.2*u

solver = ForwardEuler(f)
solver.set_initial_condition(3)
tidspunkter = np.linspace(0, 10, 51)
u, t = solver.solve(tidspunkter)

a) Nevn tre ting som må endres for at programmet skal virke mot modulen i dette kapitlet.
b) Skriv om programmet, og skriv ut sluttverdien.

📝Oppgave 6

Eksamensnivå, sjanger J…

RungeKutta4 gjengir polynomløsninger av lav grad eksakt. Det gir en presis test av at du bruker hierarkiet riktig.

a) Skriv en testfunksjon test_rk4() uten argumenter som løser $u' = 3t^2$ , $u(0) = 0$ , på $[0, 2]$ med 4 steg og sjekker sluttverdien mot den eksakte $u(2) = 2^3 = 8$ , med toleranse $10^{-10}$ .
b) Kall den og skriv ut en bekreftelse.
c) Forklar hvorfor toleransen kan være så stram.

📝Oppgave 7

Eksamensnivå, sjanger J

Kaffekoppen fra kap. 7.1 igjen: $u' = -0{,}05(u - 20)$ , $u(0) = 80$ , på $[0, 60]$ minutter.

a) Løs med RungeKutta4 og 600 steg.
b) Finn første tidspunkt under 60 grader og temperaturen etter én time.
c) Sammenlign med svaret din egen forward_euler ga i kap. 7.1 (8,1 minutter og 22,96 grader). Hvilket svar er nærmest sannheten, og hvorfor?

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.

Utlevert vedlegg (ODE-modulen)

En PDF med en ferdigskrevet Python-modul som følger ODE-oppgavene på eksamen. Den inneholder ODESolver og subklassene, og du skal bruke den, ikke skrive den. Det ene du må gjøre med vedlegget: sjekke argumentrekkefølgen i høyresiden og hva solve tar imot. Det er de to tingene som har endret seg mellom årganger.

`t_span`-tuppelen

Start- og sluttid samlet i én verdi: (0, 10) . Fordi det er ett argument, blir kallet solver.solve((0, 10), 100) med to parenteser. Skriver du solve(0, 10, 100), sender du tre argumenter og får TypeError .

Instans

Et konkret objekt laget fra en klasse. solver = RungeKutta4(f) lager én instans; du kan lage flere med forskjellige høyresider samtidig. Attributtene (self.u, self.t, \dots) tilhører instansen, så to løsere tråkker ikke i hverandres data.

Kallbar (callable)

Noe du kan sette parenteser bak. Både en vanlig funksjon og et objekt med__call__er kallbart. Derfor kan høyresiden være enten def f(t, u): eller en instans av en klasse med__call__(self, t, u) . Løseren merker ingen forskjell — den skriver bare self.f(t, u) . Klasse-varianten brukes når modellen har parametre, se kap. 7.3 .

`NotImplementedError`

En exception — altså en feil som avbryter programmet — som basisklassens advance kaster. Meningen er å si « ODESolver er ikke en løser i seg selv; bruk en subklasse». Får du den, har du instansiert basisklassen ved et uhell, eller glemt å skrive advance i din egen subklasse.

Skalar likning gir 1D-løsning

\to

Bytte metode = bytte klassenavn

Hele poenget med arv i dette hierarkiet: ForwardEuler(f) og RungeKutta4(f) brukes helt likt. Det gjør det trivielt å sammenligne metoder, og det er derfor du kan sende klassenavnet som argument til en hjelpefunksjon: def los(metode, N): solver = metode(f) .

Kostnad per steg

Antall kall til høyresiden i hvert tidssteg: ForwardEuler ett, Heun/midtpunkt to, RungeKutta4 fire. Regnestykket som avgjør valget: RK4 med 30 steg koster 120 kall og var i Eksempel 2 nesten 7 000 ganger nøyaktigere enn forlengs Euler med 30 steg. Høyere orden er nesten alltid billigere enn flere steg.

API-drift (#14)

At det utleverte grensesnittet har endret seg over tid: f(u, t) og solve(time_points) i 2020-2022, f(t, u) og solve(t_span, N) fra 2023. Feilen gir ingen feilmelding når t ikke brukes i høyresiden — den gir bare gale tall. Derfor: les vedlegget, og vær konsekvent.

`AttributeError` fra `solve`

Symptomet på glemt set_initial_condition : solve prøver å bruke noe som skulle vært satt. AttributeError: 'ForwardEuler' object has no attribute 'neq' Den nøyaktige meldingen avhenger av modulen. Mønsteret — en AttributeError idet du kaller solve — er alltid det samme.

Referanseløsning

En numerisk løsning med veldig mange steg, brukt som «sannhet» når det ikke finnes en eksakt løsning. Regel: bruk den mest nøyaktige metoden du har, med et par tusen steg, og mål de grove løsningene mot den. Sammenlign aldri en metode med seg selv på samme steglengde.

Asymptotisk orden

At ordenen gjelder når steget blir lite . Med svært få steg kan forholdene mellom feilene være ganske ulike 2 eller 16. Praktisk konsekvens: mål alltid over minst tre-fire steglengder, og se om forholdene nærmer seg det teoretiske tallet nedover i tabellen.

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.