7.4 Subklasse av ODESolver: skriv `advance`

Arv fra basisklassen i løser-hierarkiet. Klassen på innsiden av parentesen kalles basisklassen, og den nye klassen arver alle metodene dens: set_initial_condition, solve og oppstartsmetoden __init__.

Subklassen trenger derfor ingen __init__ med mindre den skal lagre noe eget, og den skal aldri skrive sin egen solve. Bruken er identisk med de utleverte metodene:

from ODESolver import ODESolver

class Metode(ODESolver):
    def advance(self):
        u, f, t, n, dt = self.u, self.f, self.t, self.n, self.dt
        return u[n] + dt*f(t[n], u[n])

def f(t, u):
    return -u

solver = Metode(f)
solver.set_initial_condition(1)
t, u = solver.solve((0, 1), 4)
print(f"{u[-1]:.6f}")

Utskrift:

0.316406

Merk hvor lite som står i klassen: én metode, to linjer. Alt det andre kom med arven.

De fem attributtene basisklassen stiller til rådighet. Et attributt er en verdi som er lagret på objektet selv.

Den robuste måten å finne steglengden på. Den regner differansen mellom to nabotidspunkter i tidsarrayet, i stedet for å stole på at modulen har lagret et attributt som heter dt.

Grunnen til å kunne den: attributtnavnene i modulen du får utlevert kan variere fra år til år. Tidsarrayet heter derimot alltid self.t, for det er halvparten av det solve returnerer.

from ODESolver import ODESolver

class EulerRobust(ODESolver):
    def advance(self):
        u, f, t, n = self.u, self.f, self.t, self.n
        dt = t[n+1] - t[n]
        return u[n] + dt*f(t[n], u[n])

def f(t, u):
    return -u

solver = EulerRobust(f)
solver.set_initial_condition(1)
t, u = solver.solve((0, 2), 8)
print(f"{u[-1]:.6f}")

Utskrift:

0.100113

Indeksen n+1 er alltid trygg her: advance kalles bare for n opp til N-1, og tidsarrayet har N+1 elementer.

Kallet som kjører basisklassens oppstartsmetode fra subklassens egen. Det skal stå som første linje i subklassens __init__, før du setter dine egne attributter.

from ODESolver import ODESolver

class EulerMedTelling(ODESolver):
    def __init__(self, f):
        super().__init__(f)
        self.antall_kall = 0

    def advance(self):
        u, f, t, n, dt = self.u, self.f, self.t, self.n, self.dt
        self.antall_kall = self.antall_kall + 1
        return u[n] + dt*f(t[n], u[n])

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

solver = EulerMedTelling(f)
solver.set_initial_condition(4)
t, u = solver.solve((0, 6), 60)
print(f"u(6) = {u[-1]:.5f}")
print("kall til hoyresiden:", solver.antall_kall)

Utskrift:

u(6) = 0.18428
kall til hoyresiden: 60

Regelen er den samme som i kap. 5.3: super().__init__(...) først, egne attributter etterpå. Trenger subklassen ingen egne attributter, skriver du ingen __init__ i det hele tatt — den arves.

Tallet 60 er verdt en tanke: forlengs Euler kaller høyresiden én gang per steg. Runge–Kutta 4 ville brukt 240 kall på samme antall steg. Det er den prisen du betaler for høyere orden.

$$g(v) = v - u_n - \Delta t\, f(t_{n+1}, v) = 0.$$

Nullpunktet $v$ er den nye verdien. Alt annet — $u_n$, $t_{n+1}$ og $\Delta t$ — er kjent når steget starter, og fungerer som konstanter i likningen.

Formen er den samme for systemer: da er $v$ en vektor, $g$ returnerer en vektor like lang, og du løser like mange likninger som du har ukjente funksjoner. Du skriver ikke noe ekstra for det tilfellet.

Startgjetningen er nesten alltid u[n], altså den forrige verdien. Den er per definisjon nær svaret når steget er lite, og det er akkurat det en nullpunktsmetode trenger for å konvergere raskt.

Nullpunktsløseren fra kap. 6.1, brukt én gang per tidssteg. Mønsteret er fast:

from scipy.optimize import root
from ODESolver import ODESolver

class BaklengsEuler(ODESolver):
    def advance(self):
        u, f, t, n, dt = self.u, self.f, self.t, self.n, self.dt

        def likning(v):
            return v - u[n] - dt*f(t[n+1], v)

        sol = root(likning, u[n])
        return sol.x

Utskrift: ingen (bare klassedefinisjonen; den brukes i Eksempel 4).

