Numerisk derivasjon#
Læringsutbytte
Etter å ha arbeidet med denne delen av emnet, skal du kunne:
forklare den teoretisk bakgrunnen for numerisk derivasjon, og forskjellen på numerisk og analytisk derivasjon
implementere framoverdifferansen, bakoverdifferansen og sentraldifferansen for numerisk derivasjon
derivere funksjoner og data
Derivasjonsbegrepet#
Derivasjon handler om endring. Den deriverte kan beskrive stigningen og forandringen i et forløp. Nærmere bestemt handler det om momentan endring, altså endringen mellom to tilstander (funksjonsverdier) over en svært liten endring i en annen tilstand:
Vi har derfor nytte av derivasjon i mange tilfeller der vi ønsker å beskrive en utvikling. Vi kan derimot ikke alltid derivere analytisk for hånd, så det er en stor verdi i å beherske numerisk derivasjon. Da kan vi enklest tilnærme den analytisk deriverte med en numerisk derivert:
Numerisk derivasjon (framoverdifferansen)
For en liten verdi av \(dx\) kan vi tilnærme den førstederiverte slik:
der vi tilnærmer grenseverdien med en svært liten verdi av dx. Her skal vi se på denne metoden og andre metoder som kan brukes til å tilnærme den deriverte numerisk.
Underveisoppgave
Bruk definisjonen ovenfor og regn ut \(f'(1)\) for \(f(x) = 2x + 2\). Sett \(dx = 1\cdot 10^{-8}\).
Løsningsforslag
def f(x):
return x**2
dx = 1E-8
x = 1
fder = (f(x + dx) - f(x))/dx
print("f'(1) =", fder)
Løsningsforslaget ovenfor viser en enkel måte å implementere den numeriske deriverte på. Det kan også være nyttig å kunne lage metoden som en funksjon:
def deriver(f, x, dx = 1E-8):
dy = f(x + dx) - f(x)
return dy/dx
Legg merke til at vi ikke deriverer symbolsk. Det betyr at vi ikke får en annen funksjon når vi deriverer en funksjon. Vi får bare funksjonsverdier. Vi må altså deriverere i gitte punkter, for eksempel \(f'(1)\) eller \(f'(-5)\). Dersom vi ønsker å visualisere den deriverte til en funksjon, må vi derfor derivere funksjonen i flere punkter. Dette kan vi gjøre vektorisert ved hjelp av arrayer:
import numpy as np
import matplotlib.pyplot as plt
def f(x):# Definerer en funksjon vi skal derivere
return 3*x**3 + x**2 - 1
x = np.linspace(-2,2,100)
y = f(x)
yder = deriver(f,x)
plt.plot(x,y,label="f(x)") # Plotter funksjonen
plt.plot(x,yder,label="f'(x)") # Plotter den deriverte funksjonen
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
Feilanalyse#
La oss nå ta en titt på hvilke verdier av \(\Delta x\) som gir best resultat. Det må vel være den verdien som ligger nærmest 0, altså en så liten verdi som mulig – eller? La oss teste dette ved å skrive ut den deriverte for ulike verdier av \(\Delta x\):
def f(x):
return 2*x**2 + x - 5
def fder_analytisk(x):
return 4*x + 1
x = 1
delta_x = [10**-i for i in range(1,18)] # liste med verdier fra 10^-18 til 10^-1
analytisk = fder_analytisk(x)
for i in range(len(delta_x)):
numerisk = deriver(f, x, delta_x[i])
feil = abs(numerisk-analytisk)/analytisk * 100
print("For delta_x =", delta_x[i],"er feilen:", feil, "%")
For delta_x = 0.1 er feilen: 4.000000000000092 %
For delta_x = 0.01 er feilen: 0.40000000000048885 %
For delta_x = 0.001 er feilen: 0.03999999998569592 %
For delta_x = 0.0001 er feilen: 0.003999999920267783 %
For delta_x = 1e-05 er feilen: 0.00040000068821655077 %
For delta_x = 1e-06 er feilen: 3.999977025159751e-05 %
For delta_x = 1e-07 er feilen: 4.010780685348436e-06 %
For delta_x = 1e-08 er feilen: 6.07747097092215e-07 %
For delta_x = 1e-09 er feilen: 8.274037099909037e-06 %
For delta_x = 1e-10 er feilen: 8.274037099909037e-06 %
For delta_x = 1e-11 er feilen: 8.274037099909037e-06 %
For delta_x = 1e-12 er feilen: 0.008890058234101161 %
For delta_x = 1e-13 er feilen: 0.07992778373591136 %
For delta_x = 1e-14 er feilen: 0.524016993585974 %
For delta_x = 1e-15 er feilen: 6.581410364015028 %
For delta_x = 1e-16 er feilen: 100.0 %
For delta_x = 1e-17 er feilen: 100.0 %
Vi ser at “store” verdier som 0.1 og 0.01 gir en del feil. Men vi ser også faktisk at nøyaktigheten er størst ved \(dx = 10^{-8}\), og at den synker både med økende og med minkende \(dx\). Og attpåtil gir \(dx \leq 10^{-16}\) null som svar! Dette gir naturlig nok en feil på 100 %, siden den analytiske verdien er 5.
Vi forventer kanskje ikke dette resultatet. Dersom vi kun ser på definisjonen av den deriverte, er det ikke spesielt logisk at det skal slå slik ut. Men det hele handler om at tall ikke er representert eksakt i en datamaskin, og når datamaskinen skal operere med svært små tall, kan det bli en liten avrundingsfeil når den regner med tallene. Denne avrundingsfeilen gjør at vi får feil dersom vi velger for små verdier av \(dx\). Dersom vi gjør en mer generell feilanalyse, viser det seg at \(10^{-8}\) er en god verdi å velge her.
Underveisoppgave
Lag et plott med feilen som funksjon av dx med utgangspunkt i programmet ovenfor. Bruk logaritmiske akser – dette får resultatene tydeligere fram. Du kan lage logaritmiske akser slik:
plt.yscale('log')
plt.xscale('log')
Løsningsforslag
import matplotlib.pyplot as plt
def f(x):
return 2*x**2 + x - 5
def fder_analytisk(x):
return 4*x + 1
x = 1
delta_x = [10**-i for i in range(1,18)] # liste med verdier fra 10^-18 til 10^-1
analytisk = fder_analytisk(x)
avvik = []
for i in range(len(delta_x)):
numerisk = deriver(f, x, delta_x[i])
feil = abs(numerisk-analytisk)/analytisk * 100
avvik.append(feil)
print("For delta_x =", delta_x[i], "er feilen:", feil, "%")
plt.plot(delta_x, avvik)
plt.yscale('log')
plt.xscale('log')
plt.xlabel('dx')
plt.ylabel('Feil (%)')
plt.show()
Andre tilnærminger#
Tilnærmingen til den deriverte som vi har sett på, tar utgangspunkt i punktene \((x, f(x))\) og \((x+dx, f(x+dx))\) for å regne ut den momentane vekstfarten (altså den deriverte) i punktet \(x\). Men vi kan like godt bruke andre punkter. Metoden vi har sett på, kalles framoverdifferansen fordi den tar utgangspunkt i punktet \(x\) og neste punkt \(x + dx\). Tilsvarende kan vi ta utgangspunkt i punktet \(x\) og forrige punkt \(x - dx\). Dette kaller vi bakoverdifferansen. Bakoverdifferansen gir samme feil som framoverdifferansen, men er teoretisk nyttig for å utlede andre metoder.
Numerisk derivasjon (bakoverdifferansen)
For en liten verdi av \(dx\) kan vi tilnærme den førstederiverte slik:
En metode som derimot gir mindre feil enn både framover- og bakoverdifferansen, er sentraldifferansen. Det er en slags kombinasjon av framover- og bakoverdifferansen, der vi tar utgangspunkt i gjennomsnittet (midtpunktet) av \(x+dx\) og \(x-dx\).
Numerisk derivasjon (sentraldifferansen)
For en liten verdi av \(dx\) kan vi tilnærme den førstederiverte slik:
Underveisoppgave
Implementer bakover- og sentraldifferansen som Python-funksjoner. Gjør en feilanalyse med de tre ulike tilnærmingene for ulike verdier av \(\Delta x\). Bruk funksjonen \(f(x) = \sin{(x)}\) og sammenlikn med den analytiske verdien av den deriverte, \(f'(x) = \cos{(x)}\)
Numerisk derivasjon av data#
Nå kommer vi til den nyttigste delen av numerisk derivasjon, nemlig derivasjon av diskrete data. Vi kan derivere på samme måte som vi gjorde med kontinuerlige funksjoner, men vi har gitt en \(dx\) som er gitt av avstanden mellom datapunktene våre. Hvis målefrekvensen er lav, blir \(dx\) høy, og motsatt. La oss se på hvordan vi kan derivere noen posisjonsdata:
import pandas as pd
# Leser og sjekker ut dataene
data = pd.read_csv("https://raw.githubusercontent.com/andreasdh/NAT3000/master/docs/datafiler/posisjon.txt")
data.head()
| tid_s | posisjon_m | |
|---|---|---|
| 0 | 0.00000 | 0.000000 |
| 1 | 0.01001 | 0.000060 |
| 2 | 0.02002 | 0.000240 |
| 3 | 0.03003 | 0.000541 |
| 4 | 0.04004 | 0.000962 |
# Plotter dataene
import matplotlib.pyplot as plt
x = data["tid_s"]
y = data["posisjon_m"]
plt.scatter(x, y)
plt.title("Bevegelse av legeme")
plt.xlabel("Tid (s)")
plt.ylabel("Posisjon (m)")
plt.show()
Husk at \(v(t) = s'(t)\). Hvis vi deriverer, får vi altså farten til legemet:
dydx = []
for i in range(len(y)-1):
dy = y[i+1]-y[i]
dx = x[i+1]-x[i]
der = dy/dx
dydx.append(der)
dydx.append(None)
plt.plot(x, y, label='Posisjon (m)')
plt.plot(x, dydx, label='Fart (m/s)')
plt.xlabel('Tid (s)')
plt.legend()
plt.show()
I løkka bruker vi samme framgangsmåte som vi gjorde da vi deriverte funksjoner, men nå tar vi verdiene fra ei liste med verdier. Vi kan kun kjøre løkka til lengden av lengden av listene minus 1, fordi vi skal ta forskjellen mellom verdier. Da blir det nemlig én verdi mindre enn i y- og x-listene. Derfor legger vi til en ekstra verdi “None” til slutt i lista, slik at alle listene blir like lange.
Som du ser, kan vi benytte numerisk derivasjon på både kontinuerlige funksjoner og diskrete data. Hovedpoenget er at vi finner ut noe om endringen i en funksjon eller i et datasett. Og desto mindre dx er, desto bedre tilnærming er denne endringen til den momentante endringen i et punkt, altså den deriverte.
Didaktisk oppgave
Kan programmering og numeriske metoder være en måte å forstå derivasjon på? I så fall hvordan?
I denne delen har du sett på enkel numerisk derivasjon, feilanalyse, ulike tilnærminger til den deriverte og derivasjon av diskrete data. Hvilke deler av dette egner seg for elever, og i hvilke(n) sammenheng(er)?
Derivasjon med numeriske biblioteker#
Selvfølgelig har noen også lagd funksjoner som kan derivere for oss. Vi finner dem i scipy-biblioteket. Funksjonen derivative finner den numerisk deriverte ved hjelp av midtpunktstilnærmingen:
from scipy.misc import derivative
import numpy as np
def f(x):
return x**2 - 4
derivert = derivative(f, 1)
print("f'(1) =", derivert)
f'(1) = 2.0
Oppgaver#
Oppgave 1
Beregn f’(1) numerisk for følgende funksjoner. Kontroller ved å derivere for hånd.
\(f(x) = x^2 - 4x + 5\)
\(f(x) = e^x\)
\(f(x) = \sqrt{\ln(x)}\)
Oppgave 2
Skriv om funksjonen deriver slik at dx har en standardverdi. Velg en standardverdi som sannsynligvis vil gi gode resultater.
def deriver(f, x, dx = 1E-8):
dy = f(x + dx) - f(x)
return dy/dx
Oppgave 3
En partikkel følger posisjonsfunksjonen \(x(t) = t^3 + \frac{1}{3}\cdot t\). Plott både posisjon, hastighet og akselerasjon for \(t = [0,10]\). Husk at \(v(t) = s'(t)\) og \(a(t) = v'(t)\).
Oppgave 4
Forklar hva som er forskjellen mellom analytisk og numerisk derivasjon.
Oppgave 5
Temperaturen \(T(t)\) (i celsius) etter \(t\) minutter til en nylig lagd te følger denne modellen, som vi har kommet fram til ved regresjon av eksperimentelle data:
a) Bruk Newtons kvotient for å tilnærme den deriverte til \(T(t)\) for 1000 jevnt fordelte verdier av \(t\) i intervallet \([0, 60]\) b) Plott resultatet fra a) sammen med grafen for \(T(t)\). Bruk merkelapper (labels og legend) for grafene. c) Hvor mye synker temperaturen med ved element nr. 42 i lista over verdier for \(t\) (vi teller med 0-te element)?
Oppgave 6
Programmet nedenfor leser av fila “heistur_kjemi_fysikk.txt” og finner fart og akselerasjon ved hjelp av numerisk derivasjon. Programmet fungerer derimot ikke helt som det skal. Rett opp feilen. Lag også en ny kolonne “fart” og en ny kolonne “akselerasjon” i dataframen.
Løsningsforslag
import pandas as pd
import matplotlib.pyplot as plt
data = pd.read_csv("heistur_kjemi_fysikk.txt")
p = data["height_m"]
t = data["time_s"]
v = [] # fart i m/s
a = [] # akselerasjon i m/s^2
# Vi trenger to løkker fordi a tar utgangspunkt i v, som derfor må lages først
for i in range(len(p)-1):
dt = t[i+1] - t[i]
fart = (p[i+1] - p[i])/dt
v.append(fart)
for i in range(len(v)-1):
dt = t[i+1] - t[i]
akselerasjon = (v[i+1] - v[i])/dt
a.append(akselerasjon)
v.append(None)
a.append(None)
a.append(None)
# Legger til verdiene i dataframen
data["fart"] = v
data["akselerasjon"] = a
plt.subplot(3,1,1)
plt.ylabel("Posisjon (m)")
plt.plot(t, p, color = "limegreen")
plt.subplot(3,1,2)
plt.ylabel("Fart (m)")
plt.ylim(-1,1)
plt.plot(t, v, color = "navy")
plt.subplot(3,1,3)
plt.ylabel("Akselerasjon (m)")
plt.plot(t, a, color = "firebrick")
plt.ylim(-1,1)
plt.xlabel("Tid (s)")
plt.tight_layout()
plt.show()
Oppgave 7
Løs puslespillet nedenfor, som skal illustrere derivasjon av diskrete data.