Modul 3

import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

Modul 3#

Populationsdynamik og Bæredygtighed#

Det er helt centralt i forbindelse med politiske beslutninger for samfundet at kunne forudsige, hvordan en given population udvikler sig fremover. Det er selvsagt vigtigt at se på et lands befolkning og kunne planlægge dagsinstitutioner, skoler, ældrepleje etc., og vurdere den økonomiske bæredygtighed for den demografiske profil. Præcis samme overvejelser er meningsfulde for dyrearter. Det er for eksempel nødvendigt at kunne beskrive dynamikken for en fiskepopulation før man går i gang med at uddele fiskekvoter. Det er samme problemstillinger blot i meget forskellig kontekst.

I denne notebook vil vi se på modellering af populationer ved brug af matricer og linear algebra. De grundlæggende ideer går tilbage til Patrick Leslie i 1945, som forelog at studere populationer ved at inddele i aldersklasser og trække på viden om klassernes individuelle fertilitets- og overlevelsesrater. Som eksempel så han på den brune rottes population. Samlet giver det et diskret, dynamisk system, som kan analyseres med egenværdier og -vektorer.

Nedenfor vil starte med små og simple modeller med udgangspunkt i kaninpopulationer og Fibonacci-tallene. Modellerne udvides, fertilitetsrater og overlevelsesrater introduceres, og vi ser på, hvordan konkurrerende arter kan modelleres. Endeligt når vi til en model for den danske befolknings udvikling med betingelser og rater baseret på data fra Danmarks statistik.

Undervejs vil vi fremhæve væsentlige emner fra lineær algebra, matematisk modellering med Lesliematricer, egenværdier og -vektorer, iteration og konvergens, og vi vil tage nogle små sidespring for at belyse beregningselementer, herunder se på kontrolstrukturer og kompleksitet.

Program#

Opvarmning: Fibonacci-tallene genbesøgt
Getting real - fødsel og overlevelse
Konkurrerende arter
Den danske befolkning

Opsætning#

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
try:
    from ipywidgets import interact, FloatSlider, Layout
except ImportError:
    import micropip
    await micropip.install("ipywidgets")
    from ipywidgets import interact, FloatSlider, Layout

import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

from matplotlib.patches import Circle
from IPython.display import clear_output

1. Opvarmning: Fibonacci-tallene genbesøgt#

Fibonacci var en italiensk matematiker i det 13. århundrede; han er også kendt som Leonardo Figlio di Bonacci. Gennem rejser til nord-Afrika stiftede han bekendtskab med det hindu-arabiske talsystem (tallene 0-9, decimaler), som han i sin bog Liber Abaci fra 1203 bragte til Europa. De nye tal afløste romertallene og bragte en revolution for handel og bankvæsenet; også dengang var matematik en kilde til magt og rigdom.

Fibonacci havde et væsentligt fokus på algoritmer, ligninger og praktisk problem-løsning, og det er helt rimeligt at se ham som en pioner inden for computational thinking. Meget længe inden computerens fødsel.

I Liber Abaci gennemfører Fibonacci et tankeeksperiment om en kaninpopulations udvikling i isolerede omgivelser. Udgangspunktet er et kaninpar og den matematiske model, at ethvert voksent par giver afkom i form af et nyt, ungt par hver måned. Unger bliver voksne på en måned, og der er ingen dødelighed.

Vi vælger at tælle hunkaniner fremfor par. Vi starter med én voksen hunkanin, og i en kontekst af en populationsdynamik ser vi på en aldersopdelt populationsmodel, hvor der skelnes mellem to aldersklasser:

\(\boldsymbol{x}_k[0]\): antal unge hunkaniner efter \(k\) måneder,
\(\boldsymbol{x}_k[1]\): antal voksne hunkaniner efter \(k\) måneder.

Til \(k=0\) har vi \(\boldsymbol{x}_0[0] = 0\) og \(\boldsymbol{x}_0[1] = 1.\) Dynamikken beskrives ved

\[\begin{align*} \boldsymbol{x}_{k+1}[0] &= \boldsymbol{x}_{k}[1] &\qquad \text{(hver voksen kanin giver anledning til en unge i næste måned)},\\ \boldsymbol{x}_{k+1}[1] &= \boldsymbol{x}_{k}[0] + \boldsymbol{x}_{k}[1] & \qquad \text{(hver unge giver anledning til en voksen i næste måned, og de nuværende voksne overlever)}. \end{align*}\]

Hvis vi introducerer matricen

\[\begin{split} A = \begin{bmatrix}0 & 1\\[4pt]1 & 1\end{bmatrix} \end{split}\]

så er populationen beskrevet ved det diskrete dynamiske system

\[\begin{align*} \boldsymbol{x}_{k+1} = A \boldsymbol{x}_k, \quad k = 0,1,2,\ldots\;\;, \quad \boldsymbol{x}_0 = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \end{align*}\]

Dermed kan vi følge, hvordan antallet af unge og voksne udvikler sig måned for måned.

Vi kan beskrive udviklingen ved at iterere matricen \(A\) så

\[ \boldsymbol{x}_1 = A \boldsymbol{x}_0, \quad \boldsymbol{x}_2 = A \boldsymbol{x}_1 = A^2 \boldsymbol{x}_0, \quad \dots, \quad \boldsymbol{x}_{k+1} = A \boldsymbol{x}_k = A^{k+1} \boldsymbol{x}_0. \]

Ser vi kun på det antallet af voksne hunkaniner \(F_k = \boldsymbol{x}_k[1]\), genfinder vi rekursionen

\[ F_{k+1} = F_k + F_{k-1}, \quad k = 1,2,3,\ldots\;\;, \]

som med begyndelsen \(F_0 = F_1 = 1\) lige præcis giver Fibonacci-tallene \(1,1,2,3,5,8,\ldots\)

Matricen \(A\) kaldes derfor ofte Fibonacci-matricen, da hver iteration svarer til én måned i modellen, og den voksne bestand følger Fibonacci-følgen.

Læg mærke til, hvordan første og anden række i \(A\) beskriver henholdsvis næste generations unger og voksne. Fortolk ligeledes søjlerne i \(A\) i en kontekst af kaninpopulationer.

Iteration og for-løkker i Python#

Løkker (loops) er en såkaldt kontrolstruktur, som tillader os at løbe en gennem en sekvens af elementer og for hvert element udføre en given handling.

Før vi går i gang med at implementere Fibonacci-matricen og iterere i Python, laver vi en simpel opvarmningsøvelse.

Læs koden nedenfor, og svar på spørgsmålet inden du kører den i Python.

Hvad forventer du at der bliver printet? Kør koden og sammenlign med din forudsigelse.

k = 1
for i in range(5):
    k = k + i
    print("i =", i, ", k =", k)

Hvad sker der hvis vi ændrer range(5) til range(1,6)?

Færdiggør kodecellen nedenfor, så der benyttes et for-loop til at beregne og printe de første \(10\) Fibonacci-tal. NB: I koden bruger vi variablerne F0 og F1 til løbende at beregne og overskrive de to seneste Fibonacci-tal. Det kan virke en smule forvirrende, fordi man kunne tro, at F0 altid repræsenterer det første Fibonacci-tal og F1 det næste. Her anvender vi dem dog kun som “pladsholdere” for de nyeste værdier i sekvensen - ikke som faste indeks.

F0 = "INDSÆT KODE HER"
F1 = "INDSÆT KODE HER"

for i in range("INDSÆT KODE HER"):
    print("INDSÆT KODE HER")
    F0, F1 = "INDSÆT KODE HER" , "INDSÆT KODE HER"

Matrix-iteration#

Vi har set, at kaninpopulationen kan beskrives med

