Differensiallikninger og kontinuerlige modeller#
Læringsutbytte
Etter å ha arbeidet med dette temaet, skal du kunne:
Forklare hva som menes med differensiallikninger, og hvordan slike likninger kan brukes til å modellere ulike systemer.
Forklare og implementere Eulers metode (Forward Euler).
Løse ordinære differensiallikninger med egne Python-funksjoner og funksjoner fra Scipy-biblioteket.
Differensiallikninger#
Vi har tidligere sett på modeller som beskriver endringer ved bestemte tidssteg, for eksempel den enkleste modellen for utviklingen i antall smittede ved et influensautbrudd:
Vi kan beskrive modellen ovenfor med ord slik: Antall smittede ved neste tidspunkt er lik antall smittede ved forrige tidspunkt pluss en endring i antall smittede. Endringen er her gitt ved en vekstfaktor ganger antall smittede ved forrige tidspunkt. Det betyr at endringen er proporsjonal med antall smittede. Dette er et eksempel på en differenslikning, der vi finner en tallfølge, som for eksempel beskriver antall smittede i en populasjon eller antall gauper i et økosystem. Differenslikninger kan brukes der vi kan nøye oss med å beskrive systemet vårt på diskrete, altså bestemte, tider. Men hva hvis vi ønsker å beskrive systemet på ethvert tidspunkt, altså den momentane endringen? Da må vi bruke den deriverte. La oss første skrive differenslikningen på en annen form:
Nå har vi et uttrykk for endringen i systemet. Delta (\(\Delta\)) betyr her endring eller forskjell. Dersom denne endringen skjer mellom to tidspunkter som er så nær hverandre som mulig (\(\Delta t\) nærmer seg null), får vi den momentane endringen. Dette skriver vi slik:
Dette er ikke lenger en differenslikning, for her beskriver vi systemet med en kontinuerlig funksjon, og ikke funksjonsverdier på bestemte tidspunkter. Vi kaller slike likninger for differensiallikninger. En differensiallikning (også kalt difflikning) er en likning som inneholder den deriverte av en funksjon. Det vil si at differensiallikninger beskriver den momentane endringen i et system, altså endringen på ethvert tidspunkt. Når vi løser en differensiallikning, får vi en funksjon, eller funksjonsverdiene til en funksjon.
Differensiallikninger er svært anvendelige. De brukes til å beskrive endring i alt fra molekyler i en kjemisk reaksjon til bakterier i en petriskål eller bevegelsen til satelitter og planeter. De brukes også til å modellere overføring av informasjon i hjernen, smitteutvikling, økonomisk vekst, inntekt og investeringer, klima, miljø og vær. De fleste difflikninger er ikke løsbare for hånd. Derfor er det svært nyttig å kunne løse dem numerisk.
Differensiallikninger
En differensiallikning er en likning som inneholder den deriverte av en funksjon. I de fleste praktiske situasjoner beskriver slike likninger sammenhengen mellom endringen, \(f'(t)\), og tilstanden, \(f(t)\), til et system ved tida \(t\).
Det vil altså si at en difflikning er en sammenheng mellom den deriverte til en funksjon, altså endringen i funksjonen, og funksjonen selv. Når vi løser en difflikning, får vi ikke et uttrykk for \(x\), som når vi løser likninger – vi får et uttrykk for \(f(x)\). Siden vi som oftest løser difflikninger på datamaskinen, får vi heller ikke et uttrykk for \(f(x)\), men vi får funksjons_verdier_.
Siden både endringen og tilstanden i et system er avhengig av hvordan systemet er til å begynne med, trenger vi alltid initialbetingelser (startbetingelser) for å løse difflikninger. For eksempel trenger vi antall smittede ved t = 0 for å kunne regne ut antall smittede ved neste tidssteg. Vi må derfor alltid oppgi initialbetingelser når vi skal simulere et system.
Eksempler på differensiallikninger#
Dersom du møter på differensiallikninger i matematikk, er de ofte satt opp som enkle uttrykk, for eksempel slik:
\(y' = 1\)
\(y' = y\)
\(y' - 2x = -1\)
Dersom likningene betyr noe praktisk, pleier vi å skrive dem litt annerledes. Her er eksempler på differensiallikninger som kan brukes til å beskrive ulike systemer:
Bakterievekst: \(B'(t) = kB(t)\)
CO\(_2\)-utslipp: \(U'(t) = aU(t)\left(1 - \frac{U(t)}{b}\right)\)
Fallende ball: \(v'(t) = -g - \frac{kv(t)^2}{m}\)
Underveisoppgave
Studer differensiallikningene ovenfor og forklar med ord hva de beskriver (både overordnet og hvert ledd for seg).
Legg merke til at modellene ikke kan bevises eller utledes. Det er matematiske modeller som vi kan ha en hypotese om at stemmer godt med virkeligheten. Hvis vi for eksempel tror at veksten til bakteriene vil flate ut, justerer vi dette ved å legge til et ledd som gir logistisk vekst. Modellene må så etterprøves ved for eksempel å sammenlikne med reelle data. Modellering er derfor en prosess der vi kontinuerlig sammenlikner modell med virkelighet for å forbedre modellen.
Newtons 2. lov#
Et annet viktig eksempel på en differensiallikning er Newtons 2. lov. Newtons 2. lov er grunnlaget for simuleringer av alt fra satellitt- og planetbaner til molekyler, som vi snart skal se eksempler på. Loven sier at summen (\(\Sigma\)) av krefter som virker på et legeme er lik produktet av massen og akselerasjonen til legemet:
Dersom vi har en modell som beskriver kreftene som virker på et legeme, kan vi derfor regne ut akselerasjonen. Og siden \(a(t) = v'(t) = s''(t)\), kan vi formulere Newtons 2. lov som en differensiallikning. Da kan vi løse denne likningen og finne farten og posisjonen til et legeme på en hvilken som helst tid, gitt en startposisjon og en startfart (initialverdier).
Vi kan til og med formulere dette som en difflikning som inneholder den andrederiverte. Dette kaller vi en andreordens difflikning:
La oss modellere et enkelt system, for eksempel en fallende kule. I modellen vår ønsker vi å ta hensyn til både gravitasjon og luftmotstand. Da er \(\sum F = - G + L\), dersom vi velger positiv retning oppover. Vi skriver også akselerasjonen som en funksjon, siden den varierer med tida. Da får vi at:
\(\sum F = ma(t)\)
\(a(t) = \frac{\sum F}{m}\)
\(a(t) = - \frac{G}{m} + \frac{L}{m}\)
Deretter må vi sette inn modeller for gravitasjonskraften og luftmotstanden. Gravitasjonen på jorda kan enkelt beskrives med modellen \(G = mg\), der \(m\) er massen i kg og \(g \approx 9.8\) m/s\(^2\), som er tyngdeakselerasjonen. En enkel modell for luftmotstand er \(L = -kv\), der k er en konstant som varierer med legemets form og massetetthet, og \(v\) er farten til legemet. Konstanten \(k\) må bestemmes eksperimentelt.
\(a(t) = \frac{mg}{m} - \frac{kv(t)}{m}\)
\(v'(t) = g - \frac{kv(t)}{m}\)
Her ser vi at den siste likningen inneholder den deriverte av farten (altså akselerasjonen). Altså er dette en differensiallikning.
Eulers metode#
Vi har altså en startverdi \(f(x_0)\) og ønsker å finne \(f(x_0 + dx)\), altså neste funksjonsverdi. I tillegg har vi et uttrykk for den deriverte av funksjonen, nemlig differensiallikningen. Du kjenner faktisk allerede til et uttrykk som inneholder en funksjon og dens deriverte, nemlig definisjonen av den deriverte. La oss bruke den til å utlede en metode for å finne \(f(x + dx)\). Vi bruker den numeriske definisjonen der vi tilnærmer grenseverdiene med en dx (\(\Delta x\)) som er så liten som mulig:
Til å begynne med kjenner vi \(f(x)\), altså \(f(x_0)\). Dette er initialbetingelsen, for eksempel startkonsentrasjonen \(c(t_0)\) i eksempelet ovenfor. Vi kjenner også et uttrykk for den deriverte, nemlig differensiallikningen. I tillegg bestemmer vi selv tidssteget dx, men husk at det verken bør være for lite eller for stort. Den eneste ukjente i den numeriske tilnærmingen til den deriverte er altså \(f(x+dx)\). Det er jo nettopp \(f(x+dx)\) vi prøver å finne, fordi det beskriver tilstanden til systemet ved neste tidssteg. Med litt enkel algebra får vi omformet uttrykket slik at det blir et uttrykk for \(f(x+dx)\). Vi ganger først med \(dx\) på begge sider:
Deretter får vi \(f(x+dx)\) aleine på høyre side og ender opp med følgende likning:
Dette er Eulers metode, eller nærmere bestemt Forward Euler. Metoden kalles dette fordi den tar utgangspunkt i framoverdifferansen til den deriverte. Den brukes til å løse differensiallikninger, det vil si å integrere den deriverte slik at vi finner funksjonsverdiene. Siden vi ofte har å gjøre med funksjoner som varierer med tid, kaller vi gjerne dx for dt.
Eulers metode er en iterativ algoritme. Vi starter derfor med \(f(x_0)\) og finner de påfølgende funksjonsverdiene slik:
Eulers metode (Forward Euler)
Vi kan finne funksjonsverdiene \(f(t_{k+1})\) ved å bruke funksjonsverdien \(f(t_k)\) og den deriverte av funksjonen ved tida \(t_k\), \(f'(t_k)\) sammen med en steglengde \(dt\) som representerer en liten \(\Delta t\).
Implementering#
Vi starter med et enkelt eksempel der vi løser difflikningen \(f'(x) = x\) (som er det samme som \(y' = x\)). Vi gjør det først med en while-løkke der vi appender til lister for å spare på og til slutt plotte verdiene. Stort sett er det sammme framgangsmåte hver gang vi løser en difflikning:
Definer initialbetingelser (for eksempel \(y_0\)).
Definer hvilke verdier vi skal evaluere funksjonen i, altså et intervall \([x_0, x_{slutt}]\). Husk at vi bare kan derivere og integrere numerisk i bestemte punkter. Vi får altså ikke funksjoner når vi integrerer en differensiallikning – vi får funksjonsverdier.
Definer en steglengde. Denne bør være “ganske” liten, som da vi så på numerisk derivasjon. Som regel er \(dx = 1\cdot 10^{-8}\) en god verdi, men du kan gjerne kjøre med for eksempel \(dx = 1\cdot 10^{-3}\), slik at det ikke skal ta så lang tid å kjøre programmet.
Lag lister eller arrayer og legg initialbetingelser og startverdier inn i disse. Disse trenger vi for å kunne plotte resultatene til slutt. Dersom du bare er interessert i slutt-tilstanden, trenger du ikke dette trinnet.
Lag ei løkke som går igjennom hele det definerte intervallet:
Regn ut neste verdi av endringen \(f'(x)\), gitt av differensiallikninga vår.
Løs difflikninga med Eulers metode (eller en annen metode, når du lærer om det).
Oppdater den uavhengige variabelen (x, t eller liknende) med tidssteget.
Skriv ut eller visualiser resultatene.
import matplotlib.pyplot as plt
y0 = 1 # intitialbetingelse for y
x0 = 0 # startverdi for x
x_slutt = 5.0 # sluttverdi for x
dx = 1E-5 # steglengde
x_liste = [x0] # Legger inn første x-verdi i en liste for å spare på verdiene
y_liste = [y0] # Legger inn initalbetingelsen for y i en liste for å spare på verdiene
# Disse variablene skal oppdateres i løkka, mens x0 og y0 holdes konstante
x = x0
y = y0
while x < x_slutt:
yder = x # Difflikningen vi skal løse
y = y + yder*dx # Eulers metode for å finne neste y-verdi
x = x + dx # Oppdaterer x-verdien med steglengden
x_liste.append(x)
y_liste.append(y)
plt.plot(x_liste,y_liste)
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()
For å rydde opp litt i programmet vårt, kunne vi også definert difflikninga som en funksjon. Dette gjør det litt enklere å finne igjen og endre likningen.
import matplotlib.pyplot as plt
import numpy as np
y0 = 1 # intitialbetingelse for y
x0 = 0 # startverdi for x
x_slutt = 5.0 # sluttverdi for x
dx = 1E-5 # steglengde
x_liste = [x0] # Legger inn første x-verdi i en liste for å spare på verdiene
y_liste = [y0] # Legger inn initalbetingelsen for y i en liste for å spare på verdiene
# Disse variablene skal oppdateres i løkka, mens x0 og y0 holdes konstante
x = x0
y = y0
# Definerer difflikningen som en funksjon
def fder(x):
return x
while x < x_slutt:
yder = fder(x) # Difflikningen vi skal løse
y = y + yder*dx # Eulers metode for å finne neste y-verdi
x = x + dx # Oppdaterer x-verdien med steglengden
x_liste.append(x)
y_liste.append(y)
plt.plot(x_liste,y_liste)
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()
Underveisoppgave
Modifiser programmet ovenfor slik at du løser difflikningen \(f'(x) = \frac{1}{y+1} - lg(x+2)\). Bruk samme initialbetingelse og løsningsintervall.
Løsningsforslag
Her er det bare å endre funksjonen slik at den representerer den nye difflikninga. Husk å importere tierlogaritmen.
import matplotlib.pyplot as plt
import numpy as np
y0 = 1 # intitialbetingelse for y
x0 = 0 # startverdi for x
x_slutt = 5.0 # sluttverdi for x
dx = 1E-5 # steglengde
x_liste = [x0] # Legger inn første x-verdi i en liste for å spare på verdiene
y_liste = [y0] # Legger inn initalbetingelsen for y i en liste for å spare på verdiene
# Disse variablene skal oppdateres i løkka, mens x0 og y0 holdes konstante
x = x0
y = y0
# Definerer difflikningen som en funksjon
def fder(x):
return 1/(y+1) - np.log10(x+2)
while x < x_slutt:
yder = fder(x) # Difflikningen vi skal løse
y = y + yder*dx # Eulers metode for å finne neste y-verdi
x = x + dx # Oppdaterer x-verdien med steglengden
x_liste.append(x)
y_liste.append(y)
plt.plot(x_liste,y_liste)
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()
Vi kan også bruke arrayer istedenfor lister. Dette gir en raskere kode, og kan være en fordel i mange tilfeller. Da appender vi ikke til lister, men bruker indeksering. Vi bruker også gjerne en for-løkke, og da må vi vite hvor mange ganger løkka skal gå. Dette kan vi regne ut ved å ta differansen mellom \(x_0\) og \(x_{slutt}\), som vi dividerer på steglengden \(dx\). Da får vi antall intervaller mellom \(x_0\) og \(x_{slutt}\) med avstand \(dx\). For å finne antall punkter, må vi derfor ta antall intervaller og legge til 1:
\(N_{punkter} = \frac{x_{slutt} - x_0}{dx} + 1\)
Underveisoppgave
Forklar hvorfor antallet punkter N kan regnes ut som ovenfor.
import matplotlib.pyplot as plt
import numpy as np
# Plotteverdier
x0 = 0 # startverdi for x
x_slutt = 5.0 # sluttverdi for x
dx = 1E-5 # steglengde
N = int((x_slutt - x0)/dx) + 1 # antall punkter
# Initialbetingelser
y0 = 1 # intitialbetingelse for y
# Differensiallikningen
def yder(y,x):
return x
# Lager arrayer for å spare på verdiene
x = np.zeros(N)
y = np.zeros(N)
y[0] = y0 # Legger inn initalbetingelsen for y
x[0] = x0
# Eulers metode
for i in range(N-1):
y[i+1] = y[i] + yder(y[i],x[i])*dx
x[i+1] = x[i] + dx
plt.plot(x,y)
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.show()
Underveisoppgave
Forklar programmet ovenfor linje for linje. Modifiser også programmet slik at det løser difflikninga \(f'(x) = 1\).
Ulike initialbetingelser#
Dersom vi endrer initialbetingelsene, får vi naturlig nok andre løsninger. Her har vi et program som regner ut og plotter ut
# Plotteverdier
x0 = 0 # startverdi for x
x_slutt = 5.0 # sluttverdi for x
dx = 1E-3 # steglengde
N = int((x_slutt - x0)/dx) + 1 # antall punkter
# Initialbetingelser
y0 = 1 # intitialbetingelse for y
# Differensiallikningen y' = 1
def yder(y,x):
return 1
# Lager arrayer for å spare på verdiene
x = np.zeros(N)
y = np.zeros(N)
x[0] = x0
y_liste = [0, 1, 2, 3]
for y0 in y_liste:
y[0] = y0
for i in range(N-1):
y[i+1] = y[i] + yder(y[i],x[i])*dx
x[i+1] = x[i] + dx
plt.plot(x,y,label=f"y0={y0}")
plt.grid()
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Det er enkelt å se at funksjonene vi får ovenfor er \(y = x + C\), der \(C\) avhenger av initialbetingelsene. Vi får derimot ikke selve funksjonsuttrykket som output, men vi får, som nevnt før, funksjonsverdiene. Vi kan også løse likningen \(y' = y + x\) på samme måte. Den analytiske løsningen for denne likningen er \(y = Ce^x - x - 1\). Resultatene blir da slik for de samme initialbetingelsene som ovenfor:
def yder(y,x):
return y + x
y_liste = [0, 1, 2, 3]
for y0 in y_liste:
y[0] = y0
for i in range(N-1):
y[i+1] = y[i] + yder(y[i],x[i])*dx
x[i+1] = x[i] + dx
plt.plot(x,y,label=f"y0={y0}")
plt.grid()
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Underveisoppgave
Løs differensiallikninga \(y' + y^3 = \sin(x)\) med Eulers metode fra \(x = 0\) til \(x = 10\). Bruk initialbetingelsen \(y(0) = 0\). Vær nøye med hvordan du definerer funksjonen.
Kontinuerlige modeller#
Modeller som er basert på differensiallikninger, er kontinuerlige. Det betyr at de har funksjonsverdier for alle verdier av x (eller tilsvarende variabler). Nå skal det sies at når vi simulerer med utgangspunkt i kontinuerlige matematiske modeller på en datamaskin, blir resultatene ikke helt kontinuerlig. Det er fordi datamaskinen ikke kan operere med uendelig mange eller uendelig små (eller store) verdier. Vi sier at vi diskretiserer problemet når vi benytter en numerisk algoritme for å løse det. Men dersom stegene er tilstrekkelig små, er det en god tilnærming til en kontinuerlig løsning.
Vi har allerede sett på diskrete modeller for blant annet CO\(_2\)-utslipp. Da tok vi utgangspunkt i følgende diskrete modell:
som er det samme som:
Vi ser at vi på venstresiden har forskjellen i utslipp fra en tid til en annen. Det kan vi skrive som \(\Delta u\), som jo er et uttrykk for (gjennomsnittlig) endring. Da får vi:
Dersom vi ønsker den momentane endringen istedenfor den gjennomsnittlige, kan vi uttrykke dette ved den deriverte:
Dermed har vi en kontinuerlig modell som beskriver den momentane endringen i et system! Vi kan da bruke denne modellen til å finne tilstanden til systemet ved ethvert lite tidssteg dt. Da kan vi bruke Eulers metode, som ovenfor. Her er et eksempel på hvordan vi kan løse modellen ovenfor:
import matplotlib.pyplot as plt
tid_slutt = 2100 # År etter 1825
tid = 1825 # Startår
a = 0.045 # Utslippsrate
u = 0.01 # Utslipp i 1825 (tonn per innbygger)
dt = 1E-5
utslipp = [u]
årstall = [tid]
while tid <= tid_slutt:
uder = a*u # Differensiallikningen
u = u + uder*dt # Eulers metode
tid = tid + dt # Går fram ett tidssteg
utslipp.append(u)
årstall.append(tid)
plt.plot(årstall, utslipp)
plt.title("Utslipp av CO$_2$ i Norge")
plt.xlabel("År")
plt.ylabel("Tonn CO$_2$ per innbygger")
plt.show()
Hvis vi sammenlikner med den diskrete modellen for det samme systemet, kan vi se at de to modellene ikke gir det samme resultatet. Dette er fordi utslippsraten a ikke har samme enhet. I det diskrete systemet representerte a vekstraten per uke, mens i den kontinuerlige modellen gir a den momentane vekstraten. De to parameterne må derfor bestemmes uavhengig av hverandre.
Underveisoppgave
Følgende modell kan beskrive størrelsen til en populasjon med bæreevnen b, vekstrate a og populasjonsstørrelse N:
Tida t er gitt i timer.
Lag et program som simulerer utviklingen i en bakteriekoloni med N(t\(_0\)) = 100, a = 0,005 og en bæreevne på 120 000 bakterier.
Hvilke forutsetninger og begrensninger ligger til grunn for modellen?
Når vil populasjonen nå bæreevnen?
Lag programmet med arrayer istedenfor lister. Husk at append ikke fungerer med arrayer, og at vi må forhåndsdefinere størrelsen til disse arrayene på forhånd. Størrelsen på disse bør være gitt ved (forklar hvorfor!):
N = int((tid_slutt-tid_start)/dt) + 1 # Antall punkter
Løsningsforslag (deloppgave 4)
import matplotlib.pyplot as plt
import numpy as np
a = 0.045 # Utslippsrate
u0 = 0.01 # Utslipp i 1825 (tonn per innbygger)
tid_slutt = 2100 # År etter 1825
tid_start = 1825 # Startår
dt = 1E-3 # tidssteg
N = int((tid_slutt-tid_start)/dt) + 1 # Antall punkter
u = np.zeros(N)
t = np.zeros(N)
u[0] = u0
t[0] = tid_start
for i in range(N-1):
uder = a*u[i] # Differensiallikningen
u[i+1] = u[i] + uder*dt # Eulers metode
t[i+1] = t[i] + dt # Går fram ett tidssteg
plt.plot(t, u)
plt.title("Utslipp av CO$_2$ i Norge")
plt.xlabel("År")
plt.ylabel("Tonn CO$_2$ per innbygger")
plt.show()
Newtons 2. lov#
Vi så at Newtons 2. lov gir oss differensiallikninger når vi omformulerer akselerasjonen som den deriverte av farten.
Vi kan altså få både en førsteordens og en andreordens difflikning ut av Newtons 2. lov. Det betyr også at vi kan finne akselerasjon og fart ut fra strekning ved å derivere, og motsatt finne fart og posisjon fra akselerasjon ved å integrere (løse difflikningen). Vi kan med andre ord skrive Newtons 2. lov slik:
La oss løse Newtons 2. lov for det enkle systemet vi så på tidligere: Når vi slipper ei kule fra en viss høyde, kan vi regne med at to krefter virker på kula, nemlig tyngdekraften og luftmotstand. For en lett ball som faller, kan vi si at summen av krefter er summen av gravitasjonskraften og luftmotstanden i \(y\)-retning på en graf. Vi definerer positiv retning oppover \(y\)-aksen.
Vi ser at vi nå har den deriverte av hastigheten, som er akselerasjonen. For å kunne løse denne likningen, må vi finne ut hva summen av kreftene er. Du husker kanskje at \(G = mg\), og at en mulig modell for luftmotstand er \(L = -kv\), der \(k\) er en konstant som blant annet er avhengig av legemets form. Da kan vi skrive uttrykket slik:
Dette er en førsteordens difflikning der vi har en funksjon \(v(t)\), og dens deriverte \(a(t)\) (det vil si \(v'(t)\)). La oss løse denne numerisk:
import matplotlib.pyplot as plt
import numpy as np
# Konstanter
k = 0.1 # Luftmotstand
g = 9.81 # Tyngdeakselerasjon i m/s^2
m = 1 # Masse i kg
v0 = 0 # Starthastighet i m/s
s0 = 0 # Startposisjon i m
# Tidsvariabler
dt = 1E-5 # Tidsintervall i s
tid_start = 0
tid_slutt = 10
N = int((tid_slutt-tid_start)/dt) + 1 # Antall punkter
# Arrayer
t = np.zeros(N)
a = np.zeros(N)
v = np.zeros(N)
s = np.zeros(N)
# Startverdier
t[0] = tid_start
v[0] = v0
s[0] = s0
for i in range(N-1):
a[i] = - g - k*v[i]/m
v[i+1] = v[i] + a[i]*dt
s[i+1] = s[i] + v[i]*dt + 0.5*a[i]*dt**2
t[i+1] = t[i] + dt
plt.plot(t,s)
plt.title('Fallende kule')
plt.xlabel('Tid (s)')
plt.ylabel('strekning (m)')
plt.show()
Underveisoppgave
Studer programmet ovenfor og forklar hvordan det fungerer linje for linje. Test programmet med andre modeller for luftmotstand. Vurder de ulike modellene opp mot hverandre.
ODE-løsere#
La oss se nå se på hvordan vi kan løse differensiallikninger med funksjoner fra scipy-biblioteket. De differensiallikningene vi har sett på, kaller vi for ODE-er (ordinary differential equations) for å skille dem fra PDE-er (partial differential equations). Vi skal ikke se på partielle differensiallikninger her, men mange av prinsippene for å løse dem er like som for ODE-er. Vi kommer til å bruke ODE som forkortelse videre.
Vi starter med enkle differensiallikninger for å illustrere de grunnleggende prinsippene. En enkel differensiallikning vi kan begynne med, er:
Denne kan vi ganske enkelt definere som en funksjon og løse med en annen funksjon som heter solve_ivp (“solve initial value problem”).
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def dy_dt(t, y):
return t - y
a = 0
b = 4
t = np.linspace(a,b,1000)
y0 = 1
y_int = solve_ivp(dy_dt, [a,b], [1], t_eval=t)
plt.xlabel("t")
plt.ylabel("y")
plt.plot(y_int.t, y_int.y[0])
plt.show()
Vi definerer her et sett med t-verdier slik at vi får et intervall å integrere over. Dernest er det viktig å ha en startbetingelse for \(y\) (husk: vi trenger den forrige \(y\)-verdien for å finne den neste). Funksjonen solve_ivp er en generell løser for differensiallikninger og tar som første parameter en funksjon av typen \(f(t,y)\) – legg merke til rekkefølgen på parameterne. Deretter legger vi inn tidsintervallet som vi skal integrere over. Dette legges inn som en liste.
Etter tidsintervallet legges startbetingelsene inn. Siden solve_ivp er en løser som kan løse systemer av differensiallikninger, må vi lage startbetingelsen som en liste. Deretter har vi en del valgfrie parametre. Vi har brukt parameteren t_eval her fordi den sier hvilke \(t\)-verdier vi skal regne ut \(y\)-verdier for. Hvis vi ikke gjør dette, får vi integralet kun evaluert i noen få punkter. Det kan være greit hvis vi for eksempel bare ønsker sluttverdien, men ikke hvis vi ønsker å plotte resultatene.
Dersom du prøver å printe ut resultatene fra solve_ivp, får du mye ulik informasjon. Derfor henter vi ut spesifikke t- og _y-_verdier ved å skrive y_int2.t, som henter ut tidsverdiene, og y_int2.y[0], som henter ut y-verdiene. Legg merke til at \(y\) kan inneholde flere elementer ettersom vi kan løse systemer av differensiallikninger. Her må vi eksplisitt be om det første elementet (element 0 med Python-språk), selv om arrayen ikke inneholder flere \(y\)-verdier.
Underveisoppgave
Løs differensiallikningen \(f'(t) = \cos(t)\) med \(f(t_0) = 0\) med solve_ivp fra \(t=0\) til \(t=\pi\). Plott den analytiske løsningen \(f(t) = \sin(t)\) i samme koordinatsystem for å sammenlikne.
Algoritmene som brukes i slike biblioteker, er ofte sammensatte algoritmer som benytter seg av flere prinsipper enn en enkelt metode. Som standard benytter solve_ivp seg av en blanding av Runge-Kutta 4 og Runge-Kutta 5, kalt RK45. Dette er gode “go to”-metoder.
Oppgaver#
Oppgave 1
Utled Eulers metode og skriv et program for hånd som implementerer denne metoden.
Oppgave 2
Løs differensiallikningene nedenfor med Eulers metode (og analytisk hvis du får til). Bruk initialbetingelse \(f(0) = 0\), og velg løsningsområde/løsningsintervall selv.
\(f'(x) = 1\)
\(f'(x) = x\)
\(f'(x) = 4 + 3\cdot f(x)\)
\(y' + 4y + 3x = 0\)
\(xy' + y = 3x^2 + 4x\)
\(f'(x) = x + y^2\)
Oppgave 3
Hvorfor må vi ha initialbetingelser for å løse en differensiallikning numerisk?
Oppgave 4
Det er alltids fint å ha definert funksjoner i et program – da kan du gjenbruke dem i andre programmer! a. Skriv et program der Eulers metode er definert som en funksjon. Funksjonen skal ta inn som argument yder, som antas å være en funksjon av to variable \(x\) og \(y\). Funksjonen skal også ta inn en initialbetingelse y0, startverdien og sluttverdien til x og steglengden dt. Funksjonene skal så returnere resultatene, altså lister eller arrayer med verdier av y og x. Tips: For at en funksjon skal returnere to verdier, kan du skrive return x,y.
b. Bruk funksjonen til å løse likninga \(y'(x) = -\dfrac{xy}{x^2+1}\) der \(y(0) = 1\) og \(x \in [0,7]\). Tips: Funksjonen returnerer to verdier:y og x. Du kan hente to eller flere verdier fra en funksjon, ved å gjøre slik:
x, y = # Her kaller du på funksjonen din med riktige parametre
Da vil variabelen x settes til å være lik den returnerte verdien for x, og variabelen y settes til å være lik den returnerte verdien for y.
c. Plott resultatet ditt sammen med den analytiske løsninga for å se om funksjonen din er implementert riktig:
Oppgave 5
Programmet under løser likninga \(f'(x) = 2f(x)\) der \(f(0) = 2\) for \(x\in[0,1]\):
import matplotlib.pyplot as plt
import numpy as np
N = 1000
a = 0
b = 1
x = np.linspace(a,b,N)
dx = x[1]-x[0]
fder = np.zeros(N)
f = np.zeros(N)
f[0] = 2
for i in range(N-1):
fder[i] = 2*f[i]
f[i+1] = f[i] + dx*fder[i]
plot(x,f)
xlabel("x")
ylabel("f(x)")
title("Plott av f(x)")
show()
Oppgave 6
Kari begynner med å ha en årsinntekt på \(I(0) = 500 000\) kroner og har nylig oppretta en sparekonto. Hun har derfor \(S(0) = 0\) kroner i sparekontoen.
Vi velger å prøve ut følgende modeller for forandringen i inntekten hennes \(I'(t)\) og sparekontoen hennes \(S'(t)\):
Lag et program som finner ut hvordan \(S(t)\) vil utvikle seg i løpet av 15 år. Visualiser resultatet.
Oppgave 7
Lag et program som regner ut posisjon, hastighet og akselerasjon til et legeme som faller. Bruk luftmotstandskoeffisienten 0.1 og masse 1 kg.
a. Benytt først en enkel modell uten luftmotstand, det vil si kun fritt fall.
b. Legg inn luftmotstand og plott i samme koordinatsystem som beregninga uten luftmotstand. Sammenlikn resultatene og kommenter svaret.
c. Varier systematisk massen og luftmotstandskoeffisienten og se hva dette gjør med resultatene.
Oppgave 8
Lag et program som modellerer en fiskepopulasjon. Du kan selv bestemme modellen du bruker og initialbetingelsene. Forklar hvorfor du brukte akkurat denne modellen.
Oppgave 9
Newtons avkjølingslov forteller oss at temperaturen \(T(t)\) etter ei tid \(t\) til et objekt i en omgivelse forandrer seg på følgende måte:
der \(k\) er er en konstant som forteller hvor raskt temperaturen forandrer seg og \(T_o\) representerer temperaturen til omgivelsen objektet befinner seg i.
a. Et varmeanlegg i en bil påvirker innetemperaturen \(T(t)\) i bilen. Temperaturen inni bilen kan uttrykkes ved:
der \(T(0) = 7^oC\).
Lag et program som bruker Eulers metode til å finne funksjonsverdiene til \(T(t)\) for \(t \in [0, 30]\). Plott resultatene.
b. Som du kanskje så i deloppgave a, tar det fryktelig lang tid før temperaturen i bilen når vanlig romtemperatur som ligger på \(20^oC\). Du finner derfor en som prøver å fikse på varmeanlegget billigst mulig. Uheldigvis greier personen å gjøre varmeanlegget verre. Etter noen målinger finner du at temperaturen forandrer seg på følgende måte:
der \(T_o(t) = 5\cdot\cos\left(\frac{\pi^2}{15}\cdot t\right) + 19\)
Utvid programmet ditt slik at det plotter bilens temperatur med det (enda mer) ødelagte varmeanlegget. Sammenlikn utfallene.
Oppgave 10
Arne skal på familiemiddag og har lovet familien sin å ta med sin smakfulle tilberedelse av en kalkun. Fem timer før middagen kommer han plutselig på middagsavtalen. Uheldigvis bruker han vanligvis fire timer på å kjøre til familien hvis han kjører omtrent 80 km i timen.
Problemet er at kalkunen må være i en ovn på 180\(^o\)C i tre timer for å bli spiseklar. Han veit at under panseret på bilen kan det bli varmere enn dette. Han bestemmer seg derfor for å prøve å tilberede kalkunen under panseret mens han kjører.
Arne har lagd følgende modell for temperaturen under panseret etter ei tid \(t\):
der \(v(t)\) er farten til bilen og \(p\) er et tilfeldig flyttall mellom 1 og 2. Verdien til \(p\) forandrer seg for hver tid \(t\)
Fra Newtons avkjølingslov vil temperaturen \(T(t)\) til kalkunen forandre seg på følgende måte:
Akselerasjonen \(a(t)\) til bilen etter ei tid \(t\) vil følge modellen
Lag et program som finner kalkunens temperatur etter fire timer under panseret. Plott temperaturen og forklar ut fra plottet om kalkunen vil være spiseklar eller ikke etter kjøreturen (I episoden Food Fables fra Mythbusters kan du se at dette faktisk er mulig!).