Calcoli scientifici con Julia/Crescita popolazione

Modello di crescita

Avendo supposto che il numero di individui di una popolazione sia una funzione continua del tempo $N(t)$ che ammette derivata continua, si ha che l'incremento della popolazione al variare del tempo può essere rappresentato dalla derivata di $N(t)$, che in un modello elementare si può supporre direttamente proporzionale al numero di individui della popolazione stessa.

Si ha pertanto la seguente equazione differenziale:

${d \over dt}N=rN(t),$

con $r$: parametro di crescita malthusiana (tasso massimo di crescita della popolazione).

Pertanto se $r$ è una costante la popolazione cresce in maniera esponenziale con pendenza dipendente da $r$.

Invece in un ambiente la cui disponibilità di risorse è limitata si può descrivere l'evoluzione della popolazione utilizzando un coefficiente $r$ che decresce all'aumentare della popolazione: il modello più semplice è $r(t):=a-bN(t)$ con $a$ e $b$ costanti. Sostituendo tale funzione nella precedente equazione differenziale si ottiene:

${\frac {dN}{dt}}=aN(t)-bN^{2}(t)$

che può essere posta nella forma:

${\frac {dN}{dt}}=aN\left(1-{\frac {N}{K}}\right)$

con $K={\frac {a}{b}}$ che è la cosiddetta popolazione massima sostenibile e $a$ uguale al parametro di crescita malthusiana. Questa è l'equazione logistica di Verhulst.

Separando le variabili si ottiene:

${\frac {1}{a}}\int \left({\frac {1}{N}}+{\frac {b}{a-bN}}\right)dN(t)=\int dt.$

Soluzione equazione differenziale

Risolvendo gli integrali, scegliendo come primitive quelle tali che $N(t_{0})=N_{0}$ e utilizzando le proprietà dei logaritmi si ottiene la soluzione:

$N(t)={\frac {k}{1+({\frac {k}{N_{0}}}-1)e^{-a(t-t_{0})}}}.$

Si nota che a causa del sovraffollamento la popolazione non cresce più in maniera esponenziale ma converge al valore asintotico $k$ indipendentemente da $N_{0}$.

La soluzione dell'equazione si può anche scrivere nelle forme:

$N(t)={kN_{0}e^{rt} \over k-N_{0}\cdot (1-e^{rt})}={kN_{0} \over N_{0}+(k-N_{0})\cdot e^{-rt}}.$

È immediato verificare che questa soluzione ha due asintoti orizzontali:

$\lim _{n\to \infty }{N(t)}=K;$

$\lim _{n\to -\infty }{N(t)}=0.$

Si ha un differente comportamento nel caso $N_{0}>k$ allora il secondo limite tenderebbe a $-\infty$, presentando anche un asintoto verticale, ma queste soluzioni non sono considerate nel modello di crescita (descrivono evidentemente una popolazione in rapida decrescita in quanto inizialmente in eccesso rispetto alle risorse presenti).

Esempio di calcolo con Julia

Supponiamo che una popolazione di conigli cresca al tasso percentuale del 30% mensile. Supponiamo che la popolazione iniziale sia composta da due conigli.
a. Quanti conigli ci saranno dopo 1 anno?
b. Dopo quanto tempo la popolazione raggiungerà i 1000 esemplari?

Innanzitutto occorre installare i pacchetti Plots e DifferentialEquations.

Da REPL digitare:

using Pkg
Pkg.add("Plots")

using Pkg
Pkg.add("DifferentialEquations")

Da Jupyter in due celle di codice digitare :

] add Plots

] add DifferentialEquations

e poi cliccare su Run. Poi eseguire il codice per risolvere l'equazione differenziale e tracciare il grafico della soluzione analitica e numerica della crescita malthusiana esponenziale (utilizzo il seguente codice rilasciato con licenza MIT https://computationalmindset.com/en/neural-networks/ordinary-differential-equation-solvers-in-julia.html):

using DifferentialEquations
using Plots

ode_fn(N,p,t) = 0.3N

an_sol(t) = 2*exp(0.3*t)

t_begin=0.0
t_end=20.0
tspan = (t_begin,t_end)
x_init=2.0

prob = ODEProblem(ode_fn, x_init, tspan)
num_sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

plot(num_sol.t, an_sol.(num_sol.t),
    linewidth=2, ls=:dash,
    title="Crescita malthusiana tramite DifferentialEquations",
    xaxis="t", yaxis="N",
    label="analitico",
    legend=true)
plot!(num_sol,
    linewidth=1,
    label="numerico")

Dopo 1 anno ci saranno 73 conigli, mentre si otterranno 600 conigli dopo 19 mesi:

num_sol(12)

73.19646891230852

using Roots
y(t)=2*exp(0.3*t)-600;
find_zero(y,18)

19.012608248854004

Nel caso in cui le risorse disponibili sono limitate la crescita non sarà più esponenziale ma raggiungerà ad esempio un valore limite di 1500 conigli tramite l'equazione logistica calcolata con Julia e visualizzata con grafico:

ode_fn(N,p,t) = 0.3N*(1-N/1500)

an_sol(t) = 1500/(1+(1500/2-1)*exp(-0.3*t))

t_begin=0.0
t_end=50.0
tspan = (t_begin,t_end)
x_init=2.0

prob = ODEProblem(ode_fn, x_init, tspan)
num_sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

plot(num_sol.t, an_sol.(num_sol.t),
    linewidth=2, ls=:dash,
    title="Equazione logistica tramite DifferentialEquations",
    xaxis="t", yaxis="N",
    label="analitico",
    legend=true)
plot!(num_sol,
    linewidth=1,
    label="numerico")