1 Πού μείναμε

Στο προηγούμενο εργαστήριο καταλήξαμε στον τύπο της εκθετικής αύξησης:

\[ \frac{dN}{dt} = r N \quad\Longrightarrow\quad N(t) = N_0\, e^{rt} \]

Με \(N_0 = 100\) νούφαρα και \(r = 0{,}1\)/ημέρα, ας δούμε τι προβλέπει το μοντέλο σε βάθος χρόνου:

N0 <- 100
r  <- 0.1

cat("Σε 100 ημέρες:", format(N0 * exp(r * 100),  big.mark = "."), "νούφαρα\n")
## Σε 100 ημέρες: 2.202.647 νούφαρα
cat("Σε 200 ημέρες:", format(N0 * exp(r * 200),  big.mark = "."), "νούφαρα\n")
## Σε 200 ημέρες: 48.516.519.541 νούφαρα
cat("Σε 365 ημέρες:", format(N0 * exp(r * 365),  big.mark = "."), "νούφαρα\n")
## Σε 365 ημέρες: 7.108019e+17 νούφαρα

Είναι παράλογο. Σε λιγότερο από έναν χρόνο, η μικρή μας λίμνη υποτίθεται ότι θα έχει περισσότερα νούφαρα από όσα κύτταρα περιέχει το ανθρώπινο σώμα.

Το πρόβλημα είναι ότι το εκθετικό μοντέλο αγνοεί τον χώρο, το φως, τα θρεπτικά συστατικά και τον ανταγωνισμό. Σε αυτό το εργαστήριο θα διορθώσουμε αυτή την έλλειψη.

2 Η ιδέα της φέρουσας ικανότητας \(K\)

Κάθε πραγματικό περιβάλλον έχει ένα όριο στον αριθμό ατόμων που μπορεί να συντηρήσει. Ονομάζουμε αυτό το όριο φέρουσα ικανότητα (carrying capacity) και το συμβολίζουμε με \(K\).

Για τη λίμνη μας, ας υποθέσουμε:

K <- 10000   # μέγιστος αριθμός νούφαρων που χωρά η λίμνη

Σημείωση επιλογής: Η τιμή \(K = 10{.}000\) είναι ενδεικτική και επιλέχθηκε για λόγους σαφήνειας στα διαγράμματα. Σε πραγματική μελέτη το \(K\) εκτιμάται από δεδομένα ή από τα φυσικά χαρακτηριστικά του οικοτόπου (επιφάνεια λίμνης, μέσο μέγεθος νούφαρου, διαθέσιμο φως κλπ).

2.1 Το \(K\) δεν είναι απλώς “ταβάνι”

Η αφελής ιδέα θα ήταν: “αφήνουμε τον πληθυσμό να αυξάνεται εκθετικά μέχρι το \(K\), και μετά τον σταματάμε.” Αυτό όμως δεν συμβαίνει στη φύση — η αύξηση επιβραδύνεται σταδιακά καθώς πλησιάζουμε στο \(K\), πολύ πριν το φτάσουμε.

Γιατί;

3 Ένας μηχανισμός: γιατί το \(K\) μειώνει τον ρυθμό

Φανταστείτε τη λίμνη μας με τα νούφαρα. Όταν ο πληθυσμός είναι μικρός, κάθε νούφαρο έχει άπλετο χώρο γύρω του και μπορεί να αναπαραχθεί ελεύθερα. Όμως καθώς ο πληθυσμός μεγαλώνει, τα νούφαρα στο εσωτερικό της συστάδας περιβάλλονται από άλλα νούφαρα και δυσκολεύονται να επεκταθούν. Μόνο τα νούφαρα στην περιφέρεια έχουν χώρο. Όσο πιο γεμάτη η λίμνη, τόσο μικρότερο το ποσοστό των νούφαρων που μπορούν πραγματικά να αναπαραχθούν.