\[\begin{split} \boldsymbol{x}_{k+1} = A \boldsymbol{x}_k = A^{k+1}\boldsymbol{x}_0, \quad A = \begin{bmatrix}0 & 1\\ 1 & 1\end{bmatrix}. \end{split}\]

hvor startpopulationen er

\[\begin{split} \boldsymbol{x}_0 = \begin{bmatrix}0\\ 1\end{bmatrix} \quad \text{(0 unge, 1 voksen).} \end{split}\]

Færdiggør kodecellen nedenfor, så populationen for de første 10 måneder udregnes ved at gentage \(\boldsymbol{x}_{k+1} = A \boldsymbol{x}_k\). Der printes antal unge, voksne og total for hvert måned.

A = "INDSÆT KODE HER"

startpopulation = "INDSÆT KODE HER"

antal_md = "INDSÆT KODE HER"

print("Måned | Unge | Voksne | Total")
population = startpopulation.copy()
for md in range(antal_md):
    print(f"{md:5d} | {population[0]:4d} | {population[1]:6d} | {population.sum():3d}")
    # Opdater populationen
    population = "INDSÆT KODE HER"

I følgende kodecelle ønsker vi at vise udviklingen af voksne kaniner i et plot. Vi bruger listen voksne_list til løbende at gemme antallet i hver iteration.

Færdiggør kodecellen nedenfor, så antallet af voksne kaniner afbildes for de første 10 måneder ved hjælp af matrix-iteration.

A = "INDSÆT KODE HER"

startpopulation = "INDSÆT KODE HER"

antal_md = "INDSÆT KODE HER"

# tome liste til antal voksne kaniner
voksne_list = []

population = startpopulation.copy()
for md in range(antal_md):
    # tilføj antal voksne kaniner til listen
    voksne_list.append(population[1])

    # Opdater populationen
    population = "INDSÆT KODE HER"
    

# Plot antal voksne kaniner
plt.figure(figsize=(7,4))
plt.plot(range(antal_md), voksne_list, 'o-', color='green')
plt.xlabel("Måned")
plt.ylabel("Antal voksne kaniner")
plt.title("Udvikling af voksne kaniner (måned 0-9)")
plt.grid(True)
plt.show()

Udnyt at \(\boldsymbol{x}_{k} = A^{k}\boldsymbol{x}_0\) til at finde populationen i måned 9 uden at iterere måned for måned.
\(A^{k}\) kan laves med np.linalg.matrix_power(A,k).

# INDSÆT KODE HER

Sammenlign dit resultat med opgaven, hvor du bestemte \(\boldsymbol{x}_{k+1}\) ved at lave gentagne \(\boldsymbol{x}_{k+1}=A\boldsymbol{x}_k\).

Egenværdier, singulærværdier og asymptotisk udvikling#

En egenværdi \(\lambda \in \mathbb{R}\) og tilhørende egenvektor \(\boldsymbol{v} \neq 0\) opfylder egenværdiproblemet

\[ A \boldsymbol{v} = \lambda \boldsymbol{v}. \]

Det betyder, at når vi anvender matricen \(A\) på vektoren \(\boldsymbol{v}\), så bliver resultatet kun en skalering af \(\boldsymbol{v}\) med \(\lambda.\)

Geometrisk er \(A v\) parallel med \(\boldsymbol{v}\), men længden ændres med faktor \(\lambda\).

En egenværdi \(\lambda\) kan findes med determinanten og opfylder ligningen

\[ \det(A - \lambda I) = 0. \]

For Fibonacci-matricen \(A\) betyder det, at

\[\begin{split} \det\begin{pmatrix}-\lambda & 1\\ 1 & 1-\lambda\end{pmatrix} = -\lambda(1-\lambda) - 1 = 0. \end{split}\]

Løs andengradsligningen for \(\lambda\).

Indse, at der findes en positiv løsning, \(\lambda = \lambda_1,\) og en negativ løsning \(\lambda = \lambda_2,\) og at \(\lambda_1 > |\lambda_2|.\)

lambda1 = "INDSÆT"
lambda2 = "INDSÆT"

Vi går nu til en geometrisk fortolkning.

Nedenfor finder du den interaktive funktion enhedsvektor_afbildning(). Det ser omfattende ud, men du behøver ikke forstå hvert trin.

def enhedsvektor_afbildning(A, slider_width='80%', figsize=(6, 6)):
    """
    Interaktiv visualisering af enhedsvektorer x og deres billede y = Ax
    for en 2x2-matrix A.
    """

    # Definer enhedsvektor som funktion af vinkel t
    def x(t): 
        return np.array([np.cos(t), np.sin(t)])  # Punkt på enhedscirklen

    # Definer billedvektor y = A@x(t)
    def y(t): 
        return A @ x(t)

    # Opret punkter for hele cirklen (fra 0 til 2π)
    ts = np.linspace(0, 2*np.pi, 361)

    # Punkter i figuren for at skabe passende aksegrænser
    circle_pts = np.column_stack([x(tt) for tt in ts]).T   # Punkter på enhedscirklen
    image_pts = np.column_stack([y(tt) for tt in ts]).T    # Punkter efter transformationen
    all_pts = np.vstack([circle_pts, image_pts])           
    lim = 1.15 * np.max(np.abs(all_pts))                   

    # Funktion der tegner vektorer for et givent t
    def plot_vectors(t):
        clear_output(wait=True)  
        fig, ax = plt.subplots(figsize=figsize)

        # plot udseende
        ax.set_aspect('equal')   
        ax.grid(True, linestyle='--', linewidth=0.4, alpha=0.5)
        ax.set_xlim(-lim, lim)
        ax.set_ylim(-lim, lim)

        # Tegn enhedscirklen
        ax.plot(circle_pts[:, 0], circle_pts[:, 1], color='lightgrey', linewidth=1.5)
        ax.add_patch(Circle((0, 0), radius=1, edgecolor='lightgrey', facecolor='none', linewidth=1.5))

        # Udvælg punkt på cirklen og dets billede
        v = x(t)     
        Av = y(t)    

        # Hjælpefunktion til at tegne pile for vektorer
        def draw_arrow(ax, vec, color, label, z):
            ax.quiver(0, 0, vec[0], vec[1], angles='xy', scale_units='xy', scale=1,
                      color=color, width=0.012, zorder=z)
            # Skriv koordinater ved spidsen af pilen
            ax.text(vec[0], vec[1], f'  {label}{tuple(np.round(vec, 2))}',
                    va='bottom', ha='left', fontsize=9, color=color)

        # Tegn vektorerne x og Ax
        draw_arrow(ax, v, 'tab:blue', 'x ', 3)
        draw_arrow(ax, Av, 'tab:red', 'Ax ', 4)

        # Titel med aktuel vinkelværdi
        ax.set_title(rf"$t={t:.2f}$")
        plt.show()

    # Slider til at ændre vinklen t interaktivt
    slider = FloatSlider(
        value=0,
        min=0,
        max=2*np.pi,
        step=0.01,
        description='t',
        continuous_update=True,
        readout_format='.2f',
        layout=Layout(width=slider_width)
    )

    # Forbind slideren med visualiseringsfunktionen
    interact(plot_vectors, t=slider)

Kør følgende eksempel.

# eksempel
A = np.array([[0, 1],
              [1, 1]])
enhedsvektor_afbildning(A)

Aflæs ved brug af enhedsvektor_afbildning() to (lineært) uafhængige vektorer \(\boldsymbol{v}=\boldsymbol{v}_1\) og \(\boldsymbol{v}=\boldsymbol{v}_2\), for hvilke \(A \boldsymbol{v}\) er parallel med \(\boldsymbol{v}\).

Disse er egenvektorerne. Bemærk, at hvis \(\boldsymbol{v}\) er en egenvektor, så er \(c\boldsymbol{v}\) det også for ethvert tal \(c\in\mathbb{R}\). Specielt er \(-\boldsymbol{v}\) også en egenvektor med samme egenværdi.

