--- title: "Εκθετική Αύξηση Πληθυσμού" subtitle: "Από τα νούφαρα της λίμνης στη διαφορική εξίσωση" author: "Εργαστήριο Πληθυσμιακής Οικολογίας — Τμήμα Βιολογίας" date: "`r format(Sys.Date(), '%d %B %Y')`" output: html_document: toc: true toc_float: true toc_depth: 3 number_sections: true theme: flatly highlight: tango code_folding: show df_print: paged pdf_document: toc: true number_sections: true latex_engine: xelatex --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, warning = FALSE, message = FALSE, fig.width = 7, fig.height = 4.5, fig.align = "center" ) # Πακέτα που θα χρειαστούμε # install.packages(c("ggplot2", "deSolve")) # τρέξτε μία φορά αν δεν τα έχετε library(ggplot2) library(deSolve) theme_set(theme_minimal(base_size = 12)) ``` # Το πρόβλημα: μια λίμνη με νούφαρα Φανταστείτε ότι φτάνουμε σε μια λίμνη και μετράμε τα νούφαρα. - **Ημέρα 0:** $N_0 = 100$ νούφαρα - **Ημέρα 1:** $N_1 = 110$ νούφαρα Επιστρέφοντας την επόμενη μέρα παρατηρούμε **αύξηση 10%**. Ο στόχος μας είναι: > Πώς, ξεκινώντας από αυτή την απλή παρατήρηση, καταλήγουμε στον *τύπο της εκθετικής αύξησης* μέσω μιας διαφορικής εξίσωσης; Το ταξίδι έχει τρεις σταθμούς: 1. **Παρατήρηση → ρυθμός μεταβολής** (διακριτή σκέψη) 2. **Διακριτό → συνεχές** (παίρνουμε το όριο και γράφουμε τη διαφορική εξίσωση $\frac{dN}{dt} = rN$) 3. **Επίλυση της διαφορικής** (χωρισμός μεταβλητών → $N(t) = N_0 e^{rt}$) --- # Από την παρατήρηση στον ρυθμό μεταβολής ## Η μεταβολή του πληθυσμού Σε χρόνο $\Delta t = 1$ ημέρα, η μεταβολή του πληθυσμού είναι: $$ \Delta N \;=\; N_1 - N_0 \;=\; 110 - 100 \;=\; 10 $$ ## Ο **κατά κεφαλήν** (per capita) ρυθμός $r$ Δεν μας ενδιαφέρει η απόλυτη μεταβολή (10 νούφαρα), γιατί αυτή εξαρτάται από το πόσα νούφαρα είχαμε στην αρχή. Μας ενδιαφέρει **πόσο μεταβλήθηκε ο πληθυσμός ανά άτομο, ανά μονάδα χρόνου**: $$ r \;=\; \frac{\Delta N}{N \cdot \Delta t} \;=\; \frac{10}{100 \cdot 1} \;=\; 0{,}1 \;=\; 10\%/\text{ημέρα} $$ Το $r$ ονομάζεται **εγγενής (intrinsic) ρυθμός αύξησης** του πληθυσμού. Είναι η βασική παράμετρος όλης της κλασικής πληθυσμιακής δυναμικής. ```{r per-capita-rate} # Υπολογισμός του r από τα δεδομένα N0 <- 100 N1 <- 110 dt <- 1 # μέρες r <- (N1 - N0) / (N0 * dt) r ``` --- # Από το διακριτό στο συνεχές ## Διακριτή εξίσωση Αν υποθέσουμε ότι αυτός ο ρυθμός παραμένει σταθερός, τότε σε κάθε χρονικό βήμα $\Delta t$: $$ \Delta N \;=\; r \, N \, \Delta t \quad\Longleftrightarrow\quad \frac{\Delta N}{\Delta t} \;=\; r N $$ ### παράδειγμα τιμών για r=0.1 (10%) | N | ΔΝ/Δt | | :-------- | :-------- | | 100 | 10 | | 200 | 20 | | 1000 | 100 | | 5000 | 500 | Αυτή είναι μια **διακριτή** προσέγγιση: μετράμε τον πληθυσμό σε διακριτά βήματα (κάθε μέρα). ## Το όριο $\Delta t \to 0$ Όμως τα νούφαρα δεν περιμένουν την αυγή για να αυξηθούν — αυξάνονται *συνεχώς*. Παίρνοντας το όριο όταν το χρονικό βήμα γίνεται απειροελάχιστο: $$ \lim_{\Delta t \to 0} \frac{\Delta N}{\Delta t} \;=\; \frac{dN}{dt} $$ Καταλήγουμε στη **διαφορική εξίσωση της εκθετικής αύξησης**: $$ \boxed{\;\frac{dN}{dt} \;=\; r\,N\;} $$ > **Διαισθητική ερμηνεία:** Ο ρυθμός μεταβολής του πληθυσμού σε κάθε χρονική στιγμή είναι **ανάλογος** του τρέχοντος μεγέθους του πληθυσμού. Όσα περισσότερα νούφαρα, τόσο πιο γρήγορα εμφανίζονται καινούρια. --- # Επίλυση της διαφορικής εξίσωσης Θα λύσουμε τη $\frac{dN}{dt} = rN$ με τη μέθοδο του **χωρισμού μεταβλητών**. ## Βήμα 1 — Χωρίζουμε τις μεταβλητές $$ \frac{dN}{N} \;=\; r\,dt $$ (Όλα τα $N$ στο αριστερό μέλος, όλα τα $t$ στο δεξί.) ## Βήμα 2 — Ολοκληρώνουμε και τα δύο μέλη $$ \int \frac{dN}{N} \;=\; \int r\,dt \quad\Longrightarrow\quad \ln(N) \;=\; r\,t + C $$ όπου $C$ είναι σταθερά ολοκλήρωσης. ## Βήμα 3 — Λύνουμε ως προς $N$ $$ N(t) \;=\; e^{rt + C} \;=\; e^{C} \cdot e^{rt} $$ ## Βήμα 4 — Βρίσκουμε τη σταθερά από την αρχική συνθήκη Στο $t = 0$ έχουμε $N(0) = N_0$, οπότε: $$ N_0 \;=\; e^{C} \cdot e^{0} \;=\; e^{C} \quad\Longrightarrow\quad e^{C} = N_0 $$ ## Ο τελικός τύπος $$ \boxed{\;N(t) \;=\; N_0\, e^{r t}\;} $$ Αυτός είναι ο **τύπος της εκθετικής αύξησης**. --- # Υλοποίηση στην R ## Αναλυτική λύση Ας σχεδιάσουμε την $N(t) = 100 \cdot e^{0{,}1\,t}$ για 60 ημέρες: ```{r analytical-solution} N0 <- 100 r <- 0.1 t <- seq(0, 60, by = 0.1) N_analytical <- N0 * exp(r * t) df_analytical <- data.frame(t = t, N = N_analytical) ggplot(df_analytical, aes(x = t, y = N)) + geom_line(linewidth = 1, colour = "steelblue") + labs( title = "Εκθετική αύξηση πληθυσμού νούφαρων", subtitle = expression(paste(N(t) == N[0] * e^{r*t}, ", ", N[0] == 100, ", ", r == 0.1)), x = "Χρόνος (ημέρες)", y = "Πληθυσμός N(t)" ) ``` Σε 60 ημέρες ο πληθυσμός φτάνει τα `r round(N0 * exp(r * 60))` νούφαρα — μη ρεαλιστικό για μια πεπερασμένη λίμνη! Αυτό είναι το **πρόβλημα** της εκθετικής αύξησης που θα αντιμετωπίσουμε αργότερα με το λογιστικό μοντέλο. ## Αριθμητική λύση με `deSolve` Πολύ συχνά στη βιολογία οι διαφορικές εξισώσεις **δεν λύνονται αναλυτικά**. Τότε χρησιμοποιούμε αριθμητικούς ολοκληρωτές. Για το συγκεκριμένο μοντέλο, που το ξέρουμε ήδη αναλυτικά, αυτό μας δίνει μια ευκαιρία να ελέγξουμε αν η αριθμητική μέθοδος συμφωνεί. ```{r numerical-solution} # 1. Ορίζουμε τη διαφορική εξίσωση ως συνάρτηση exponential_model <- function(t, state, parameters) { with(as.list(c(state, parameters)), { dN <- r * N # dN/dt = rN list(dN) }) } # 2. Παράμετροι, αρχική συνθήκη, χρόνος parameters <- c(r = 0.1) state <- c(N = 100) times <- seq(0, 60, by = 0.1) # 3. Επίλυση out <- ode(y = state, times = times, func = exponential_model, parms = parameters) out_df <- as.data.frame(out) head(out_df) ``` ## Σύγκριση αναλυτικής και αριθμητικής λύσης ```{r comparison-plot} df_compare <- data.frame( t = times, Αναλυτική = N0 * exp(r * times), Αριθμητική = out_df$N ) ggplot(df_compare, aes(x = t)) + geom_line(aes(y = Αναλυτική), linewidth = 1.2, colour = "steelblue") + geom_point(aes(y = Αριθμητική), data = df_compare[seq(1, nrow(df_compare), by = 50), ], colour = "tomato", size = 2) + labs( title = "Αναλυτική (γραμμή) έναντι αριθμητικής (κουκίδες) λύσης", x = "Χρόνος (ημέρες)", y = "Πληθυσμός N(t)" ) # Μέγιστο σφάλμα max(abs(df_compare$Αναλυτική - df_compare$Αριθμητική)) ``` Το σφάλμα είναι αμελητέο — οι δύο λύσεις συμπίπτουν. --- # Μια λεπτή αλλά **πολύ σημαντική** παρατήρηση Παρατηρήστε κάτι: ξεκινήσαμε λέγοντας *"αύξηση 10% την ημέρα"*. Αυτή είναι μια **διακριτή** δήλωση. Όμως μετά γράψαμε $r = 0{,}1$ μέσα σε **συνεχή** διαφορική εξίσωση. Είναι τα δύο όντως ίσα; ## Ο διακριτός τύπος Αν κάθε μέρα ο πληθυσμός πολλαπλασιάζεται με το $1{,}10$: $$ N(t) \;=\; N_0 \cdot (1+r)^{t} \;=\; 100 \cdot 1{,}1^{t} $$ ## Σύγκριση των δύο μοντέλων ```{r discrete-vs-continuous} days <- 0:30 N_discrete <- 100 * (1 + 0.1)^days N_continuous <- 100 * exp(0.1 * days) df_dc <- data.frame( day = rep(days, 2), N = c(N_discrete, N_continuous), model = factor(rep(c("Διακριτό N0·(1+r)^t", "Συνεχές N0·exp(rt)"), each = length(days))) ) ggplot(df_dc, aes(x = day, y = N, colour = model)) + geom_line(linewidth = 1) + geom_point(size = 1.5) + labs( title = "Διακριτή έναντι συνεχούς εκθετικής αύξησης", subtitle = "Με την ίδια τιμή r = 0,1 τα δύο μοντέλα **δεν** ταυτίζονται", x = "Ημέρες", y = "Πληθυσμός", colour = "" ) + theme(legend.position = "bottom") ``` ```{r dc-table} data.frame( Ημέρα = c(1, 5, 10, 20, 30), Διακριτό = round(100 * 1.1^c(1, 5, 10, 20, 30), 1), Συνεχές = round(100 * exp(0.1 * c(1, 5, 10, 20, 30)), 1) ) ``` ## Γιατί διαφέρουν; Στο **διακριτό** μοντέλο, "10% αύξηση/ημέρα" σημαίνει: στο τέλος της μέρας ο πληθυσμός είναι 110% του αρχικού. Στο **συνεχές** μοντέλο με $r = 0{,}1$, ο *στιγμιαίος* ρυθμός είναι 10%, αλλά επειδή τα νεογέννητα νούφαρα **παράγουν αμέσως κι αυτά νέα νούφαρα** (συνεχής ανατοκισμός!), στο τέλος της μέρας έχουμε: $$ N(1) \;=\; 100 \cdot e^{0{,}1} \;\approx\; 110{,}52 $$ δηλαδή λίγο παραπάνω από 10% αύξηση. ## Πώς συνδέονται; Αν θέλουμε το συνεχές μοντέλο να δίνει **ακριβώς** 10% αύξηση ανά μέρα, χρειαζόμαστε: $$ e^{r_{\text{cont}}} = 1{,}1 \quad\Longrightarrow\quad r_{\text{cont}} \;=\; \ln(1{,}1) \;\approx\; 0{,}0953 $$ ```{r r-conversion} log(1.1) # στο R, log() είναι ο φυσικός λογάριθμος (ln) ``` > **Στην οικολογική πράξη**, τις περισσότερες φορές δουλεύουμε στο συνεχές μοντέλο και θεωρούμε $r$ τον στιγμιαίο ρυθμό. Η διαφορά είναι μικρή για μικρά $r$, αλλά καλό είναι να ξέρετε ότι **υπάρχει**. --- # Λογαριθμικός άξονας: το διαγνωστικό εργαλείο του βιολόγου Παίρνοντας λογάριθμο στον τύπο $N(t) = N_0 e^{rt}$: $$ \ln N(t) \;=\; \ln N_0 \;+\; r\,t $$ Δηλαδή ο $\ln N$ είναι **γραμμική** συνάρτηση του $t$ με κλίση ίση με $r$. Αυτό είναι εξαιρετικά χρήσιμο πρακτικά: > **Αν τα δεδομένα σας σε ημι-λογαριθμικό διάγραμμα ($\ln N$ ως προς $t$) σχηματίζουν ευθεία, τότε ο πληθυσμός αυξάνεται εκθετικά**, και η κλίση της ευθείας σας δίνει το $r$. ```{r log-plot} ggplot(df_analytical, aes(x = t, y = N)) + geom_line(linewidth = 1, colour = "steelblue") + scale_y_log10() + labs( title = "Εκθετική αύξηση σε ημι-λογαριθμικό άξονα", subtitle = "Η ευθεία γραμμή υποδηλώνει εκθετική αύξηση", x = "Χρόνος (ημέρες)", y = "log10(N(t))" ) ``` ## Εκτίμηση του $r$ από δεδομένα με γραμμική παλινδρόμηση ```{r estimate-r} # Προσομοιωμένα δεδομένα με λίγο θόρυβο set.seed(42) sim_days <- 0:30 sim_N <- 100 * exp(0.1 * sim_days) * exp(rnorm(length(sim_days), 0, 0.05)) # Γραμμική παλινδρόμηση log(N) ~ t fit <- lm(log(sim_N) ~ sim_days) summary(fit)$coefficients cat("Εκτιμώμενο r =", round(coef(fit)[2], 4), " (πραγματικό r = 0.1)\n") cat("Εκτιμώμενο N0 =", round(exp(coef(fit)[1]), 2), " (πραγματικό N0 = 100)\n") ``` --- # Χρόνος διπλασιασμού Μια χρήσιμη ποσότητα: σε πόσο χρόνο διπλασιάζεται ο πληθυσμός; $$ 2 N_0 = N_0 e^{r T_2} \;\Longleftrightarrow\; T_2 = \frac{\ln 2}{r} $$ ```{r doubling-time} T2 <- log(2) / r cat("Χρόνος διπλασιασμού των νούφαρων:", round(T2, 2), "ημέρες\n") ``` Άρα τα νούφαρα διπλασιάζονται περίπου κάθε **7 μέρες**. --- # Ασκήσεις 1. **Παίξτε με το $r$.** Σχεδιάστε στο ίδιο διάγραμμα την εκθετική αύξηση για $r = 0{,}05$, $r = 0{,}1$ και $r = 0{,}2$. Πώς αλλάζει ο χρόνος διπλασιασμού; 2. **Αρνητικό $r$.** Τι συμβαίνει αν $r < 0$; Σχεδιάστε την περίπτωση $r = -0{,}05$ και ερμηνεύστε βιολογικά. 3. **Εκτίμηση από πραγματικά δεδομένα.** Σας δίνονται οι παρακάτω μετρήσεις: ```r day <- c(0, 3, 7, 14, 21, 28) N <- c(50, 73, 109, 218, 437, 870) ``` Εκτιμήστε το $r$ με γραμμική παλινδρόμηση στο $\ln N$ και προβλέψτε τον πληθυσμό στις 60 ημέρες. 4. **Ο μαθηματικός κίνδυνος.** Με $r = 0{,}1$ και $N_0 = 100$ νούφαρα, σε πόσες μέρες τα νούφαρα θα ξεπεράσουν τον αριθμό των ατόμων στη Γη ($\approx 8 \times 10^{9}$); Τι μας λέει αυτό για τους *περιορισμούς* του εκθετικού μοντέλου; --- # Επόμενο βήμα: ο νόμος επιστρέφει στην πραγματικότητα Το εκθετικό μοντέλο **δεν μπορεί** να ισχύει για πάντα. Σε μια πραγματική λίμνη υπάρχουν περιορισμοί: χώρος, φως, θρεπτικά συστατικά. Σε κάποιο σημείο ο πληθυσμός θα *κορεστεί*. Στο επόμενο μάθημα θα τροποποιήσουμε τη διαφορική εξίσωση ώστε να συμπεριλάβει τη **φέρουσα ικανότητα** $K$ του περιβάλλοντος: $$ \frac{dN}{dt} \;=\; r N \left(1 - \frac{N}{K}\right) $$ Αυτή είναι η **λογιστική εξίσωση** — και είναι το θέμα του επόμενου εργαστηρίου. --- # Στοιχεία συνεδρίας R ```{r session-info} sessionInfo() ```