Αυτό είναι ένα παράδειγμα μηχανισμού. Στην πραγματικότητα μπορούν να συντρέχουν πολλοί λόγοι:

  • Χώρος / φως: τα νούφαρα στο εσωτερικό σκιάζονται.
  • Θρεπτικά: ο ανταγωνισμός για ριζικά θρεπτικά συστατικά εντείνεται.
  • Παθογόνα: σε πυκνές συστάδες οι ασθένειες εξαπλώνονται γρήγορα.
  • Τοξικά παραπροϊόντα: συσσωρεύονται όταν ο πληθυσμός πυκνώνει.

Όλα αυτά μοιράζονται το ίδιο ποιοτικό χαρακτηριστικό: ο κατά κεφαλήν ρυθμός αναπαραγωγής μειώνεται καθώς ο πληθυσμός πλησιάζει το \(K\).

Άρα το \(K\) δεν είναι απλώς ένα όριο που χτυπάμε επάνω. Το \(K\) τροποποιεί τον ίδιο τον ρυθμό μεταβολής.

4 Από τη φυσική διαίσθηση στη μαθηματική μορφή

4.1 Τι θέλουμε από τον “πραγματοποιημένο” κατά κεφαλήν ρυθμό

Ορίζουμε τον πραγματοποιημένο (realized) κατά κεφαλήν ρυθμό αύξησης ως μια συνάρτηση \(r_{\text{eff}}(N)\) του τρέχοντος πληθυσμού. Δύο φυσικές απαιτήσεις:

  1. Όταν \(N\) είναι μικρό, υπάρχει άπλετος χώρος → \(r_{\text{eff}}(N) \to r\) (ο εγγενής ρυθμός).
  2. Όταν \(N \to K\), ο χώρος εξαντλείται → \(r_{\text{eff}}(N) \to 0\).

4.2 Η απλούστερη συνάρτηση που ικανοποιεί αυτά: γραμμική

Η πιο απλή συνάρτηση που πιάνει τα δύο άκρα είναι μια ευθεία γραμμή που πέφτει από \(r\) (στο \(N=0\)) στο \(0\) (στο \(N=K\)):

\[ \boxed{\;r_{\text{eff}}(N) \;=\; r\left(1 - \frac{N}{K}\right)\;} \]

Σημείωση επιλογής μοντέλου (σημαντική): Η γραμμική μορφή είναι μια μοντελοποιητική παραδοχή, όχι λογική αναγκαιότητα. Είναι η απλούστερη συνάρτηση που πιάνει τη σωστή ποιοτική συμπεριφορά. Άλλες δυνατές μορφές (θα τις δούμε στο τέλος):

  • θ-λογιστική: \(r_{\text{eff}}(N) = r\bigl(1 - (N/K)^{\theta}\bigr)\)
  • Gompertz: \(r_{\text{eff}}(N) = r\,\ln(K/N)\)

Όλες δίνουν παρόμοια σιγμοειδή συμπεριφορά, αλλά με διαφορετική γεωμετρία.

N_grid <- seq(0, K, length.out = 200)
df_pc <- data.frame(
  N        = rep(N_grid, 2),
  rate     = c(rep(r, length(N_grid)),                # σταθερό r (εκθετικό)
               r * (1 - N_grid/K)),                   # γραμμικά μειούμενο (λογιστικό)
  Μοντέλο  = factor(rep(c("Εκθετικό:  r σταθερό",
                          "Λογιστικό: r·(1-N/K)"),
                        each = length(N_grid)))
)