Aflæs tilsvarende de tilhørende egenværdier \(\lambda=\lambda_1\) og \(\lambda=\lambda_2.\)

Definer dem i nedenstående kodecelle ved at bruge np.array().

v1 = "INDSÆT VÆRDIER HER"
lamb1 = "INDSÆT VÆRDIER HER"

v2 = "INDSÆT VÆRDIER HER"
lamb2 = "INDSÆT VÆRDIER HER"

Sammenlign de aflæste egenværdier med de eksakte værdier \(\lambda_1,\lambda_2\) fundet ovenfor.

Egenvektorer hørende til en egenværdi \(\lambda\) kan findes ved at løse ligningen

\[ (A-\lambda I )\boldsymbol{v} = 0. \]

I nedenstående kode bruger vi vores indsigt fra SVD til at finde en enhedsvektor \(\boldsymbol{w}_1,\) som (omtrent) er en egenvektor.

M = A - lambda1 * np.eye(A.shape[0])

# Beregn SVD af M
U, S, Vh = np.linalg.svd(M)

# Vi husker, at de singulære værdier i S er aftagende, så det sidste må afspejle nulrummet for M
w1 = Vh.T[:,1]

Sammenlign den aflæste \(\boldsymbol{v}_1\) med fundne vektor \(\boldsymbol{w}_1\).

print(w1)
print(v1)
print(S)

I Python kan egenværdierne og deres tilhørende egenvektorer bestemmes med np.linalg.eig(A).

lamb, V = np.linalg.eig(A)
print("Vektor med egenværdier:\n",lamb)
print("Matrix med egenvektorer som søjler:\n",V)

Sammenlign egenværdierne med de teoretiske værdier, som du tidligere har bestemt.

Til vores anvendelser ønsker vi egenværdierne i faldende rækkefølge.

# 1. Sortér efter faldende egenværdi
idx = np.argsort(lamb)[::-1]
lamb = lamb[idx]
V = V[:, idx]

# definerer egenværdier og egenvektorer
lamb1 = lamb[0]
lamb2 = lamb[1]
v1 = V[:,0]
v2 = V[:,1]
print("lambda1 = ",lamb1)
print("lambda2 = ",lamb2)
print("v1 = ",v1)
print("v2 = ",v2)

Asymptotisk udvikling#

Vi kan skrive startvektoren som en linearkombination af egenvektorerne \(v_1\) og \(v_2\), så lad

\[ \boldsymbol{x}_0 = c_1 \boldsymbol{v}_1 + c_2 \boldsymbol{v}_2. \]

Når vi anvender matricen \(A\) på \(\boldsymbol{x}_0\), får vi pga. linearitet, at

\[ \boldsymbol{x}_1 = A \boldsymbol{x}_0 = A(c_1 \boldsymbol{v}_1 + c_2 \boldsymbol{v}_2) = c_1 A \boldsymbol{v}_1 + c_2 A \boldsymbol{v}_2 = c_1 \lambda_1 \boldsymbol{v}_1 + c_2 \lambda_2 \boldsymbol{v}_2. \]

Anvendes \(A\) igen, fås

\[ \boldsymbol{x}_2 = A \boldsymbol{x}_1 = A(c_1 \lambda_1 \boldsymbol{v}_1 + c_2 \lambda_2 \boldsymbol{v}_2) = c_1 \lambda_1 A \boldsymbol{v}_1 + c_2 \lambda_2 A \boldsymbol{v}_2 = c_1 \lambda_1^2 \boldsymbol{v}_1 + c_2 \lambda_2^2 \boldsymbol{v}_2. \]

Gentager vi processen, får vi for den \(k.\) iteration, at

\[ \boldsymbol{x}_k = c_1 \lambda_1^k \boldsymbol{v}_1 + c_2 \lambda_2^k \boldsymbol{v}_2. \]

Brug egenvektor-dekompositionen til at beregne \(\boldsymbol{x}_k\) for \(k=0,\dots,K\) via

\[\boldsymbol{x}_k = c_1 \lambda_1^k \boldsymbol{v}_1 + c_2 \lambda_2^k \boldsymbol{v}_2.\]

Du skal altså bestemme \(c_1,c_2\) fra \(\boldsymbol{x}_0\) og illustrer \(\boldsymbol{x}_k[0]\) og \(\boldsymbol{x}_k[1]\) som funktion af \(k\) ved at færdiggøre koden nedenfor.

I det følgende ønsker vi at bruge egenvektor-dekompositionen til at beregne \(\boldsymbol{x}_k\) for \(k=0,\dots,K\) via

\[ \boldsymbol{x}_k = c_1 \lambda_1^k \boldsymbol{v}_1 + c_2 \lambda_2^k \boldsymbol{v}_2. \]

Først kræver det, at vi bestemmer konstanterne \(c_1\) og \(c_2\) for \(\boldsymbol{x}_0\).

Brug, at vektoren \(\begin{bmatrix}c_1\\c_2\end{bmatrix}\) løser ligningssystemet

\[\begin{split}\begin{bmatrix}\boldsymbol{v}_1&\boldsymbol{v}_2\end{bmatrix}\begin{bmatrix}c_1\\c_2\end{bmatrix}=\boldsymbol{x}_0\end{split}\]

Vink: Løsningen til et lineært ligningssystem \(A\boldsymbol{x}=\boldsymbol{b}\) kan bestemmes ved x=np.linalg.solve(A,b).
Vink: Søljevektorer \(\boldsymbol{x}\) og \(\boldsymbol{y}\) kan samles til en matrix \(M=\begin{bmatrix}\boldsymbol{x}&\boldsymbol{y}\end{bmatrix}\) ved M=np.column_stack([x,y]).

startpopulation = np.array([0., 1.])
# Saml systemmatrix
V = "INDSÆT KODE HER"

# Løs lingingssystemet
c = "INDSÆT KODE HER"

c1 = "INDSÆT KODE HER"
c2 = "INDSÆT KODE HER"

Følgende celle genererer et plot for antal unge og voksne kaniner for \(k=0,\ldots,K\).

Færdiggør cellen så \(\boldsymbol{x}_k\) beregnes via egenværdi-dekompositionen i hver iteration.

A = np.array([[0., 1.],
              [1., 1.]])

K = "INDSÆT KODE HER" # vælg fx 10 (beregn k=0..K)

# Beregn x_k vha formel
punkter = []
for k in range(K+1):
    xk = "INDSÆT KODE HER"
    punkter.append(xk)

punkter = np.array(punkter) # shape (K+1, 2)

# Plot komponenterne som funktion af k
plt.figure(figsize=(8,5))
plt.plot(range(K+1), punkter[:,0], 'o-', label='x[0] (unge)')
plt.plot(range(K+1), punkter[:,1], 'o-', label='x[1] (voksne)')
plt.xlabel('Måned (k)')
plt.ylabel('Antal kaniner')
plt.title('Udvikling af hver aldersklasse over tid')
plt.legend()
plt.grid(True)
plt.show()

Tjek at de værdier du får for \(\boldsymbol{x}_k\) er de samme, som dem du tidligere har bestemt. (Hvis du bruger egenvektorer og egenværdier som du har aflæst fra plottet, får du nok ikke nøjagtigt de samme tal.)

Hvis vi dividerer begge sider af egenvektor-dekompositionen med \(\lambda_1^k\), fås

\[ \frac{\boldsymbol{x}_k}{\lambda_1^k} = c_1 \boldsymbol{v}_1 + c_2 \left(\frac{\lambda_2}{\lambda_1}\right)^k \boldsymbol{v}_2. \]

