Logistisk vekst

Vi kopierer funksjonen fra regresjonsanalyseverktøyet til grafikkfeltet. Så bruker vi CAS til å analysere funksjonen. Her må vi finne det punktet der grafen er brattest. Det vil være i vendepunktet til grafen til $$ g$$.

I linje 3 og 4 kontrollerer vi at nullpunktet til den dobbeltderiverte faktisk er et vendepunkt ved å sjekke om fortegnet skifter ved nullpunktet, noe det gjør. Ørretbestanden vokste raskest nesten 5 år etter at kalkingen startet, og da vokste bestanden med 2 740 ørreter per år.

I oppgave FM-23 viser vi at vi alltid har ett vendepunkt for en logistisk funksjon, og at funksjonen er brattest der, så vi trenger ikke å sjekke dette flere ganger.

Når $$ x$$ betyr antall år etter 1970, får vi at $$ x=0$$ betyr året 1970, $$ x=1$$ betyr 1971 og så videre. Vi legger derfor tallene fra og med 0 til og med 10 i én kolonne i regnearket i GeoGebra og antall løvetannplanter i neste kolonne. Så velger vi "Regresjonsanalyse" og regresjonsmodellen logistisk.

Åpne bilde i et nytt vindu

CC BY-SA

Bilde: Bjarne Skurdal

Bildet viser at en logistisk modell passer bra med tallene. Den modellen som passer best, er

$$ n\left(x\right)=\frac{1675,39}{1+46,3 e^{-0,58x}}$$

Når $$ x\to \infty $$, vil $$ e^{-0,58x}\to 0$$ og nevneren gå mot 1. Funksjonen går derfor mot tallet i telleren.

Antall løvetannplanter på lang sikt, eller det vi kaller bæreevnen, vil derfor være 1 675 etter modellen.

Vi kopierer funksjonen fra regresjonsanalyseverktøyet til grafikkfeltet. Så bruker vi CAS til å analysere funksjonen.

Etter modellen økte antallet løvetannplanter mest i det sjette året etter 1970, det vil si i 1976. Da økte antallet løvetannplanter med 244 planter per år.

Forslag til kode:

1        # importerer nødvendige bibliotek
2from scipy.optimize import curve_fit
3import numpy as np
4import matplotlib.pyplot as plt
5
6        # lager funksjonen som beskriver modellen
7def modell(x,B,a,k):
8  return B/(1+a*np.e**(-k*x))
9  
10        # legger inn måledataene i lister
11x_verdier = [0,1,2,3,4,5,6,7,8,9,10]
12y_verdier = [3,15,80,201,300,501,731,915,1131,1350,1490]
13
14        # bruker metoden curve_fit og legger resultatene i to lister
15konstanter,kovarians = curve_fit(modell,x_verdier,y_verdier, p0 = [max(y_verdier),1,1])
16        # henter ut konstantene fra lista konstanter
17B, a, k = konstanter
18
19        # lager utskrift av funksjonen
20print(f"Funksjonen blir n(x) = {B:.0f}/(1+{a:.2f}e^(-{k:.2f}x)).")
21
22        # plotter data
23plt.plot(x_verdier,y_verdier,'.', label = "Målinger")
24 
25        # plotter modell
26x_array = np.linspace(min(x_verdier),max(x_verdier),300)
27y_array = modell(x_array,B,a,k)
28plt.plot(x_array,y_array,"brown", label = "Modell")
29plt.grid(True)
30plt.legend(bbox_to_anchor=(0.6,1))
31plt.xlabel("$x$, antall år etter 1970") # tittel på x-aksen
32plt.ylabel("$y$, antall løvetannplanter")
33plt.show()

Koden gir følgende utskrift:

"Funksjonen blir n(x) = 1675/(1+46.30e^(-0.58x))."

Kommentar til koden: I linje 12 har vi lagt til p0 = [max(y_verdier),1,1] i kommandoen curve_fit. Hvis vi ikke gjør det, feiler regresjonen. Tallene max(y_verdier), 1 og 1 er startverdier for de tre variablene B, a og k, som programmet bruker i de numeriske beregningene for å komme fram til de endelige verdiene for de tre variablene. Jo riktigere startverdiene er, jo større er sjansen for at de numeriske beregningene blir vellykket.

