6.1 Numerisk derivasjon, integrasjon og nullpunkt

$$f'(x) \approx \frac{f(x+h) - f(x)}{h}$$

med en liten, fast $h$. Nevneren er h, og de to funksjonsverdiene er i x og x+h. Feilen er omtrent proporsjonal med $h$: halverer du $h$, halveres feilen. Metoden kalles førsteordens av samme grunn.

from math import sin, cos, pi

h = 1e-4
x = pi/3
print((sin(x+h) - sin(x))/h)
print(cos(x))

Utskrift:

0.4999566978958203
0.5000000000000001

$$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$$

med funksjonsverdier på hver side av $x$. Nevneren er $2h$, ikke $h$ — det er den vanligste feilen i formelen. Feilen er omtrent proporsjonal med $h^2$: tideler du $h$, blir feilen hundre ganger mindre. Metoden er andreordens, og den er derfor standardvalget i faget når oppgaven ikke sier noe annet.

from math import sin, cos, pi

h = 1e-4
x = pi/3
print((sin(x+h) - sin(x-h))/(2*h))
print(cos(x))

Utskrift:

0.49999999916616744
0.5000000000000001

Bokas standardform for den numeriske deriverte som funksjon: den tar funksjonen f som argument, punktet x, og en steglengde med standardverdi. Kroppen er én linje med den sentrerte formelen, og funksjonen returnerer ett tall.

from math import sin, pi

def diff(f, x, h=1e-5):
    return (f(x+h) - f(x-h))/(2*h)

print(diff(sin, pi/3))
print(diff(lambda x: x**3, 2.0))

Utskrift:

0.4999999999921733
12.00000000021184

Samme formel pakket i en klasse, slik at instansen selv blir en funksjon: __init__ husker funksjonen og steglengden, og __call__ regner ut den deriverte i et punkt. Klasseapparatet står i kap. 5.1; her er poenget bare at eksamen har brukt begge formene, og at de gir samme tall.

En instans er ett konkret objekt laget av klassen, og en instans er kallbar når den kan brukes med parenteser som en funksjon — det er __call__ som gjør den kallbar.

from math import sin, pi

class Diff:
    def __init__(self, f, h=1e-5):
        self.f = f
        self.h = h

    def __call__(self, x):
        f, h = self.f, self.h
        return (f(x+h) - f(x-h))/(2*h)

dsin = Diff(sin)
print(dsin(pi/3))

Utskrift:

0.4999999999921733

$$\int_a^b f(x)\,dx \approx \Delta x\left(\frac{f(a)}{2} + \frac{f(b)}{2} + \sum_{k=1}^{N-1} f(a + k\,\Delta x)\right), \qquad \Delta x = \frac{b-a}{N}$$

Endepunktene teller halvt, alle de indre punktene teller helt. Legg merke til grensene i summen: range(1, N) — den går fra 1 til $N-1$, altså bare de indre punktene. Feilen er proporsjonal med $\Delta x^2$, så en dobling av $N$ firedobler nøyaktigheten.

def trapes(f, a, b, N):
    dx = (b - a)/N
    s = 0.5*f(a) + 0.5*f(b)
    for k in range(1, N):
        s += f(a + k*dx)
    return s*dx

print(trapes(lambda x: x**2, 0.0, 1.0, 10))

Utskrift:

0.3350000000000001

