Cresimento Logístico com Python

Crescimento Logístico

O crescimento logístico leva em conta que a população tem um valor máximo \(M\). Quando a população se aproxima da capacidade máxima, os recursos tornam-se menos abundantes e a taxa de crescimento começa a diminuir. Se \(y=y(t)\) representa o número de indivíduos no instante \(t\), uma relação simples que exibe esse comportamento é quando \[\frac{dy}{dt}=ky(M-y)\] Vamos usar o python para nos ajudar a resolver o seguinte exemplo:

Exemplo

Biólogos colocaram em um lago 400 peixes e estimaram a capacidade de suporte como 10.000. O número de peixes triplicou no primeiro ano. Encontre uma expressão para o tamanho da população de peixes depois de \(t\) anos.

Solução: Vamos trabalhar com as populações em unidades de mil, assim \(y(0)=0.4\) e \(y(1)=1.2\). Primeiramente, vamos resolver a EDO usando o método das equações separáveis. \[y'(t)=ky(t)(10-y(t))\Rightarrow \frac{y'(t)}{y(10-y(t))}=k\Rightarrow \int\frac{y'(t)}{y(10-y(t))}dt =kt+c\Rightarrow \int\frac{dy}{y(10-y)}=kt+c.\]

Com isso o problema se resume em determinar a integral \[\int\frac{dy}{y(10-y)}.\]

Note que \[\int\frac{dy}{y(10-y)}=\frac{1}{10}\int\frac{1}{y}+\frac{1}{10-y}dy=\frac{1}{10}(\log|y|-\log|10-y|)+c_1\]

Com isso, \[\log\left|\frac{y}{10-y}\right|=kt+c_2\Rightarrow \frac{y}{10-y}=Ce^{kt}\]

Agora, vamos usar a biblioteca sympy do python para substituir o dado inicial \(y(0)=0.4\) a fim de obter \(C\).

import sympy as sp

#variáveis
y,k,C=sp.symbols('y k C',real=True)
t=sp.symbols('t',positive=True)

#dados do problema
y0=sp.Rational(400,1000)
y1=3*y0

eq1=sp.Eq(y/(10-y),C*sp.exp(k*t))
sol=sp.solve(eq1.subs([(t,0),(y,y0)]),C)
c0=sol[0]
print('C=',c0)

C= 1/24

Neste caso obtemos que \(C=\frac{1}{24}\). Substituindo na equação obtida ficamos com

eq2=eq1.subs(C,c0)
eq2

\(\displaystyle \frac{y}{10 - y} = \frac{e^{k t}}{24}\)

Agora vamos encontrar \(k\) usando o segundo dado \(y(1)=1.2\).

sol2=sp.solve(eq2.subs([(t,1),(y,y1)]),k)
k0=sol2[0]
print('k=',k0)

k= log(36/11)

Vamos substituir na equação e simplificar:

eq3=eq2.subs(k,k0)
sp.simplify(eq3)

\(\displaystyle \frac{\left(\frac{36}{11}\right)^{t}}{24} = - \frac{y}{y - 10}\)

Agora vamos isolar \(y\), encontrando finalmente a função buscada.

sol3=sp.solve(eq3,y)
Y=sol3[0]
sp.Eq(y,Y)

\(\displaystyle y = \frac{10 \left(\frac{36}{11}\right)^{t}}{\left(\frac{36}{11}\right)^{t} + 24}\)

Vamos agora esboçar o gráfico desta função no intervalo \([0,10]\).

sp.plot(Y,(t,0,10))

Note que se dividirmos o numerador e o denominador por \(\left(\frac{36}{11}\right)^t\), podemos reescrever a função da seguinte forma:

\[y(t)=\frac{10 }{1 + 24\left(\frac{36}{11}\right)^{-t}},\ t\geq 0.\] Por fim, vamos fazer previsões:

Exemplo

Determine em quanto tempo a população chegará a 5 mil.

Queremos encontrar o tempo \(T\) para o qual \(y(T)=5\), isto é, resolver a equação:

\[5=\frac{10 }{1 + 24\left(\frac{36}{11}\right)^{-t}}.\]

sol3=sp.solve(sp.Eq(Y,5),t)
T=sol3[0]
print('Instante em que a população chega a 5000: ', T)
print('Valor aproximado: ', T.evalf(), 'anos' )
print('Convertendo os decimais em meses', (T.evalf()-2)*12 )

Instante em que a população chega a 5000:  -log(24**(1/log(11/36)))
Valor aproximado:  2.68049122365053 anos
Convertendo os decimais em meses 8.16589468380641

Assim, temos que a população chega a 5.000 peixes em aproximadamente 2 anos e 8 meses.