ggplot(df_pc, aes(x = N, y = rate, colour = Μοντέλο)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  labs(
    title = "Κατά κεφαλήν ρυθμός αύξησης ως συνάρτηση του πληθυσμού",
    subtitle = "Η ουσιαστική διαφορά μεταξύ εκθετικού και λογιστικού μοντέλου",
    x = "Πληθυσμός N",
    y = expression(r[eff](N))
  ) +
  theme(legend.position = "bottom")

5 Η λογιστική εξίσωση

Ο συνολικός ρυθμός μεταβολής του πληθυσμού είναι:

\[ \frac{dN}{dt} \;=\; r_{\text{eff}}(N) \cdot N \;=\; r\,N\left(1 - \frac{N}{K}\right) \]

Αυτή είναι η λογιστική εξίσωση (Verhulst, 1838).

\[ \boxed{\;\frac{dN}{dt} \;=\; r\,N\left(1 - \frac{N}{K}\right)\;} \]

6 Ποιοτική ανάλυση πριν λύσουμε

Ένα από τα πιο χρήσιμα εργαλεία στα δυναμικά συστήματα είναι να εξετάζουμε τη συμπεριφορά πριν λύσουμε αναλυτικά. Πολύ συχνά μάλιστα η λύση δεν υπάρχει αναλυτικά — αλλά η ποιοτική εικόνα μάς αρκεί.

6.1 Σημεία ισορροπίας

Όπου \(\dfrac{dN}{dt} = 0\) ο πληθυσμός μένει στάσιμος. Λύνουμε:

\[ r\,N\left(1 - \frac{N}{K}\right) = 0 \quad\Longrightarrow\quad N^* = 0 \;\text{ ή }\; N^* = K \]

  • \(N^* = 0\): άδεια λίμνη (πληθυσμός εξαφανισμένος).
  • \(N^* = K\): γεμάτη λίμνη (κορεσμός).

6.2 Ευστάθεια σημείων ισορροπίας

  • Αν ξεκινήσουμε κοντά στο \(0\) (π.χ. \(N=1\)): \(\frac{dN}{dt} > 0\), ο πληθυσμός μεγαλώνει — απομακρύνεται από το 0. Άρα \(N^*=0\) είναι ασταθές.
  • Αν ξεκινήσουμε κοντά στο \(K\) (π.χ. \(N=0{,}9K\)): \(\frac{dN}{dt} > 0\), ο πληθυσμός πλησιάζει το \(K\).
  • Αν ξεκινήσουμε πάνω από \(K\) (π.χ. \(N=1{,}5K\)): \(\frac{dN}{dt} < 0\), ο πληθυσμός μειώνεται προς το \(K\).

Άρα \(N^*=K\) είναι ευσταθές — είναι ένας ελκυστής του συστήματος.

6.3 Φασικό διάγραμμα: \(\dfrac{dN}{dt}\) ως προς \(N\)

N_grid <- seq(-500, K * 1.4, length.out = 300)
dN_dt  <- r * N_grid * (1 - N_grid/K)
df_phase <- data.frame(N = N_grid, dN = dN_dt)

# Σημεία ισορροπίας
eq_unstable <- data.frame(N = 0, dN = 0)
eq_stable   <- data.frame(N = K, dN = 0)

# Βέλη ροής στον άξονα N: το πρόσημο του (xend - x) δίνει την κατεύθυνση κίνησης
arrows_df <- data.frame(
  x    = c(K/4,        3*K/4,        1.2*K),
  xend = c(K/4 + 800,  3*K/4 + 800,  1.2*K - 800),
  y    = 0, yend = 0
)

ggplot(df_phase, aes(x = N, y = dN)) +
  geom_hline(yintercept = 0, colour = "grey50") +
  geom_line(linewidth = 1, colour = "darkgreen") +
  geom_point(data = eq_unstable, aes(x = N, y = dN),
             shape = 21, fill = "white", colour = "darkgreen",
             size = 4, stroke = 1.2) +
  geom_point(data = eq_stable, aes(x = N, y = dN),
             colour = "darkgreen", size = 4) +
  geom_segment(data = arrows_df,
               aes(x = x, xend = xend, y = y, yend = yend),
               arrow = arrow(length = unit(0.25, "cm")),
               colour = "tomato", linewidth = 0.9,
               inherit.aes = FALSE) +
  annotate("text", x = 0,    y = -150, label = "N* = 0\n(ασταθές)",  vjust = 1) +
  annotate("text", x = K,    y = -150, label = "N* = K\n(ευσταθές)", vjust = 1) +
  annotate("text", x = K/2,  y = max(dN_dt) + 30,
           label = "μέγιστος ρυθμός\nστο N = K/2", vjust = 0) +
  labs(
    title = "Φασικό διάγραμμα της λογιστικής εξίσωσης",
    subtitle = "Τα βέλη δείχνουν την κατεύθυνση μεταβολής του N",
    x = "N", y = "dN/dt = r"
  )

Η παραβολή αναπαριστά όλη τη δυναμική του συστήματος σε μια εικόνα:

  • Όπου η καμπύλη είναι πάνω από τον άξονα → \(\frac{dN}{dt} > 0\) → ο πληθυσμός αυξάνεται.
  • Όπου η καμπύλη είναι κάτω από τον άξονα → ο πληθυσμός μειώνεται.
  • Στις ρίζες έχουμε σημεία ισορροπίας.
  • Η κορυφή της παραβολής (στο \(N = K/2\)) είναι ο μέγιστος ρυθμός αύξησης — εκεί ο πληθυσμός μεγαλώνει πιο γρήγορα. Αυτό είναι και το σημείο καμπής της χρονικής τροχιάς.

7 Αναλυτική λύση (προαιρετικό για μαθηματικά διατεθειμένους)

Για όσους θέλουν να δουν τη λύση: ξεκινάμε με χωρισμό μεταβλητών:

\[ \frac{dN}{N\,(1 - N/K)} = r\,dt \quad\Longleftrightarrow\quad \frac{K\,dN}{N\,(K - N)} = r\,dt \]

Παραγοντική ανάπτυξη (partial fractions):

\[ \frac{K}{N(K-N)} = \frac{1}{N} + \frac{1}{K-N} \]

Άρα:

\[ \int\!\left(\frac{1}{N} + \frac{1}{K-N}\right) dN = \int r\,dt \;\Longleftrightarrow\; \ln N - \ln(K-N) = rt + C \]

\[ \ln\!\left(\frac{N}{K-N}\right) = rt + C \;\Longrightarrow\; \frac{N}{K-N} = A\,e^{rt}, \quad A = \frac{N_0}{K - N_0} \]

Λύνοντας ως προς \(N\) καταλήγουμε στον κλειστό τύπο της λογιστικής λύσης:

\[ \boxed{\;N(t) \;=\; \dfrac{K}{1 + \left(\dfrac{K - N_0}{N_0}\right) e^{-rt}}\;} \]

Έλεγχος (sanity checks):

  • Στο \(t = 0\): \(N(0) = K \big/ \left(1 + \tfrac{K-N_0}{N_0}\right) = K \cdot \tfrac{N_0}{K} = N_0\). ✓
  • Στο \(t \to \infty\): το \(e^{-rt} \to 0\), άρα \(N \to K\). ✓
  • Αν \(N_0 \ll K\), για μικρά \(t\) ο όρος \(\tfrac{K-N_0}{N_0} \approx \tfrac{K}{N_0}\) είναι τεράστιος, και η συμπεριφορά είναι σχεδόν εκθετική (γιατί;). ✓

8 R: αριθμητική λύση με deSolve

Τώρα που ξέρουμε τι περιμένουμε, ας προσομοιώσουμε:

logistic_model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N * (1 - N / K)
    list(dN)
  })
}