Da \(|\lambda_2 / \lambda_1| < 1\), vil det andet led gå mod \(0\) for \(k\to\infty\), og dermed fås

\[ \frac{\boldsymbol{x}_k}{\lambda_1^k} \to c_1 \boldsymbol{v}_1 \quad \text{for} \; k\to \infty. \]

Det betyder, at vektoren \(\boldsymbol{x}_k\) vokser med omtrent \(\lambda_1^k\), og retningen af \(\boldsymbol{x}_k\) nærmer sig \(\boldsymbol{v}_1\).

Hvis vi definerer den demografiske profil som

\[\boldsymbol{d}_k = \frac{\boldsymbol{x}_k}{\|\boldsymbol{x}_k\|},\]

så vis med papir og blyant, at når \(\boldsymbol{v}_1\) er en enhedsvektor, så er

\[\boldsymbol{d}_k \to \boldsymbol{v}_1 \quad \text{for } k \to \infty.\]

Vis numerisk, at \(\boldsymbol{d}_k\) går mod \(\boldsymbol{v}_1\) for store \(k\) ved at færdiggøre koden nednefor.
Brug np.linalg.norm() til at finde længden af en vektor. Vink: Hvis grafen ikke går mod 0, så prøv at sammenligne med \(-\boldsymbol{v}_1\).

K = 'INDSÆT KODE HER'

# Beregn d_k = x_k / ||x_k|| for k = 0..K vha. egenvektor-formlen
dk_list = []
for k in range(K+1):
    xk = 'INDSÆT KODE HER'
    dk = 'INDSÆT KODE HER'
    dk_list.append(dk)
dk_arr = np.array(dk_list)

# Beregn afstand til reference og plottet (log-skala)
ks = np.arange(K+1)
afstand = np.linalg.norm(dk_arr - v1, axis=1)

plt.figure(figsize=(6,4))
plt.semilogy(ks, afstand, 'o-')
plt.xlabel('k')
plt.ylabel(r'$\|d_k - \hat v_1\|$')
plt.title('Afstand mellem $d_k$ og $\hat v_1$')
plt.grid(True)
plt.show()

2. Getting real - fødsel og overlevelse#

Indtil nu har vi arbejdet med Fibonacci-matricen, hvor alle elementer var \(1\) og populationen fulgte Fibonacci-følgen. Elementerne i matricen kan fortolkes som overlevelsesrater og fødselsrater. Vi generaliserer nu modellen ved at indføre.

\(f_{2}\): fødselsrate - forventet antal unger per voksen per md. (sandsynlighed for afkom),
\(s_{1}\): overlevelsesrate fra ung til voksen (sandsynlighed for at en ung overlever til at blive voksen næste md),
\(s_{2}\): overlevelsesrate for voksne (sandsynlighed for at en voksen overlever til næste md).

Den nye transitionsmatrix bliver:

\[\begin{split} A = \begin{bmatrix} 0 & f_{2} \\ s_{1} & s_{2} \end{bmatrix}. \end{split}\]

Vis ved at se på den karakteristiske ligning \(\det(A-\lambda I) = 0\) med papir og blyant, at hvis \(f_2, s_1, s_2\) alle er positive, da har matricen to reelle egenværdier, \(\lambda_1>0\) og \(\lambda_2<0\), og at \(\lambda_1 > |\lambda_2|.\)

I følgende opgave arbejder vi ud fra samme fremgangsmåde som i første opgave, hvor vi blot nu undersøger en anden transitionsmatrix. Du kan genbruge størstedelen af koden, som vi brugte i første opgave.

Opgaven går derfor ud på følgende:

Vælg nogle værdier for \(f_{2} \in (0.2, 0,7), s_{1}, s_{2} \in (0.7,1)\) og opbyg \(A\). Overvej, hvad gode biologiske værdier kunne være (brug evt. en chatbot til at hjælpe dig). Husk at definere \(A\) i en kodecelle.

Vælg en startpopulation. Husk at definere den i en kodecelle.

Brug enhedsvektor_afbildning(A) eller np.linalg.eig(A) til at finde egenvektorerne \(\boldsymbol{v}_1\) og \(\boldsymbol{v}_2\) og egenværdier \(\lambda_1\) og \(\lambda_2\) som tidligere.

Brug egenvektor-dekompositionen til at bestemme \(x_k\) for \(k=0,\dots,K\) (brug formlen \(\boldsymbol{x}_k = c_1\lambda_1^k \boldsymbol{v}_1 + c_2\lambda_2^k \boldsymbol{v}_2\)) og plot graferne for både \(\boldsymbol{x}_k[0]\) og \(\boldsymbol{x}_k[1]\).

Vis numerisk, at \(\boldsymbol{d}_k\) stadig går mod \(\boldsymbol{v}_1\) for store \(k\).

'INDSÆT KODE HER'

Når du har fået ovenstående til at fungere, kan du afprøve forskellige værdier - både for transitionsmatricen og for startpopulationen.

Flere aldersklasser#

Det viser sig, at en kaninpopulation med fordel kan inddeles i tre aldersklasser efter samme princip.
Inddelingen er:

\(\boldsymbol{x}_k[0]\): antal kaninunger i måned \(k\)
\(\boldsymbol{x}_k[1]\): antal unge voksne i måned \(k\)
\(\boldsymbol{x}_k[2]\): antal voksne i måned \(k\)

Det skyldes, at unge voksne har en relativ lav fødselsrate (ca. den halve af voksne) og overlevelsesraten i naturen fra ung voksen til voksen er relativt lav. Naturen er barsk, der er et stort tab.

Overgangene modelleres vha. transitionmatricen

\[\begin{split} A = \begin{bmatrix} 0 & f_2 & f_3\\[4pt] s_1 & 0 & 0\\[4pt] 0 & s_2 & s_3 \end{bmatrix}, \end{split}\]

hvor

\(f_2\) og \(f_3\) er fødselsrater: forventet antal unger per individ i aldersklasser hhv. 2 og 3 (månedligt).
\(s_1\) er overlevelsesraten fra unge til unge-voksne (klasse \(0 \rightarrow 1\)).
\(s_2\) er overlevelsesraten fra unge-voksne til voksne (klasse \(1 \rightarrow 2\)).
\(s_3\) er overlevelsesraten for voksne → voksne (klasse \(2 \rightarrow 2\)).

Vælg værdier for \(f_2,f_3,s_1,s_2,s_3\) og for startpopulationen og definér dem i kodecellen nedenfor. Find igen (fx ved hjælp af generativ AI) biologisk rimelige værdier.

f2 = 'INDSÆT KODE HER'
f3 = 'INDSÆT KODE HER'
s1 = 'INDSÆT KODE HER'
s2 = 'INDSÆT KODE HER'
s3 = 'INDSÆT KODE HER'
startpopulation = 'INDSÆT KODE HER'

Det viser sig, at sådan en matrix vi nu har defineret altid har en positiv egenværdi, som er dominerende. Og den tilhørende egenvektor kan vælges, sådan indgangene er ikke-negative. Det er en variant af Perron-Frobenius sætning.

Definér matricen \(A\) og brug lamb, V = np.linalg.eig(A) til at finde den største og dominerende egenværdi (lamb) og tilhørende egenvektorer (V).

A = 'INDSÆT KODE HER'

lamb, V = 'INDSÆT KODE HER'

# Sortérer så dominerende egenværdi kommer først
ind = np.argsort(np.abs(lamb))[::-1]
lamb = lamb[ind]
V = V[:, ind]

Vi kan som tidligere bruge egenværdi-dekompositionen på \(\boldsymbol{x}_0\), så

\[ \boldsymbol{x}_0 = c_1 \boldsymbol{v}_1 + c_2 \boldsymbol{v}_2 + c_3 \boldsymbol{v}_3. \]