Tre detaljer avgjør om det virker:

- likning defineres inne i advance. Den må se u[n], t[n+1] og dt, og de er forskjellige for hvert steg. En funksjon definert inne i en annen ser den ytre funksjonens variabler.
- Tiden er t[n+1], ikke t[n]. Det er hele forskjellen på baklengs og forlengs Euler. Skriver du t[n], har du laget en tungvint variant av forlengs Euler.
- Svaret er sol.x. root returnerer et resultatobjekt, ikke et tall. sol.x er alltid et array — også når likningen har én ukjent — og det er nettopp den formen basisklassen vil ha.

Dokumentasjonen for root leveres som vedlegg på eksamen, så signaturen trenger du ikke pugge. Innpakningen må du kunne.

Bokas navn på oppgavetypen «utvid det utleverte løser-hierarkiet med en ny metode». Den er topp-differensierende og teller 3–5 poeng, ofte som siste deloppgave på det avsluttende settet.

Kjennetegnet i oppgaveteksten: du får en formel for ett steg og blir bedt om å implementere den «ved arv» eller «som en subklasse». Svaret er alltid en klasselinje, en advance og ingenting mer.

I løser-hierarkiet er arbeidsdelingen fast: basisklassen eier tidsløkka, subklassen eier formelen for ett steg.

En metode der den nye verdien regnes ut av verdier som allerede er kjent. Forlengs Euler, Heun, midtpunkt, Ralston og Runge–Kutta 4 er alle eksplisitte.

Fordelen er at hvert steg er én utregning. Ulempen er en stabilitetsgrense: for stive problemer må steget være svært lite, ellers sprenger løsningen.

En metode der den nye verdien står på begge sider av likhetstegnet, slik som baklengs Euler: $u_{n+1} = u_n + \Delta t f(t_{n+1}, u_{n+1})$.

Hvert steg krever at du løser en likning — i denne boka med root. Prisen er regnearbeid; gevinsten er stabilitet uansett steglengde.

Et problem der noe skjer svært raskt sammenlignet med tidsrommet du er interessert i — for eksempel $u' = -15(u - 20)$ studert over to tidsenheter.

En eksplisitt metode må da ta absurd små steg for å holde seg i ro. En implisitt metode slipper unna med store steg. Kjennetegnet i en oppgave er en stor negativ rate.

Mellomregninger inne i advance: høyresiden evaluert i punkter du foreløpig bare gjetter deg fram til. Antallet prøvesteg er omtrent metodens orden, og hvert prøvesteg koster ett kall til høyresiden.

Navnene k1, k2, k3 er standard. Bruk dem ikke om modellparametre samtidig — kall parametrene noe annet, som ka og kb.

Halvér steget og se hva feilen ganges med. Faktor 2 betyr orden 1, faktor 4 orden 2, faktor 8 orden 3 og faktor 16 orden 4.

Tallene treffer sjelden eksakt for grove steg; du ser på tendensen over flere halveringer, ikke på ett enkelt forholdstall.

En metode som bruker mer enn den siste verdien — for eksempel Adams–Bashforth med to steg, som også trenger u[n-1] og høyresiden der.

Fordi self.u inneholder hele historien, er det lett å slå opp bakover. Det aller første steget må gjøres med en énstegsmetode, siden u[-1] ville pekt på siste rad i arrayet.

Det arrayet root starter fra. I en implisitt ODE-metode er u[n] det naturlige valget: den forrige verdien ligger nær den nye når steget er lite.

En dårlig startgjetning kan gi sol.success lik False, og da er svaret i sol.x ikke til å stole på.

En funksjon definert inne i en annen funksjon eller metode. Den ser den ytre funksjonens lokale variabler, og det er nettopp derfor likning(v) defineres inne i advance: den trenger u[n], t[n+1] og dt, som skifter for hvert steg.

Definerer du likning utenfor, må du sende alle disse verdiene inn som argumenter — det virker, men blir tungvint.

Resultatet når advance glemmer return: metoden gir None, basisklassen lagrer det i et flyttallsarray, og alle radene fra og med den andre blir nan.

Symptomet er farlig fordi ingen exception kastes. Ser du nan i et løsningsarray, er den første mistanken en manglende return.

Sensorprinsippet som gjelder hele Del 7: riktig klasselinje, riktig advance-signatur og riktig oppbygging av formelen gir uttelling selv når en koeffisient eller en indeks er feil.

Å skrive plettfri kode uten å kunne kjøre den er nesten umulig, og oppgavene er poengsatt med det i tankene. Skriv derfor strukturen først, og fyll inn detaljene etterpå.