parameters <- c(r = 0.1, K = 10000)
state      <- c(N = 100)
times      <- seq(0, 200, by = 0.5)

out_log <- as.data.frame(
  ode(y = state, times = times, func = logistic_model, parms = parameters)
)

head(out_log)

8.1 Σύγκριση αναλυτικής και αριθμητικής λύσης

N0 <- 100
K  <- 10000
r  <- 0.1

logistic_analytical <- function(t, N0, r, K) {
  K / (1 + ((K - N0) / N0) * exp(-r * t))
}

df_log <- data.frame(
  t          = times,
  Αναλυτική  = logistic_analytical(times, N0, r, K),
  Αριθμητική = out_log$N
)

ggplot(df_log, aes(x = t)) +
  geom_line(aes(y = Αναλυτική),  linewidth = 1.2, colour = "steelblue") +
  geom_point(aes(y = Αριθμητική),
             data = df_log[seq(1, nrow(df_log), by = 25), ],
             colour = "tomato", size = 1.8) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95, y = K * 1.04,
           label = "K = 10.000", colour = "grey40", hjust = 1) +
  labs(
    title = "Λογιστική αύξηση: αναλυτική (γραμμή) vs αριθμητική (κουκίδες)",
    x = "Χρόνος (ημέρες)",
    y = "Πληθυσμός N(t)"
  )