Definér \(\boldsymbol{v}_1\), \(\boldsymbol{v}_2\) og \(\boldsymbol{v}_3\) i kodecellen nedenfor.

v1 = 'INDSÆT KODE HER'
v2 = 'INDSÆT KODE HER'
v3 = 'INDSÆT KODE HER'

Bestem koefficienterne \(c_1,c_2,c_3\) i egenværdi-dekompositionen af \(\boldsymbol{x}_0\) med np.linalg.solve(V, startpopulation).

# INDSÆT KODE HER

Brug egenvektor-dekompositionen til at bestemme \(\boldsymbol{x}_k\) for \(k=0,\dots,K\).

# INDSÆT KODE HER

Undersøg den demografiske profil \(\boldsymbol{d}_k\).

# INDSÆT KODE HER

Eksperimentér med valget af raterne og forsøg at svare på følgende:
Hvordan skal vi regulere raterne, sådan at vi har en stabil population?
Vink: Overvej tilfældende \(\lambda_1 <1,\) \(\lambda_1 =1\) eller \(\lambda_1 >1\).

3. Konkurrerende arter#

Nu introducerer vi en ny konkurrent - en anden kanin-art - der sættes ud i samme område. Arterne vil blande sig, parre sig på kryds og tværs og få afkom. Vi ønsker at undersøge, hvordan populationerne udvikler sig, som tiden går. Vi ser stadig på hunner, på to klasser (unger og voksne) og beskriver med \(\boldsymbol{x}_k = \begin{bmatrix}\boldsymbol{x}_k[0]\\ \boldsymbol{x}_k[1]\end{bmatrix}\) den oprindelige art som tidligere og med \(\boldsymbol{y}_k = \begin{bmatrix}\boldsymbol{y}_k[0]\\ \boldsymbol{y}_k[1]\end{bmatrix}\) den nye art.

De to arter blander sig og får afkom på tværs. Vi modellerer efter det princip, at afkom, hvor mindst en forælder er af den nye art, tilhører den nye art.

Dette giver en udvidet dynamik:

\[\begin{split} \boldsymbol{z}_{k+1} = \begin{bmatrix} \boldsymbol{x}_{k+1}[0] \\ \boldsymbol{x}_{k+1}[1] \\ \boldsymbol{y}_{k+1}[0] \\ \boldsymbol{y}_{k+1}[1] \end{bmatrix} = \begin{bmatrix} 0 & a_{12} & 0 & 0 \\ a_{21} & a_{22} & 0 & 0 \\ 0 & a_{32} & 0 & a_{34} \\ 0 & 0 & a_{43} & a_{44} \end{bmatrix} \begin{bmatrix} \boldsymbol{x}_k[0] \\ \boldsymbol{x}_k[1] \\ \boldsymbol{y}_k[0] \\ \boldsymbol{y}_k[1] \end{bmatrix} = A \boldsymbol{z}_k \end{split}\]

Selektion og if-elif-else udsagn i Python#

Inden vi kigger på populationsudviklingen for de to arter, vil vi fokusere på en kontrolstruktur, som hedder selektion.

Det er en kontrolstruktur, som er helt central når man programmerer. Basalt set handler det om at sætte en betingelse i programmet, sådan, at hvis betingelsen er opfyldt, så foretager vi en handling, og hvis betingelsen ikke er opfylkdt, så foretager vi en anden handling (eller ingen handling).

Grundformen er:

if første_betingelse:
    # kør dette hvis første betingelse er sand
elif anden_betingelse:
    # ellers hvis anden betingelse er sand
else:
    # ellers kør dette

Man kan udelade både elif og else efter behov.

Fuldfør kodecellen nedenfor, som ud fra et heltal x gør følgende:

Hvis x er større end 0, udskriv "Tallet er positivt".

Hvis x er mindre end 0, udskriv "Tallet er negativt".

Ellers (hvis x er lig 0), udskriv "Tallet er nul"

x  = np.random.standard_normal() # Vi vælger et tilfældigt tal efter en standard normalfordeling

if x>0:
    print("Tallet er positivt , x = ",x)
elif "INDSÆT KODE HER":
    print("Tallet er negativt, x = ",x)
else:
    print("Tallet er nul")

Test at kodecellen virker for både positive x, negative x og hvis x er 0.

Et biologisk princip hedder “Competitive Exclusion Principle”. Det handler om, at to arter, som konkurrerer og formerer sig baseret på samme ressourcer ikke kan sameksistere. En af arterne vil overtage, og den anden vil uddø. Det er en af årsagerne til, at arter udvikler nicheadskillelse ved hjælp af evolution.

Vi vil se dette princip udfolde sig for vores to kaninpopulationer.

Vi er nu klar til at undersøge udviklingen af de to konkurrerende arters populationer. I kodecellen nedenfor er der valgt værdier for \(a_{12},a_{21},a_{22},a_{32},a_{34},a_{43},a_{44}\) og for en startpopulation for \(z_0\). De to populationer udvikler sig som tidligere, men det nye er, at en del af afkommet fra den oprindelige art tilhører den nye art.

# Rater
a12 = 0.4    # x-voksne -> x-unger
a21 = 0.55   # x-unger -> x-voksne
a22 = 0.6    # x-voksne -> x-voksne

a32 = 0.2    # x-voksne -> y-unger
a34 = 0.4    # y-voksne -> y-unger
a43 = 0.6    # y-unger -> y-voksne
a44 = 0.6   # y-voksne -> y-voksne

# Startpopulation
startpopulation = np.array([2000, 3000, 1, 2])

I de tidligere opgaver har vi udnyttet egenværdi-dekompositionen til iteration. I resten af opgaverne vil vi igen se på den simple matrix-iteration \(\boldsymbol{z}_{k+1}=A\boldsymbol{z}_k\)

Implementér transitionsmatricen \(A\) i kodecellen nedenfor.

A = 'INDSÆT KODE HER'

Færdiggør kodecellen nedenfor, så der gøres føgelende i hver iteration for \(k=0, \dots, K\):

beregner \(\boldsymbol{z}_{k+1}\) .

beregner total population for \(\boldsymbol{x}\), \(\boldsymbol{y}\) og den samlede total population (både \(\boldsymbol{x}\) og \(\boldsymbol{y}\)).

benytter et if-statement til at printe de år, hvor art \(\boldsymbol{y}\) udgør mere end halvdelen af populationen.

Bemærk: Vi definerer i koden en matrix \(Z\), så den \(k\). række er \((Z)_k=\boldsymbol{z}_k\). I Python får vi så \(\boldsymbol{z}_k\) som Z[k] og elementerne \(\boldsymbol{z}_k[i]\) som Z[k,i] for \(i\in\{0,1,2,3\}\).

# Iterér over K år
K = 'INDSÆT KODE HER'
Z = np.zeros((K+1, 4))
Z[0] = startpopulation

for k in range(K):
    Z[k+1] = 'INDSÆT KODE HER'
    
    # Beregn total
    total_x = Z['INDSÆT KODE HER', 'INDSÆT KODE HER'] + Z['INDSÆT KODE HER', 'INDSÆT KODE HER'] # x_{k+1}[0] + x_{k+1}[1]
    total_y = Z['INDSÆT KODE HER', 'INDSÆT KODE HER'] + Z['INDSÆT KODE HER', 'INDSÆT KODE HER'] # y_{k+1}[0] + y_{k+1}[1]
    total = 'INDSÆT KODE HER' # # x_{k+1}[0] + x_{k+1}[1] + y_{k+1}[0] + y_{k+1}[1]
    
    # Tjek om art y udgør mere end halvdelen
    if 'INDSÆT KODE HER':
        print(f"Måned {k+1}: Art Y udgør >50% (andel = {total_y/total:.2f})")
        #break
    else:
        print(f"Måned {k+1}: Art Y <50% (andel = {total_y/total:.2f})")