max(y_verdier) er et godt valg for startverdien til $$ B$$ fordi den plukker ut den største $$ y$$-verdien i målingene, og vi vet at i en logistisk modell er $$ B$$ den øvre grenseverdien for funksjonen. For de to andre parametrene a og k setter vi bare startverdien lik 1, som er den automatiske startverdien for parametre ved bruk av curve_fit.

Vi skriver inn tallene i tabellen inn i regnearkdelen i GeoGebra og regner ut verdiene for $$ t$$. Så velger vi verktøyet "Regresjonsanalyse" og regresjonsmodellen "Logistisk".

Vi ser at grafen som kommer opp, passer godt til punktene. Den logistiske funksjonen som passer best med tallene i tabellen, er

$$ N\left(t\right)=\frac{4 988}{1+9,10 e^{-0,203t}}$$

Vi kopierer grafen (inkludert punktene) fra regresjonsverktøyet over i grafikkfeltet.

Her må vi basere oss på modellen. Når $$ t\to \infty $$, vil $$ e^{-0,203t}\to 0$$, nevneren i funksjonen, gå mot 1, og antallet hjort vil gå mot tallet i telleren. Det er dette vi kaller bæreevnen for hjortebestanden.

Bæreevnen er 4 988 hjort etter modellen, noe vi runder av til 5 000.

Vi må finne vendepunktet til funksjonen, og vi bruker CAS.

Hjortebestanden vokste mest omtrent 11 år etter 1973, det vil si i 1984, og da vokste bestanden med 253 dyr per år.

Vi får at hjortebestanden var på halvparten av bæreevnen da den vokste raskest.

Vil det være slik for alle logistiske modeller, tror du?

Vi får at grafen til $$ f$$ vokser mest når $$ x=\frac{\ln a}{k}$$, og da vokser den med $$ \frac{1}{4}Bk$$ per enhet på $$ x$$-aksen. Det betyr at for alle logistiske modeller kan vi finne når veksten er størst, og hvor stor den er, ved å bruke resultatene i linje 2 og linje 4.

For å sjekke at nullpunktet til den dobbeltderiverte faktisk er et vendepunkt, kan vi ikke løse ulikheten $$ f^{\text{'}\text{'}}\left(x\right)>0$$ slik vi vanligvis ville ha gjort fordi GeoGebra ikke vet at konstantene $$ a, k$$ og $$ B$$ er positive slik forutsetningen er. Det vi kan gjøre, er å bruke den trippelderiverte funksjonen. Et vendepunkt er et ekstremalpunkt for den deriverte. Da kan vi bruke dobbeltderiverttesten på den deriverte (altså den trippelderiverte til $$ f$$) for å sjekke om det er et toppunkt eller et bunnpunkt for den deriverte. Vi ser av linje 3 at den trippelderiverte er negativ, som betyr at grafen til den deriverte vender den hule siden ned. Dette viser at den deriverte har et toppunkt der, og vi har vist at veksten er størst der.

Siden $$ B$$ er den verdien funksjonen går mot når $$ x\to \infty $$, det vi kaller bæreevnen, er antallet individer når funksjonen vokser raskest, lik halve bæreevnen.

Vi starter med å derivere funksjonen to ganger.