cat("Μέγιστο σφάλμα:", format(max(abs(df_log$Αναλυτική - df_log$Αριθμητική)),
                              scientific = TRUE), "\n")
## Μέγιστο σφάλμα: 6.769672e-04

9 Η σιγμοειδής καμπύλη και τα τρία στάδια

Παρατηρήστε προσεκτικά τη χρονική τροχιά. Έχει τρία διακριτά στάδια:

stages <- data.frame(
  t = times,
  N = logistic_analytical(times, N0, r, K)
)

# Σημείο καμπής: N = K/2
t_inflection <- (1/r) * log((K - N0) / N0)

ggplot(stages, aes(x = t, y = N)) +
  geom_line(linewidth = 1.1, colour = "steelblue") +
  geom_hline(yintercept = K,   linetype = "dashed", colour = "grey50") +
  geom_hline(yintercept = K/2, linetype = "dotted", colour = "grey50") +
  geom_vline(xintercept = t_inflection, linetype = "dotted", colour = "tomato") +
  annotate("point", x = t_inflection, y = K/2, size = 3, colour = "tomato") +
  annotate("text",  x = t_inflection + 5, y = K/2,
           label = "Σημείο καμπής (N=K/2)\nμέγιστος ρυθμός", hjust = 0) +
  annotate("text", x =  20, y = 1500, label = "1. Σχεδόν εκθετικό",     hjust = 0) +
  annotate("text", x =  60, y = 5500, label = "2. Μεταβατική φάση\n(επιτάχυνση → επιβράδυνση)", hjust = 0) +
  annotate("text", x = 130, y = 9300, label = "3. Κορεσμός προς K",       hjust = 0) +
  labs(
    title = "Η σιγμοειδής (S-shaped) λογιστική καμπύλη",
    x = "Χρόνος (ημέρες)",
    y = "N(t)"
  )

cat("Σημείο καμπής στο t =", round(t_inflection, 1), "ημέρες (όπου N = K/2 =", K/2, ")\n")
## Σημείο καμπής στο t = 46 ημέρες (όπου N = K/2 = 5000 )
  1. Σχεδόν εκθετική φάση (μικρά \(N\)): ο όρος \(1 - N/K\) είναι κοντά στο 1, άρα \(\frac{dN}{dt} \approx rN\).
  2. Σημείο καμπής (\(N = K/2\)): εδώ ο ρυθμός \(\frac{dN}{dt}\) μεγιστοποιείται.
  3. Φάση κορεσμού (\(N \to K\)): ο όρος \(1 - N/K \to 0\), άρα ο ρυθμός φθίνει εκθετικά προς το μηδέν.

10 Σύγκριση: εκθετικό vs λογιστικό

Με ίδιο \(r\), ίδιο \(N_0\), αλλά διαφορετικό \(K\):

df_compare <- data.frame(
  t          = times,
  Εκθετικό   = N0 * exp(r * times),
  Λογιστικό  = logistic_analytical(times, N0, r, K)
)

# Long format για ggplot
df_long <- data.frame(
  t       = rep(times, 2),
  N       = c(df_compare$Εκθετικό, df_compare$Λογιστικό),
  Μοντέλο = factor(rep(c("Εκθετικό", "Λογιστικό"), each = length(times)))
)