I if-statementet er der følgende linje:

# break

Prøv at fjern # og kør kodecellen igen.

Forklar ud fra dette, hvad et break-statement gør og brug dette til at finde det mindste år, hvor \(y\) udgør over halvdelen af den samlede population.

Du kan bruge nedenstående funktion til at se udvikling af de to arters population. Hvis du har fjernet # i ovenstående kodecelle fra break-statementet, skal du skrive det igen.

def plot_totalpopulation(Z):
    år = np.arange(Z.shape[0])
    total_x = Z[:,0] + Z[:,1]
    total_y = Z[:,2] + Z[:,3]

    fig, ax = plt.subplots(1, 2, figsize=(12,5), sharey=True)

    # Art X
    ax[0].semilogy(år, total_x, 'o-', color='tab:blue')
    ax[0].set_xlabel('Måneder')
    ax[0].set_ylabel('Antal individer')
    ax[0].set_title('Art X (totalpopulation)')
    ax[0].grid(True)

    # Art Y
    ax[1].plot(år, total_y, 'o-', color='tab:orange')
    ax[1].set_xlabel('Måneder')
    ax[1].set_title('Art Y (totalpopulation)')
    ax[1].grid(True)

    plt.tight_layout()
    plt.show()


plot_totalpopulation(Z)

Hvor mange måneder går der cirka før art \(\boldsymbol{x}\) udgør mindre end 10% af den samlede population?

I den følgende kodecelle definerer vi en funktion dominating_eigenvalue, der returnerer den dominerende egenværdi og den tilhørende egenvektor for en given matrix \(A\).

Prøv at forudsige den demografiske profil ved hjælpe af den dominerende egenværdi og -vektor.
Vil den oprindelige art altid uddø?

def dominating_eigenvalue(A):
    """
    Beregner den dominerende egenværdi og tilhørende egenvektor for matrix A.
    """
    eigenvalues, eigenvectors = np.linalg.eig(A)

    # Finder indeks for den største egenværdi (absolut værdi)
    idx = np.argmax(np.abs(eigenvalues))

    # største egenværdi-par
    largest_eigenvalue = eigenvalues[idx]
    largest_eigenvector = eigenvectors[:, idx]

    # Normaliserer til norm 1 og positiv sum
    largest_eigenvector = largest_eigenvector / np.sum(largest_eigenvector)
    largest_eigenvector = largest_eigenvector / np.linalg.norm(largest_eigenvector)

    # tager realdel
    largest_eigenvalue = np.real(largest_eigenvalue)
    largest_eigenvector = np.real(largest_eigenvector)

    return largest_eigenvalue, largest_eigenvector

# Print results
largest_eigenvalue, largest_eigenvector = dominating_eigenvalue(A)
print("Largest eigenvalue:", largest_eigenvalue)
print("Corresponding eigenvector:", largest_eigenvector)

4. Den danske befolkning#

Vi udvider nu vores aldersstrukturerede modeller til en større population. Vi vælger at betragte hele den danske befolkning. Populationen beskrives ved en vektor \(\boldsymbol{x} \in \mathbb{R}^{100}\), hvor hvert element \(\boldsymbol{x}[i]\) svarer til antallet af personer i alderen \(i\). Vi betrager altså kun mennesker i alderen 0-99 år, men lader \(\boldsymbol{x}[99]\) beskrive aldrene \(\geq\) 99 år. Dynamikken

\[ \boldsymbol{x}_{k+1} = A \boldsymbol{x}_k + \boldsymbol{b} \]

styres af en Leslie-matrix \(A \in \mathbb{R}^{100 \times 100}\).

Transitionsmatricen \(A\) (der som nævnt er en Leslie-matrix) har følgende form:

Første række af \(A\) indeholder fødselsrater \(f_{15}, \dots, f_{51} > 0\), altså aldersklasser hvor folk typisk får børn.
Underdiagonal indeholder overlevelsesrater \(s_i > 0\), som angiver sandsynligheden for at personer i aldersklasse \(i\) overlever til aldersklasse \(i+1\).
Alle andre elementer i \(A\) er nul.

Dvs. \(A\) har formen

\[\begin{split} A = \begin{bmatrix} f_1 & f_2 & f_3 & \cdots & f_n \\ s_1 & 0 & 0 & \cdots & 0 \\ 0 & s_2 & 0 & \cdots & 0 \\ \vdots & & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & s_{n-1} & 0 \end{bmatrix}. \end{split}\]

\(\boldsymbol{b}\) kan bruges til at tilføje migration eller andre eksterne tilførsler.

Befolkningsdata#

I csv-filen ‘populationsdata.csv’ har vi samlet befolkningsdata fra den danske befolkning i 2024. Tallene er alle fra Danmarks Statistik. I følgende kodecelle importeres filen, og de øverste rækker i datasættet printes.

Kør cellen, og se hvad datasættet indeholder.

df = pd.read_csv('Data/populationsdata.csv')
print("Kolonner i datafilen:\n",df.head())

I første omgang er vi interesserede i kollonnen Befolkning, der er data for aldersfordelingen af den danske befolkning. Det er denne data vi vil bruge som vores startpopulation \(\boldsymbol{x}_0\).

Vi kan definere en kollonne som en \(\mathbb{R}^{100}\) vektor i et Numpy-array.

# Konverter til NumPy arrays
alder = df['Alder'].to_numpy()
befolkning = df['Befolkning'].to_numpy()

print("Størrelsen af datasættet:\n", len(alder))

Følgende celle generer et plot over aldersfordelingen i Danmark.

# Plot startpopulation
plt.figure(figsize=(8,4))
plt.bar(alder, befolkning)
plt.xlabel("Alder")
plt.ylabel("Antal personer")
plt.title("Aldersfordeling i startpopulationen (0-99 år)")
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.show()

Valg af eget datasæt (Valgfri)#

Denne sektion er en valgfri øvelse, hvis man ønsker at se, hvordan man kan importere et datasæt til Python. Her er det nødvendigt at du arbejder på JupyterLite-siden, hvor du kan uploade en CSV-fil ved at “drag-and-droppe”.

Hvordan vil modelleringen af populationen se ud, hvis du tager udgangspunkt i din hjemkommunes aldersfordeling som startpopulation \(\boldsymbol{x}_0\)?

Opgaven er som følger:

Åbn følgende link til Danmarks Statistik.

Find din hjemkommune under kategorien OMRÅDE.

Vælg Markér alle under kategorien ALDER.

Klik på VIS TABEL.

Skift fra Excel (*.xlsx) til Matrix (*.csv) og klik på den blå pil lige ved siden af for at downloade CSV-filen.

Læg CSV-filen i samme directory som denne Jupyter-notebook.

Ændre INDSÆT KODE HER til filnavnet på CSV-filen i nedenstående kodecelle.

# Indlæs CSV
df = pd.read_csv("INDSÆT KODE HER", sep=";", quotechar='"', encoding="latin1", header=None)

# Hent startpopulation for 0-99 år
# Kolonne 0: alder, kolonne 1: antal personer
# Skipper første række "Alder i alt"
startpopulation = df.iloc[1:101, 1].astype(int).to_numpy()  # 0-99 år

# Plot startpopulation
plt.figure(figsize=(10,5))
plt.bar(alder, startpopulation)
plt.xlabel("Alder")
plt.ylabel("Antal personer")
plt.title("Aldersfordeling i startpopulationen (0-99 år)")
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.show()

Leslie-matricen: En simpel model#

Vi ønsker at undersøge, hvordan befolkningen udvikler sig over tid ved hjælp af Leslie-matricen. Vi starter med en forenklet model, hvor fødselsrater og overlevelsesrater antages at være uniforme.

