Modul 3 Elever

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

Modul 3 Elever#

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 rotte. 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. 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#

  1. Opvarmning: Fibonacci-tallene genbesøgt

  2. Getting real - fødsel og overlevelse

  3. 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 antal af voksne hunkaniner \(F_k = \boldsymbol{x}_k[1]\), genfinder vi rekursionen

\[ F_{k+1} = F_k + F_{k-1}, \quad k = 1,2,\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 voksen-bestanddel 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.
Du kan bruge @ til at lave matrix-vektor produkter. Du kan bruge np.array() til at definere \(A\) og \(x_0\).

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|.\) Indsæt værdierne i kodecellen nedenfor:

lambda1 = "INDSÆT KODE HER"
lambda2 = "INDSÆT KODE HER"

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)

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

To vektorer \(\boldsymbol{v}_1\) og \(\boldsymbol{v}_2\) er lineært uafhængige, hvis

\[ \alpha \boldsymbol{v}_1 + \beta \boldsymbol{v}_2 = \boldsymbol{0} \ \Rightarrow\ \alpha = 0 \ \text{og} \ \beta = 0. \]

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.

Asymptotisk udvikling#

Vi kan skrive startvektoren som en linearkombination af egenvektorerne \(\boldsymbol{v_1}\) og \(\boldsymbol{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. \]

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{split} \begin{bmatrix}c_1\\c_2\end{bmatrix} \end{split}\]

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ærdidekompositionen i hver iteration.
Vink: Man kan beregne \(a^b\) i Python ved a**b.

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 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. \]

Hint

Brug udtrykket

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

og faktoriser \(\lambda_1^k\) ud af både tæller og nævner i

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

Udnyt til sidst, at

\[ \left(\frac{\lambda_2}{\lambda_1}\right)^k \to 0 \quad \text{når } k\to\infty. \]

Vis numerisk, at \(\boldsymbol{d}_k\) går mod \(\boldsymbol{v}_1\) for store \(k\) ved at færdiggøre koden nedenfor.
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(r'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 ikke-nul 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|.\)

Hint

Den karakteristiske ligning kommer fra

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

For matricen

\[\begin{split} A = \begin{bmatrix} 0 & f_2 \\ s_1 & s_2 \end{bmatrix} \end{split}\]

får du

\[\begin{split} \det \begin{bmatrix} -\lambda & f_2 \\ s_1 & s_2 - \lambda \end{bmatrix} = (-\lambda)(s_2-\lambda) - f_2 s_1 = 0. \end{split}\]

Omskriv, så du får en andengradsligning i \(\lambda\) - og brug fortegnene på \(f_2, s_1, s_2\) til at argumentere for fortegnene på rødderne samt hvilken der har størst absolutværdi.

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:

  1. 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.

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

  3. 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.

  4. 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]\).

  5. 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.

3. 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.

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

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

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.

Vi ønsker at bruge den nuværende aldersfordeling af den danske befolkning som \(\boldsymbol{x}_0\). Derfor skal du på arbejde JupyterLite-siden, hvor du kan uploade en CSV-fil ved at “drag-and-droppe”.

Opgaven er som følger:

  1. Åbn følgende link til Danmarks Statistik.

  2. Vælg Markér alle under kategorien ALDER.

  3. Klik på VIS TABEL.

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

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

  6. Æ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

Kør kodecellen nedenfor for at visualisere startpopulationen:

# Antal aldersklasser
k = 100

# x-akse: alder 0-99
aldre = np.arange(k)

# Plot startpopulation
plt.figure(figsize=(10,5))
plt.bar(aldre, 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()

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] = 'INDSÆT KODE HER'

# Overlevelsesrater
for i in range(k-1):
    A[i+1, i] = '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.

Kør cellen.