ggplot(df_long, aes(x = t, y = N, colour = Μοντέλο)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  scale_y_log10() +
  annotate("text", x = max(times) * 0.95, y = K * 1.3,
           label = "K = 10.000", colour = "grey40", hjust = 1) +
  labs(
    title = "Εκθετική vs λογιστική αύξηση (λογαριθμικός άξονας)",
    subtitle = "Στις πρώτες ημέρες είναι σχεδόν ίδιες — μετά αποκλίνουν δραματικά",
    x = "Χρόνος (ημέρες)",
    y = "N(t) [λογαριθμική κλίμακα]"
  ) +
  theme(legend.position = "bottom")

Στις πρώτες \(\approx 30\) ημέρες οι δύο καμπύλες είναι σχεδόν αξεχώριστες. Αυτό μας δίνει μια προειδοποίηση: στα πρώιμα στάδια ενός πληθυσμιακού φαινομένου, δεν μπορούμε εύκολα να ξεχωρίσουμε αν αυξάνεται εκθετικά ή λογιστικά — και τα δύο μοντέλα ταιριάζουν εξίσου καλά. Η διαφορά αναδύεται μόνο όταν πλησιάζουμε στο \(K\).

11 Τι συμβαίνει αν ξεκινήσουμε από πάνω από το \(K\);

times2 <- seq(0, 100, by = 0.5)

scenarios <- data.frame(
  t   = rep(times2, 3),
  N   = c(logistic_analytical(times2, N0 = 100,   r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 5000,  r = 0.1, K = 10000),
          logistic_analytical(times2, N0 = 15000, r = 0.1, K = 10000)),
  N0  = factor(rep(c("N0 = 100",  "N0 = 5.000",  "N0 = 15.000 (overshoot)"),
                   each = length(times2)))
)

ggplot(scenarios, aes(x = t, y = N, colour = N0)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = K, linetype = "dashed", colour = "grey50") +
  annotate("text", x = max(times2) * 0.95, y = K * 1.04,
           label = "K = 10.000", colour = "grey40", hjust = 1) +
  labs(
    title = "Λογιστική αύξηση: τρία διαφορετικά αρχικά σενάρια",
    subtitle = "Όλα συγκλίνουν στο K, ανεξάρτητα από την αρχική συνθήκη",
    x = "Χρόνος (ημέρες)",
    y = "N(t)"
  ) +
  theme(legend.position = "bottom")

Παρατηρήστε ότι η περίπτωση \(N_0 = 15.000 > K\) μειώνεται προς το \(K\). Αυτό βγαίνει αυτόματα από την εξίσωση: αν \(N > K\), τότε \(1 - N/K < 0\), άρα \(\frac{dN}{dt} < 0\).

Αυτό μοντελοποιεί π.χ. την περίπτωση που εισάγουμε υπερβολικό αριθμό νούφαρων σε μια λίμνη — ο πληθυσμός θα πέσει στο φέρον επίπεδο.

12 Εκτίμηση παραμέτρων από δεδομένα

Στην πράξη οι βιολόγοι έχουν δεδομένα και θέλουν να εκτιμήσουν τα \(r\) και \(K\). Αυτό γίνεται με μη-γραμμική παλινδρόμηση (nls() στην R).

# Προσομοιώνουμε δεδομένα με θόρυβο
set.seed(123)
sample_days <- seq(0, 150, by = 10)
true_N      <- logistic_analytical(sample_days, N0 = 100, r = 0.1, K = 10000)
observed_N  <- true_N * exp(rnorm(length(sample_days), 0, 0.08))

# Προσαρμογή του λογιστικού μοντέλου
fit <- nls(observed_N ~ K / (1 + ((K - N0)/N0) * exp(-r * sample_days)),
           start = list(K = 8000, r = 0.05, N0 = 200))

summary(fit)$coefficients
##        Estimate   Std. Error   t value     Pr(>|t|)
## K  1.023666e+04 241.96331612 42.306681 2.599395e-15
## r  9.074180e-02   0.01263879  7.179627 7.158844e-06
## N0 1.601025e+02  91.19055361  1.755692 1.026589e-01
fit_curve <- data.frame(
  t = times,
  N = predict(fit, newdata = data.frame(sample_days = times))
)

ggplot() +
  geom_point(data = data.frame(t = sample_days, N = observed_N),
             aes(x = t, y = N), size = 2.5, colour = "tomato") +
  geom_line(data = fit_curve, aes(x = t, y = N),
            linewidth = 1, colour = "steelblue") +
  geom_hline(yintercept = coef(fit)["K"],
             linetype = "dashed", colour = "grey40") +
  annotate("text", x = max(times) * 0.95,
           y = coef(fit)["K"] * 1.05,
           label = paste0("Εκτίμηση K = ", round(coef(fit)["K"])),
           colour = "grey30", hjust = 1) +
  labs(
    title = "Προσαρμογή λογιστικού μοντέλου σε δεδομένα",
    subtitle = paste0("Εκτιμήσεις: r = ", round(coef(fit)["r"], 3),
                      ", K = ",  round(coef(fit)["K"]),
                      ", N0 = ", round(coef(fit)["N0"])),
    x = "Χρόνος (ημέρες)",
    y = "Πληθυσμός N(t)"
  )

Πρακτική σημασία: Τα nls() εκτιμώμενα \(\hat{r}\) και \(\hat{K}\) έρχονται με τυπικό σφάλμα — άρα έχουμε στατιστικά διαστήματα εμπιστοσύνης για τις παραμέτρους. Αυτό είναι κρίσιμο για να συγκρίνουμε δύο πληθυσμούς ή δύο πειραματικές συνθήκες.

13 Πέρα από τη λογιστική: άλλες λειτουργικές μορφές

Όπως αναφέραμε, η γραμμική μορφή \(r_{\text{eff}}(N) = r(1 - N/K)\) είναι μία επιλογή. Δύο σημαντικές εναλλακτικές:

13.1 θ-λογιστική

\[ \frac{dN}{dt} = r N \left(1 - \left(\frac{N}{K}\right)^{\theta}\right) \]

  • \(\theta = 1\) → κλασική λογιστική.
  • \(\theta > 1\) → ο πληθυσμός μένει κοντά στο \(K\) για περισσότερο χρόνο πριν επιβραδύνει (π.χ. \(r\)-στρατηγικές).
  • \(\theta < 1\) → η επιβράδυνση ξεκινά νωρίς (π.χ. ευαίσθητοι σε ανταγωνισμό).

13.2 Gompertz

\[ \frac{dN}{dt} = r N \ln\!\left(\frac{K}{N}\right) \]

Χρησιμοποιείται κυρίως σε αύξηση όγκων (π.χ. καρκινικοί όγκοι, ανάπτυξη οργανισμών). Επιβραδύνει νωρίτερα από τη λογιστική.

gompertz <- function(t, N0, r, K) {
  K * exp(log(N0/K) * exp(-r * t))
}
theta_logistic <- function(t, N0, r, K, theta) {
  # Αριθμητική λύση
  parms <- c(r = r, K = K, theta = theta)
  func  <- function(t, state, parms) {
    with(as.list(c(state, parms)),
         list(r * N * (1 - (N/K)^theta)))
  }
  out <- ode(y = c(N = N0), times = t, func = func, parms = parms)
  as.data.frame(out)$N
}

t_vec <- seq(0, 200, by = 0.5)
df_alt <- data.frame(
  t        = rep(t_vec, 4),
  N        = c(logistic_analytical(t_vec, 100, 0.1, 10000),
               gompertz(t_vec,             100, 0.1, 10000),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 0.5),
               theta_logistic(t_vec,       100, 0.1, 10000, theta = 3)),
  Μοντέλο  = factor(rep(c("Λογιστικό (θ=1)",
                          "Gompertz",
                          "θ-λογιστικό (θ=0.5)",
                          "θ-λογιστικό (θ=3)"),
                        each = length(t_vec)),
                    levels = c("Λογιστικό (θ=1)",
                               "θ-λογιστικό (θ=0.5)",
                               "θ-λογιστικό (θ=3)",
                               "Gompertz"))
)