$$ \begin{array}{rcl}f\left(x\right) & =& \frac{B}{1+a·e^{-kx}}\\ f^{\text{'}}\left(x\right) & =& \frac{\left(1+a·e^{-kx}\right)·0-B·a·\left(-k\right)e^{-kx}}{{\left(1+a·e^{-kx}\right)}^2}\\ & =& \frac{akB{e}^{-kx}}{{\left(1+a·e^{-kx}\right)}^2}\\ f^{\text{'}\text{'}}\left(x\right) & =& \frac{{\left(1+a·e^{-kx}\right)}^2·akB·\left(-k\right)e^{-kx}-akB{e}^{-kx}·2\left(1+a·e^{-kx}\right)·a\left(-k\right)e^{-kx}}{{\left(1+a·e^{-kx}\right)}^4}\\ & =& \frac{-a{k}^{2}B{e}^{-kx}{\left(1+a·e^{-kx}\right)}^2+2a^2{k}^{2}B{e}^{-kx}\left(1+a·e^{-kx}\right)e^{-kx}}{{\left(1+a·e^{-kx}\right)}^4}\\ & =& \frac{a{k}^{2}B{e}^{-kx}\cancel{\left(1+a·e^{-kx}\right)}\left(-\left(1+a·e^{-kx}\right)+2a{e}^{-kx}\right)}{{\left(1+a·e^{-kx}\right)}^{\cancel{4} 3}}\\ & =& \frac{a{k}^{2}B{e}^{-kx}\left(a{e}^{-kx}-1\right)}{{\left(1+a·e^{-kx}\right)}^3}\end{array}$$

Kommentar: I tredje linje i utregningen av $$ f^{\text{'}\text{'}}$$ har vi satt uttrykket $$ a{k}^{2}B{e}^{-kx}\left(1+a·e^{-kx}\right)$$ utenfor parentes.

Når funksjonen er mål på en bestand, er konstantene $$ a, k$$ og $$ B$$ positive. Det betyr for eksempel at nevneren i den dobbeltderiverte alltid er positiv.

Så må vi finne nullpunktene til den dobbeltderiverte.

$$ \begin{array}{rcl}{f}^{\text{'}\text{'}}\left(x\right) & =& 0\\ \frac{a{k}^{2}B{e}^{-kx}\left(a{e}^{-kx}-1\right)}{{\left(1+a·e^{-kx}\right)}^3} & =& 0\\ a{k}^{2}B{e}^{-kx}\left(a{e}^{-kx}-1\right) & =& 0\\ a{e}^{-kx}-1 & =& 0\\ a{e}^{-kx} & =& 1\\ \ln \left(e^{-kx}\right) & =& \ln \frac{1}{a}=\ln a^{-1}=-\ln a\\ -kx & =& -\ln a\\ x & =& \frac{\ln a}{k}\end{array}$$

I løsningen har vi enkelt kvittet oss med nevneren og uttrykket foran parentesen i telleren siden begge disse ikke kan være null. Vi må nå sjekke at løsningen faktisk er $$ x$$-koordinaten til et vendepunkt. Det er slitsomt å skulle derivere én gang til for hånd slik vi gjorde enkelt med CAS. Vi prøver å løse ulikheten $$ f^{\text{'}\text{'}}\left(0\right)>0$$ slik vi vanligvis gjør i stedet.

$$ $ \begin{array}{rcl}{f}^{\text{'}\text{'}}\left(x\right) & >& 0\\ \frac{a{k}^{2}B{e}^{-kx}\left(a{e}^{-kx}-1\right)}{{\left(1+a·e^{-kx}\right)}^3} & >& 0\\ a{k}^{2}B{e}^{-kx}\left(a{e}^{-kx}-1\right) & >& 0\\ a{e}^{-kx}-1 & >& 0\\ a{e}^{-kx} & >& 1\\ \ln \left(e^{-kx}\right) & >& \ln \frac{1}{a}=\ln a^{-1}=-\ln a\\ -kx & >& -\ln a\\ x & <& \frac{\ln a}{k}\\ & & \end{array}$$$

Kommentar til løsningen: Siden de tre konstantene $$ a, k$$ og $$ B$$ alle er større enn null, vet vi at vi kan dele med dem uten å måtte snu ulikhetstegnet i overgangen mellom tredje og fjerde linje.

Vi får altså at grafen vender sin hule side opp på venstresiden av nullpunktet og motsatt på høyre side, som vi forventet. Da vet vi at løsningen er $$ x$$-koordinaten til et vendepunkt, og videre at grafen vil stige raskest her.

Til slutt må vi regne ut $$ f^{\text{'}}\left(\frac{\ln a}{k}\right)$$ og $$ f\left(\frac{\ln a}{k}\right)$$.

$$ \begin{array}{rcl}{f}^{\text{'}}\left(\frac{\ln a}{k}\right) & =& \frac{akB{e}^{-k{\displaystyle \frac{\ln a}{k}}}}{{\left(1+a{e}^{-k{\displaystyle \frac{\ln a}{k}}}\right)}^2}=\frac{akB{e}^{\ln \left(a^{-1}\right)}}{{\left(1+a{e}^{\ln \left(a^{-1}\right)}\right)}^2}\\ & =& \frac{akB{a}^{-1}}{{\left(1+a{a}^{-1}\right)}^2}=\frac{kB}{{\left(1+1\right)}^2}\\ & =& \frac{kB}{4}\end{array}$$

$$ \begin{array}{rcl}f\left(\frac{\ln a}{k}\right) & =& \frac{B}{1+a{e}^{-k{\displaystyle \frac{\ln a}{k}}}}\\ & =& \frac{B}{1+a{a}^{-1}}\\ & =& \frac{B}{2}\end{array}$$