def population_over_tid(A, startpopulation, K, b = None, plot=True):
    """
    Beregner population over K år givet transitionmatrix A og startpopulation.
    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] = A @ X[år] + b

    # generér plots hvis ønsket
    if plot:
        # plot total population over tid
        total_population = X.sum(axis=1)
        plt.figure(figsize=(8,5))
        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=(10,5))
        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, startpopulation, K='INDSÆT KODE HER')

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.

Hint

Tænk på, hvordan værdierne i en Leslie-matrix er organiseret:

  • Den første række indeholder fertilitetsrater for hver aldersklasse.

  • De andre rækker indeholder overlevelsesrater, som viser sandsynligheden for at gå fra én aldersklasse til den næste.

Overvej, hvordan du kan hente disse værdier i en Python-matrice med A[i,j] og bruge dem i løkkerne for fi og sj.

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.

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\).

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)

Undersøg den dominerende egenværdi vha. dominating_eigenvalue.

Svarer fortolkningerne til hinanden? 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.

Ny fødselsrate som parabel#

Vi vil nu eksperimentere med en ikke-uniform fødselsrate, som varierer med alderen. 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 32 og nul i alderen 14 og 52.

Brug denne parabel til modellering af fødselsraten til at opdatere Leslie-matrixen i koden nedenfor, hvor du igen får

  • Et plot for total befolkning over tid.

  • Et plot, der viser, hvordan populationen fordeler sig på aldersklasser efter \(K\) år.

Du kan selv eksperimentere med, hvad toppunktsværdien skal være.

Hint

Tænk på fødselsraten som en parabel, der er nul ved de yngste og ældste fertile aldersklasser og har et maksimum ved den alder, hvor fertiliteten er størst.

\[ f(i) = f_{\max} \cdot \frac{(i - i_\text{min})(i_\text{max} - i)}{(i_\text{top} - i_\text{min})(i_\text{max} - i_\text{top})} \]

hvor \(i_\text{min}\) og \(i_\text{max}\) er de yderste fertile aldre, \(i_\text{top}\) er toppen, og \(f_{\max}\) er den maksimale rate.

# Antal aldersklasser
k = 100

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

# Fødselsrater: udfyld parablen her
for i in range(15, 52):
    # beregn A[0, i] som parabel med top i 32 og nul i 14 og 52
    A[0, i] = 'INDSÆT KODE HER'

# Overlevelsesrater
for i in range(k-1):
    A[i+1, i] = 0.95

# Beregn population over K år
X = population_over_tid(A, startpopulation, K=100)

Ny dødsrate som lineær funktion#

Vi vil nu udvide vores model af den danske befolkning yderligere.

  • Fødselsraten skal være en parabel, som vi lavede tidligere (top ved 32 år, nulpunkt ved 14 og 50 år).

  • Nyt tillæg: Dødsraten skal være en lineær funktion med laveste dødsrate ved 0 år, stigende mod ældre aldersklasser.

Indfør den nye dødsrate i koden nedenfor. Du vil igen få:

  • Et plot, der viser den samlede befolkning over tid.

  • Et plot, der viser aldersfordelingen af populationen efter \(K\) år.

Du kan selv eksperimentere med dødsraten i år 0 og år 99.

Hint

En lineær funktion kan skrives generelt som

\[ y = y_\text{max} - (y_\text{max}-y_\text{min}) \frac{x}{x_\text{max}}. \]

Formlen sikrer, at overlevelsesraten falder lineært fra høj værdi ved unge til lav værdi ved gamle.
Du kan bruge denne lineære formel inde i løkken for A[i+1, i].

# Antal aldersklasser
k = 100

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

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

# Lineær overlevelsesrate
for i in range(k-1):
    A[i+1, i] = 'INDSÆT KODE HER'

# Beregn population over K år
X = population_over_tid(A, startpopulation, K=100)

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.

Kør cellen.

print("Direkte metode:\n",direkte_total)
print("Iterativ metode:\n",iterativ_total)
print("Sammenligning:\nDirekte/Iterativ operationer:",direkte_total / iterativ_total)
Contents