7.1 ODE-løser som funksjon: forward Euler

$$u'(t) = f(t, u(t)), \qquad u(t_0) = u_0$$

Funksjonen $f$ kalles høyresiden: den tar tidspunktet $t$ og den nåværende verdien $u$, og returnerer stigningstallet der. I Python er høyresiden en helt vanlig funksjon med to argumenter, og rekkefølgen er t først:

def f(t, u):
    return -0.1*(u - 20)

Legg merke til at t ikke behøver å brukes i kroppen. Argumentet skal likevel være med, fordi løseren alltid kaller f(t, u) med to argumenter.

Startverdien er data, ikke en del av f. Den sendes til løseren separat.

Lager et array med jevnt fordelte tall fra a til b, der begge endepunktene er med. Det tredje argumentet er antall PUNKTER, ikke antall steg.

Skal intervallet deles i N like store steg, trenger du derfor N+1 punkter:

import numpy as np
print(np.linspace(0, 1, 5))

Utskrift:

[0.   0.25 0.5  0.75 1.  ]

Fire steg på 0,25 — fem punkter. Dette er kilden til den hyppigste feilen i faget, og den har eget nummer i feilkatalogen: #1 — av-én-feilen på antall punkter.

Lager et nytt array med samme form og samme datatype som t, fylt med nuller. Brukes til å allokere løsningsarrayet, så det garantert blir like langt som tidsarrayet:

import numpy as np
t = np.linspace(0, 1, 4)
u = np.zeros_like(t)
print(u, u.shape)

Utskrift:

[0. 0. 0. 0.] (4,)

Alternativet np.zeros(N+1) gjør det samme når du vet tallet. Fordelen med zeros_like er at formen ikke kan komme i utakt med t.

Steglengden lest ut av tidsarrayet, i stedet for regnet ut på nytt. Den er alltid riktig når t kommer fra np.linspace, og du unngår å blande antall punkter og antall steg.

import numpy as np
t = np.linspace(0, 3, 31)
print(t[1] - t[0], 3/30, 3/31)

Utskrift:

0.1 0.1 0.0967741935483871