Definer Leslie-matricen. Vælg selv værdier for:

Uniform fødselsrate for aldersklasser 15-51 år.

Uniform overlevelsesrate for alle aldersklasser.

# Definér Leslie-matrix
k = 100
A = np.zeros((k, k))

# Fødselsrater
for i in range(15, 52):
    A[0, i] = 0.05 #'INDSÆT KODE HER'

# Overlevelsesrater
for i in range(k-1):
    A[i+1, i] = 0.95 #'INDSÆT KODE HER'

I følgende kodecelle defineres en funktion population_over_tid, der simulerer en populations udvikling over tid givet en transitionmatrix, og en startpopulation.

Funktionen generer også to figurer:

total befolkning over tid.
aldersfordelingen efter \(K\) år.

Færdiggør funktionen.

def population_over_tid(A, startpopulation, K, b = None, plot=True):
    """
    Beregner population over K år givet transitionmatrix A, startpopulation og migrationsvektor b.
    Returnerer en matrix med population for hver aldersklasse over tid.
    """
    k = A.shape[0]
    X = np.zeros((K+1, k))
    X[0] = startpopulation

    # sætter b til nulvektor hvis ikke givet
    if b is None:
        b = np.zeros(k)

    # iterationstep
    for år in range(K):
        X[år+1] = 'INDSÆT KODE HER'

    # generér plots hvis ønsket
    if plot:
        # plot total population over tid
        total_population = X.sum(axis=1)
        plt.figure(figsize=(8,4))
        plt.plot(range(K+1), total_population, 'o-')
        plt.xlabel('År')
        plt.ylabel('Total befolkning')
        plt.title('Udvikling af den danske befolkning')
        plt.grid(True)
        plt.show()

        # plot aldersfordeling efter K år
        plt.figure(figsize=(8,4))
        plt.bar(range(k), X[K])
        plt.xlabel('Alder')
        plt.ylabel('Antal personer')
        plt.title(f'Aldersfordeling efter {K} år')
        plt.show()

    return X

Kør funktionen for din matrix \(A\), og se plottene. Vælg selv antal år \(K\).

_ = population_over_tid(A, befolkning, K=100)

Hvordan er den langvarige udvikling af totalbefolkningen? Stigende/faldende/stabil?

For at kunne diskutere den langvarige udvikling af befolkningen introducerer vi begrebet nettoreproduktionsrate.

Nettoreproduktionsrate#

Nettoreproduktionsraten (NRR) er det gennemsnitlige antal afkom, som en nyfødt individ forventes at producere i løbet af sin levetid. I praksis modelleres der som regel kun hunner, da de bestemmer den effektive reproduktion i populationen. Genkald formen af en Leslie-matrix

\[\begin{split} A = \begin{bmatrix} f_1 & f_2 & f_3 & \cdots & f_n \\ s_1 & 0 & 0 & \cdots & 0 \\ 0 & s_2 & 0 & \cdots & 0 \\ \vdots & & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & s_{n-1} & 0 \end{bmatrix}. \end{split}\]

For en \(n\times n\) Leslie-matrix beregnes nettoreproduktionsraten \(R_0\) som

\[ R_0 = \sum_{i=1}^{n} f_i \, l_i , \]

hvor \(f_i\) er fertilitetsraten for aldersklasse \(i\), og \(l_i\) er sandsynligheden for at overleve til begyndelsen af aldersklasse \(i\). Overlevelsessandsynligheden \(l_i\) defineres rekursivt som

\[ l_1 = 1, \qquad l_2 = s_1, \qquad l_3 = s_1 s_2, \qquad \dots, \qquad l_i = \prod_{j=1}^{i-1} s_j , \]

hvor \(s_j\) er sandsynligheden for at overleve fra aldersklasse \(j\) til \(j+1\). Dermed kan \(R_0\) også skrives eksplicit som

\[ R_0 = f_1 + f_2 s_1 + f_3 s_1 s_2 + \dots + f_n s_1 s_2 \cdots s_{n-1} =\sum_{i=1}^n f_i\left( \prod_{j=1}^{i-1}s_j \right). \]

I følgende kodecelle definerer vi en funktion, der beregner nettoreproduktionsraten for en given Leslie-matrix.

Færdiggør funktionen.

def NRR(A):
    """
    Beregner nettoreproduktionsraten R0 for Leslie-matrix A.
    """
    n = A.shape[0]
    R0 = 0.0

    # Beregn R0
    for i in range(n):
        # fødselsrate for aldersklasse i
        fi = "INDSÆT KODE HER"

        # overlevelsesrate fra fødsel til aldersklasse i
        li = 1.0
        for j in range(i):
            # overlevelsesrate for aldersklasse j til j+1
            sj = "INDSÆT KODE HER"
            li *= sj
        R0 += fi * li

    return R0

Fortolkningen af \(R_0\) er, at

hvis \(R_0>1\) forventes populationen af vokse.
hvis \(R_0=1\) er populationen stabil.
hvis \(R_0<1\) forventes populationen at falde.

Brug funktionen til at undersøge nettoreproduktionsraten. Svarer fortolkningen til hvad du kunne se af plottene?

"INDSÆT KODE HER"

Vi har de samme fortolkninger for den dominerende egenværdi, og om den er over/under 1.

Undersøg den dominerende egenværdi vha. dominating_eigenvalue som defineret tidligere.

Svarer fortolkningerne til hindanden? Er \(R_0=\lambda_{max}\)?

"INDSÆT KODE HER"

Prøv at gå tilbage til hvor Leslie-matricen blev defineret. Kan du finde værdier for fødsel- eller overlevelsesraterne, hvor nettoreproduktionsraten er tæt på 1? Vis plottene igen, og se den langvarige udvikling.

Udvidelse af Leslie-modellen#

For bedre at afspejle virkeligheden vil vi nu eksperimentere med en ikke-uniform fødsels- og overlevelsesrate, som varierer med alderen.

Første udvidelse er, at i stedet for at vælge en konstant rate for alle fertile aldersklasser, lader vi raten følge en parabel, der har top i alder 33 og nul i alderen 14 og 52.

Brug denne parabel til modellering af fødselsraten til at opdatere Leslie-matrixen i koden nedenfor. Lad f_max være værdien af fødselsraten i toppunktet. Brug fortsat uniform overlevelsesrate.

# Definér Leslie-matrix
k = 100
A_new = np.zeros((k, k))

# toppunkt af fødselsrater
f_max = "INDSÆT KODE HER"

# Fødselsrater
for i in range(15, 52):
    A_new[0, i] = 'INDSÆT KODE HER'

# Overlevelsesrater
for i in range(k-1):
    A_new[i+1, i] = 'INDSÆT KODE HER'

Undersøg udviklingen af befolkningen vha. af population_over_tid og nettoreproduktionsraten.

R0 = NRR(A_new)
print("Nettoreproduktionsrate R0 for udvidede model:", R0)

_ = population_over_tid(A_new, befolkning, K=100)

Giver den nye fødselsrate en forskel fra figurerne fra den simple model?

Vi vil nu udvide vores model yderligere: Lad overlevelsesraten være en lineær funktion med højeste rate ved 0 år, faldende mod ældre aldersklasser.

Indfør den nye overlevelsesrate ovenfor, hvor A_new defineres. Vælg selv eksperimentere med overlevelsesraten i år 0 og år 99.

Prøv at finde værdier for f_max og overlevelsesraten for aldrene 0 år og 99 år så Leslie-matricen har nettoreproduktionsrate \(R_0\approx 1\).
Sammenlign igennem plots med den simple model med \(R_0\approx 1\).

Databaseret populationsmodellering#

Vi ønkser nu at anvende vores datasæt til at skabe en Leslie-matrix, der afspejler den danske population bedst muligt. Vi laver i følgende celle Numpy-arrays for de restende kolonner i datasættet.