ggplot(df_alt, aes(x = t, y = N, colour = Μοντέλο)) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 10000, linetype = "dashed", colour = "grey50") +
  labs(
    title    = "Διαφορετικές μοντελοποιήσεις του κορεσμού",
    subtitle = "Όλα τα μοντέλα έχουν ίδιο r=0.1, K=10.000, N0=100",
    x = "Χρόνος (ημέρες)", y = "N(t)"
  ) +
  theme(legend.position = "bottom")

Παιδαγωγικό σχόλιο: Το ότι όλα αυτά τα μοντέλα δίνουν διαφορετικές καμπύλες είναι το πιο σημαντικό μήνυμα. Σε πραγματικά δεδομένα η επιλογή του μοντέλου δεν προκύπτει αυτόματα — εξαρτάται από τη βιολογία του οργανισμού και πρέπει να ελεγχθεί στατιστικά (π.χ. AIC).

14 Ασκήσεις

  1. Ευαισθησία στο \(K\). Σχεδιάστε λογιστικές καμπύλες με \(r=0{,}1\), \(N_0=100\) και \(K \in \{1.000,\ 5.000,\ 10.000,\ 50.000\}\) στο ίδιο διάγραμμα. Τι αλλάζει στη μορφή της καμπύλης;

  2. Χρόνος για να φτάσουμε στο \(90\%\) του \(K\). Λύστε αναλυτικά: για ποιο \(t\) ισχύει \(N(t) = 0{,}9 K\); Επαληθεύστε αριθμητικά.

  3. Σημείο καμπής. Δείξτε ότι ο ρυθμός \(\frac{dN}{dt}\) μεγιστοποιείται στο \(N = K/2\). (Υπόδειξη: παραγωγίστε τον \(\frac{dN}{dt}\) ως προς \(N\) και βρείτε τη ρίζα.)

  4. Από τα δεδομένα στο \(r\). Σας δίνονται μετρήσεις:

    day <- c(0, 10, 20, 40, 60, 90, 120, 180)