$$x_{\text{ny}} = x - \frac{f(x)}{f'(x)}$$

Løkka kjører til funksjonsverdien er liten nok: while abs(f(x)) > eps:. Metoden er svært rask når startgjetningen er brukbar — antall riktige siffer omtrent dobles per runde. Den krever at du kan gi den deriverte som en egen funksjon, og den kan bomme helt hvis $f'(x)$ er nær null.

def newton(f, df, x, eps=1e-12):
    while abs(f(x)) > eps:
        x = x - f(x)/df(x)
    return x

print(newton(lambda x: x**2 - 2, lambda x: 2*x, 1.0))

Utskrift:

1.4142135623730951

$$x_{\text{ny}} = x_1 - f(x_1)\,\frac{x_1 - x_0}{f(x_1) - f(x_0)}$$

Metoden trenger to startverdier og ingen derivert, og den er nesten like rask som Newton. Legg merke til at oppdateringen bare bruker de to siste verdiene — dette er samme mønster som en differenslikning (kap. 6.2), og det løses med en løkke som skyver verdiene videre.

def sekant(f, x0, x1, eps=1e-12, maks=100):
    n = 0
    while abs(f(x1)) > eps and n < maks:
        x_ny = x1 - f(x1)*(x1 - x0)/(f(x1) - f(x0))
        x0, x1 = x1, x_ny
        n += 1
    return x1

print(sekant(lambda x: x**2 - 2, 1.0, 2.0))

Utskrift:

1.4142135623730954

Biblioteksfunksjonen som løser $f(x) = 0$, også når $x$ har flere komponenter. Tre ting må sitte:

- fun tar ett array-argument og returnerer ett array — én komponent per likning.
- Startgjetningen x0 er også et array.
- Løsningen ligger i sol.x, som er et array. sol.success sier om metoden lyktes.

Dokumentasjonen for root leveres alltid som vedlegg på eksamen, så du trenger ikke pugge signaturen — men du må kunne pakke funksjonen din riktig, altså array inn og array ut. Det er nettopp det oppgavene tester.

import numpy as np
from scipy.optimize import root

def fun(x):
    return np.array([x[0]**2 - 2])

sol = root(fun, np.array([1.0]))
print(sol.x, sol.success)

Utskrift:

[1.41421356] True

Eksamensoppgaven «skriv en kort funksjon som regner ut noe tilnærmet, og bruk den»: numerisk derivert, integral eller nullpunkt. Vekten er 2 til 8 poeng. Derivasjon er i 4 av 8 sett (50 %), nullpunkt i 3 av 8 (38 %, voksende), og integrasjon i 1 av 8 (13 %).

En førsteordens metode har feil omtrent proporsjonal med steglengden: halverer du steget, halveres feilen. En andreordens metode har feil proporsjonal med steget i annen: halverer du steget, blir feilen en firedel. Forlengs differanse er førsteordens, sentrert differanse og trapesmetoden er andreordens.

Standardmåten å vise at en metode virker: kjør den for flere steglengder, regn feilen mot det eksakte svaret, og se på faktoren mellom radene. Faktor 10 per tideling av $h$ betyr førsteordens, faktor 100 betyr andreordens. Har du ikke noe eksakt svar, halverer du steget og viser at svaret slutter å endre seg i de sifrene du bryr deg om.

Feilen i en differanseformel synker med h bare til et punkt. Når f(x+h) og f(x-h) blir nesten like, mister differansen gyldige siffer, og feilen begynner å vokse igjen. For sentrert differanse ligger vendepunktet rundt $h = 10^{-6}$, og derfor er h = 1e-5 bokas standardvalg.

Standardverdien gjør at diff(f, x) kan kalles uten å oppgi steglengden, men at diff(f, x, 1e-3) fortsatt virker. Default-argumenter står etter de vanlige argumentene i signaturen. Merk feil #21 i feilkatalogen: en default-verdi må ikke være en liste eller annen muterbar verdi.

Endepunktene teller halvt, alle indre punkter helt, og hele summen ganges med $\Delta x$. I kode: s = 0.5*f(a) + 0.5*f(b), deretter for k in range(1, N): s += f(a + k*dx), til slutt return s*dx. Grensen range(1, N) treffer nøyaktig de indre punktene.

Newton uten derivert: linjen gjennom de to siste punktene brukes i stedet for tangenten. Trenger to startverdier, og oppdateringen er x1 - f(x1)*(x1 - x0)/(f(x1) - f(x0)). De to verdiene skyves videre med x0, x1 = x1, x_ny. Nesten like rask som Newton, og uunnværlig når den deriverte ikke finnes på formel.

Skal du finne når $A(t) = B(t)$, setter du opp differansen $g(t) = A(t) - B(t)$ og finner nullpunktet til $g$. Alle tre metodene i kapitlet gjør det samme med $g$; det eneste du gjør, er å sette opp differansen riktig.

Sett løsningen tilbake i den opprinnelige likningen og se om venstre side blir null (eller at de to størrelsene blir like). Dette er den eneste måten å oppdage at en løser har konvergert mot et annet nullpunkt enn du var ute etter, og det gir poeng i seg selv i løsningsforslag.

Alle metodene i kapitlet tar funksjonen som skal behandles, som argument: diff(f, x), trapes(f, a, b, N), newton(f, df, x). Da må du sende inn selve funksjonen (f, sin eller en lambda), ikke et funksjonskall og ikke et uttrykk. Å sende f(x) i stedet for f er å sende ett tall — en klassisk flervalgsfelle.

En grov metode krever slakk toleranse. Sentrert differanse med h = 1e-5 treffer på rundt $10^{-11}$, så tol = 1e-6 er trygt. Forlengs differanse med samme h treffer bare på rundt $10^{-5}$, og da vil tol = 1e-6 få testen til å feile selv om koden er riktig.