dødstal = df['Dødstal'].to_numpy()
fertilitet = df['Fertilitet'].to_numpy()
indvandring = df['Indvandring'].to_numpy()
udvandring = df['Udvandring'].to_numpy()

Kollonnerne beskriver følgende for hver aldersklasse:

Dødstal: Antal dødsfald.
Fertilitet: Fødte børn per 1000 kvinder.
Indvandring: Antal indvandrede personer.
Udvandring: Antal udvandrede personer.

I Leslie-matricen ønsker vi fødselsraten per voksen. Overlevelsesraten til Leslie-matricen kan findes ud fra dødstallene og befolkningstallene.

Færdiggør koden i følgende celle så vi får fødsels- og overlevelsesrater på den ønskede form.
Vink: Antag at halvdelen af befolkningen er kvinder.

# fertilitet er pr 1000 kvinder, skal omregnes til pr person
fødselsrate = "INDSÆT KODE HER"

# overlevelsesrate beregnes fra befolkning og dødstal
overlevelsesrate = "INDSÆT KODE HER"

I følgende celler vises fødselraten, overlevelsesraten og befolkningstallet.

Kør cellen.

# Fødselsrate plot
plt.figure(figsize=(6,3))
plt.plot(alder, fødselsrate, 'o-')
plt.title('Fødselsrate')
plt.ylabel('Rate')
plt.xlabel('Alder')
plt.show()

# Overlevelsesrate plot
plt.figure(figsize=(6,3))
plt.plot(alder, overlevelsesrate, 'o-')
plt.title('Overlevelsesrate')
plt.xlabel('Alder')
plt.ylabel('Rate')
plt.show()

Ser det ud, som man kunne forvente?

Hvor godt stemmer det overens med de rater vi brugte i de forrige Leslie-matricer?

Vi kan nu samle Leslie-matricen \(A\).

k = 100
A = np.zeros((k, k))
for i in range(k):
    A[0, i] = fødselsrate[i]
for i in range(k-1):
    A[i+1, i] = overlevelsesrate[i]
A[-1,-1] = overlevelsesrate[-1]  # Sidste aldersklasse er selvabsorberende

Undersøg igen udviklingen af befolkingen vha. population_over_tid og NRR.

"INDSÆT KODE HER"

Er populationen stabil eller voksende/faldende?

Hvordan har aldersfordelingen ændret sig i forhold til startpopulationen?

Forestil dig nu, at vi kan lave en kampagne der skal få danskerne til at få flere/færre børn. Kan vi opnå en nettoreproduktionsrate tæt på 1?

I følgende celle definerer vi en ny Leslie-matrix A_new, hvor fødselsraterne ganges med en faktor f_faktor. Kan du finde en værdi så \(R_0\approx 1\)?

f_faktor = 'INDSÆT KODE HER' 

A_new = A.copy()
A_new[0] = f_faktor * A_new[0] 

R0 = NRR(A_new)
print("Nettoreproduktionsrate R0 med ændret fertilitet:", R0)

Vis plottene igen for den nye Leslie-matrix, og se den langvarige udvikling.

_ = population_over_tid(A_new, befolkning, K=1000)

Migration#

I stedet for at ændre fødsels- og overlevelsesrater i Leslie-matricen kan befolkningens udvikling og påvirkes af eksterne påvirkninger, dvs. migration. Indtil videre har vi ikke taget højde for migrationsleddet \(\boldsymbol{b}\) i

\[ \boldsymbol{x}_{k+1}=A\boldsymbol{x}_k+\boldsymbol{b}. \]

Definer i følgende celle vektoren \(\boldsymbol{b}\) ud fra kollonnerne indvandring og udvandring og vis migrationsdataen i det givne plot.

b = indvandring-udvandring#"INDSÆT KODE HER"

plt.figure()
plt.plot(indvandring, label='Indvandring')
plt.plot(udvandring, label='Udvandring')
plt.plot(b, label='b')
plt.title('Migration')
plt.xlabel('Alder')
plt.ylabel('Antal personer')
plt.legend()
plt.show()

Hvis der findes en ligevægt \(\boldsymbol{x}_{lim}\) så \(\boldsymbol{x}_k\to\boldsymbol{x}_{lim}\) for \(k\to\infty\), så gælder at

\[ \boldsymbol{x}_{lim}=A\boldsymbol{x}_{lim}+\boldsymbol{b}\implies(I-A)\boldsymbol{x}_{lim}=\boldsymbol{b}. \]

Altså kan vi finde \(\boldsymbol{x}_{lim}\) som løsningen af et lineært ligningssystem.

Løs ligningssystemet for \(\boldsymbol{x}_{lim}\).
Vink: \(A\boldsymbol{x}=\boldsymbol{b}\) kan løses med np.linalg.solve(A,b).
Vink: \(n\times n\) identitesmatricen skabes med np.eye(n).

xlim = "INDSÆT KODE HER" 

Vis ligevægts befolkningen i det følgende plot.

plt.figure(figsize=(8,4))
plt.bar(alder, xlim)
plt.xlabel("Alder")
plt.ylabel("Antal personer")
plt.title("Aldersfordeling i startpopulationen (0-99 år)")
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.show()

Kan du verificere dette i simuleringer? Vis befolkningsudviklingen med population_over_tid(A, startpopulation, K, b) og giv migrationsvektoren som input b.

"INDSÆT KODE HER"

Opnår vi ligevægt?

Sammenlign med det tilsvarende plot uden migrationsled.

Beregningskompleksitet af matrix-iteration#

Vi har tidligere i opgaven undersøgt populationen ved hjælp af Leslie-matricen og beregnet populationen over tid ved iterativt at anvende

\[ \boldsymbol{x}_{k+1} = A \boldsymbol{x}_k = A (A (... (A \boldsymbol{x}_0) ...)). \]

Vi har altså ikke beregnet \(A^{k+1} \boldsymbol{x}_0\) direkte, men gentaget matrix-vektor-multiplikationen \(k+1\) gange.

I denne opgave skal vi estimere, hvor mange grundlæggende beregninger (multiplikationer og additioner) der kræves for en \(100 \times 100\) Leslie-matrix.

Brug \(k = 50\)
Dimension: \(A\) er \(100 \times 100\)

Vi definerer først \(k\) og dimensionen \(n\).

k = 50
n = 100

Overvej for matrix-vektor produktet:

Hvor mange multiplikationer bruges?

Hvor mange additioner bruges?

Vink: For dot-produktet mellem to vektorer i \(\mathbb{R}^n\) bruges \(n\) multiplikationer og \(n-1\) additioner.

Definer resultatet i følgende celle.
Vink: I Python skrives \(x^a\) som x**a.

mv_mul = "INDSÆT KODE HER"
mv_add = "INDSÆT KODE HER"
mv_ops = mv_mul + mv_add

Overvej det samme for et matrix-matrix produkt og definer resultatet.

mm_mul = "INDSÆT KODE HER"
mm_add = "INDSÆT KODE HER"
mm_ops = mm_mul + mm_add

Hvor mange hhv. matrix-matrix og matrix-vektor produkter anvendes til at regne \((A^{k+1})\boldsymbol{x}_0\) direkte?

Brug dette til at finde det totale antal operationer i udregningen og definer resultatet.

direkte_total = "INDSÆT KODE HER"

Overvej på samme måde antallet af hhv. matrix-matrix og matrix-vektor produkter, der bruges til den iterative metode \(\boldsymbol{x}_{k+1}=A(A(...A(A\boldsymbol{x}_0)))\).

Definer det totale antal operationer.

iterativ_total = "INDSÆT KODE HER"

Følgende celle printer de to resultater og forholdet mellem dem.