De to første er like, den tredje er den klassiske av-én-feilen (#1). Merk symbolkollisjonen: i numerisk derivasjon het steget h (kap. 6.1), her heter det dt fordi det er et steg i tid. Det er samme slags størrelse.

En metode har konvergensorden $p$ hvis feilen oppfører seg som $\Delta t^p$ når steget blir lite. Praktisk konsekvens, og det du måler i en tabell:

- orden 1: halvér $\Delta t$ $\to$ feilen halveres (faktor 2)
- orden 2: halvér $\Delta t$ $\to$ feilen firedeles (faktor 4)
- orden 4: halvér $\Delta t$ $\to$ feilen deles på 16

Forlengs Euler har orden 1. Det er derfor den er billig per steg, men trenger mange steg for god nøyaktighet.

$$u_{n+1} = u_n + \Delta t \cdot f(t_{n+1}, u_{n+1})$$

Den ukjente står på begge sider — metoden er implisitt. For den spesielle høyresiden $f = -\lambda u$ kan du løse for $u_{n+1}$ med algebra:

$$u_{n+1} = u_n - \Delta t \lambda u_{n+1} \quad \Longrightarrow \quad u_{n+1}(1 + \Delta t \lambda) = u_n \quad \Longrightarrow \quad u_{n+1} = \frac{u_n}{1 + \Delta t \lambda}$$

I kode blir det én linje:

u[n+1] = u[n]/(1 + dt*lam)

Metoden har også orden 1, men den er stabil for alle steglengder: nevneren er alltid større enn 1, så løsningen kan aldri vokse. For en generell høyreside må du løse en likning i hvert steg — det gjør vi i kap. 7.4.

$$k_1 = f(t_n, u_n), \qquad k_2 = f(t_n + \Delta t,\; u_n + \Delta t k_1), \qquad u_{n+1} = u_n + \frac{\Delta t}{2}(k_1 + k_2)$$

I kode er det tre linjer inne i løkka:

k1 = f(t[n], u[n])
k2 = f(t[n] + dt, u[n] + dt*k1)
u[n+1] = u[n] + dt*(k1 + k2)/2

Utskrift: ingen (fragmentet hører inne i løkka i eksempelet under).

Metoden har orden 2: halvert steg gir en firedel av feilen. Prisen er to kall til f per steg i stedet for ett.

Python-funksjonen som gir stigningstallet i et initialverdiproblem. Tar to argumenter — tid først, tilstand etterpå — og returnerer $u'$ i det punktet.

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

Boka og alle sett fra 2023 og senere bruker denne rekkefølgen. Argumentet t skal med selv om kroppen ikke bruker det.

Avstanden i tid mellom to nabopunkter i løsningen. Regnes ut som t[1] - t[0] eller T/N, aldri T/(N+1).

Mindre dt gir mindre feil, men flere steg. For forlengs Euler er feilen omtrent proporsjonal med dt.

Bokas viktigste tellergel. Deler du et intervall i N like store steg, er det N+1 punkter fordi startpunktet også teller.

Konsekvens: np.linspace(0, T, N+1), np.zeros(N+1) og for n in range(N). Feil her er feilkode #1, den hyppigste i faget.

En metode der neste verdi regnes direkte ut av verdier du allerede har. Forlengs Euler og Heun er eksplisitte.

Fordel: én enkel oppdateringslinje. Ulempe: steget må være lite nok, ellers sprenger løsningen (se stabilitet).

En metode der den ukjente u[n+1] står på begge sider av oppdateringen, slik som baklengs Euler.

Da må hvert steg løse en likning. For $f = -\lambda u$ kan du løse den med algebra; for en generell høyreside trenger du en likningsløser (kap. 7.4).

$$\Delta t < \frac{2}{\lambda}.$$

Brytes det, veksler verdiene fortegn og vokser. Baklengs Euler er stabil for alle steglengder.

Et problem der en av prosessene går mye raskere enn de andre, typisk stor $\lambda$. Da tvinger stabilitetskravet en eksplisitt metode til å bruke veldig små steg selv om løsningen er kjedelig.

Implisitte metoder er svaret på stive problemer.

Startverdien settes utenfor løkka, og løkka fyller bare inn u[1] … u[N].

Settes den etter løkka, regner løkka med nuller hele veien og du får et array av nuller med startverdien limt på. Det er en av de klassiske strukturfeilene i denne sjangeren.

En løser returnerer både tidspunktene og løsningen, som en tuppel. Kalleren pakker ut med to navn:

t, u = forward_euler(f, 3, 30, 1)
print(len(t), len(u))

Utskrift:

31 31

Uten t kan du verken plotte eller finne ut når noe skjer.

Standardmålet på hvor god en numerisk løsning er: det største avviket over alle punktene.

Både u og eksakt er arrays, så subtraksjonen skjer elementvis og np.abs virker på hele arrayet. Bruk dette, ikke feilen i ett tilfeldig punkt.

Siste element i løsningsarrayet, altså verdien i sluttiden T. Negativ indeks teller fra slutten.

Mange oppgaver spør nettopp «hva er verdien ved $t = T$» — svaret er u[-1], ikke u[N-1].

At f sendes inn til løseren, er det som gjør løseren gjenbrukbar: samme funksjon løser alle initialverdiproblemer.

Send funksjonsnavnet uten parenteser (forward_euler(f, ...)). Skriver du f(t, u) i kallet, sender du et tall, og løseren krasjer med TypeError når den prøver å kalle det.

Slik måler du orden uten teori: kjør med N og 2N steg og se på forholdet mellom feilene.

Forhold ca. 2 betyr orden 1, ca. 4 betyr orden 2, ca. 16 betyr orden 4. Metoden fungerer bare hvis du har en eksakt løsning å måle mot.

Mellomresultatene i en to-stegs metode: k1 er stigningstallet der du står, k2 stigningstallet i et prøvepunkt.

Heun bruker gjennomsnittet av de to. Midtpunktsmetoden bruker k2 alene, målt et halvt steg fram. Begge har orden 2.

En verdi $u^*$ med $f(t, u^*) = 0$ gir stigningstall null: løsningen står stille hvis den starter der.

Rask fornuftssjekk av en ODE-besvarelse: regn ut f(0, u0). Er den null, skal løsningen være konstant — og er den det, virker koden.

Blander du de to (range(N+1)), treffer siste tilordning u[N+1], som ikke finnes: IndexError.