N   <- c(120, 280, 620, 2400, 5800, 9100, 9750, 9970)

    Προσαρμόστε το λογιστικό μοντέλο με nls() και υπολογίστε τα διαστήματα εμπιστοσύνης για τα \(\hat{r}\) και \(\hat{K}\) (confint(fit)).

  5. Εννοιολογική. Στη συζήτηση αναφέραμε διάφορους μηχανισμούς που μπορούν να δικαιολογήσουν τη μείωση του κατά κεφαλήν ρυθμού (χώρος, φως, θρεπτικά, παθογόνα). Για τα νούφαρα στη λίμνη μας, ποιον μηχανισμό θεωρείτε πιο σημαντικό και γιατί; Πώς θα τον επαληθεύατε πειραματικά;

  6. Πέρα από τη μέση τιμή. Όλα τα μοντέλα εδώ είναι ντετερμινιστικά — δεν υπάρχει τύχη. Σε πραγματικούς πληθυσμούς όμως υπάρχει στοχαστικότητα (γεννήσεις/θάνατοι είναι τυχαία γεγονότα). Δοκιμάστε να προσομοιώσετε στοχαστική λογιστική αύξηση: σε κάθε χρονικό βήμα, προσθέστε στον \(\frac{dN}{dt}\) έναν τυχαίο όρο rnorm(1, 0, σ). Πώς αλλάζει η συμπεριφορά για διαφορετικά \(\sigma\);

15 Επόμενο βήμα

Μέχρι τώρα έχουμε εξετάσει έναν μοναδικό πληθυσμό. Στην πραγματικότητα όμως οι πληθυσμοί αλληλεπιδρούν μεταξύ τους:

  • Θηρευτής–θήραμα: εξισώσεις Lotka–Volterra
  • Ανταγωνισμός μεταξύ ειδών: επεκτεταμένη λογιστική με δύο ή περισσότερες εξισώσεις
  • Συμβίωση: αμοιβαία ωφέλιμες αλληλεπιδράσεις

Αυτά είναι συστήματα διαφορικών εξισώσεων και θα είναι το θέμα του επόμενου εργαστηρίου.

16 Στοιχεία συνεδρίας R

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Athens
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] deSolve_1.42  ggplot2_4.0.3
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.51         rlang_1.2.0       
##  [5] xfun_0.57          otel_0.2.0         S7_0.2.2           jsonlite_2.0.0    
##  [9] labeling_0.4.3     glue_1.8.1         htmltools_0.5.9    sass_0.4.10       
## [13] scales_1.4.0       rmarkdown_2.31     grid_4.6.1         evaluate_1.0.5    
## [17] jquerylib_0.1.4    fastmap_1.2.0      yaml_2.3.12        lifecycle_1.0.5   
## [21] compiler_4.6.1     RColorBrewer_1.1-3 farver_2.1.2       digest_0.6.39     
## [25] R6_2.6.1           bslib_0.10.0       tools_4.6.1        withr_3.0.2       
## [29] gtable_0.3.6       cachem_1.1.0