# Crecimiento ilimitado de la población

Un modelo elemental del crecimiento de una población se basa en la hipótesis de que:

## <font color =red>La velocidad de crecimiento de la población es proporcional al tamaño de la población</font>

1. Observe que la razón de cambio de una población sólo depende del tamaño de ésta. 

- En particular, las limitaciones del espacio o recursos no tienen efecto. 

- Esta hipótesis es razonable para pequeñas poblaciones en grandes entornos, por ejemplo, los primeros brotes de moho en una pieza de pan o los primeros colonizadores de Estados Unidos.

Como la hipótesis es tan simple, esperamos que el modelo también lo sea, Las cantidades implicadas son

- $t$ = tiempo (variable inependiente).
- $P$ = población (variable dependiente)
- $k$ = constante de proporcionalidad (parámetro) entre la tasa de crecimiento de la población y el tamaño de ésta.

** El parámetro $k$ suele llamarse "coeficiente de velocidad de crecimiento".** 

Las unidades para esas cantidades dependen de la aplicación.

- Si estamos modelando el crecimiento de moho en el pan, entonces $t$ podría medirse en días y $P(t)$ sería el área cubierta por el moho o bien el peso del moho. 

- Si estamos hablando de la población europea en Estados Unidos, entonces $t$ probablemente se medirá en años y $P(t)$ en millones de personas. En este caso haríamos corresponder $t=0$ a cualquier tiempo que quisiéramos. El año 1790 (el año del primer censo) es una opción conveniente.


## Expresemos ahora nuestra hipótesis usando esta notación.

La tasa de crecimiento de la población $P$ es la derivada $\frac{dP}{dt}$. Puesto que ésta es proporcional a la población, se expresa como el producto, $kP$, de la población $P$ y la constante $k$ de proporcionalidad. 

Por consiguiente, nuestra hipótesis se expresa por la ecuación diferencial

$$ \frac{dP}{dt} = kP $$

En otras palabras, la razón de cambio de $P$ es proporcional a $P$.
<font color=blue>** Éste es nuestro primer ejemplo de una ecuación diferencial.** </font>

Asociada con ella hay varios adjetivos que describen su tipo. En particula, se trata de una ecuación de **primer orden** porque contiene sólo primeras derivadas de la variable dependiente, y es una **ecuación differencial ordinaria** porque no contiene derivadas parciales. 

Hemos escrito esta ecuación diferencial usando la notación de Leibniz, es decir, $\frac{dP}{dt}$, y es la que tenderemos a usar. Sin embargo, hay muchas otras maneras de expresar la misma ecuación diferencial.

En particula, también podríamos escribirla como $P'=kP$ o como $\dot{P} = kP$. la notación "punto" suele utilizarse cuando la variable independiente es el tiempo $t$.


# ¿Qué predice el modelo?

Más importante que los adjetivos o cómo se escribe la ecuación es preguntar qué nos dice acerca de la situación que se está modelando. Como $\frac{dP}{dt}=kP$ para alguna constante $k$ y $\frac{dP}{dt}=0$ si $P=0$.

Entonces la función constante $P(t)=0$ es una solución de la ecuacioón diferencial. A este tipo especial se le denomina **solución de equilibrio ** porque es constante para siempre. En términos del modelo de población, corresponde a una especie que es no existente.

Si $P(t_0) \neq 0$ en algún tiempo $t_0$, entonces en el tiempo $t=t_0$

$$ \frac{dP}{dt} = kP(t_0) \neq 0 $$

En consecuencia, la población no es constante. Si $k>0$ y $P(t_0)>0$, tenemos

$$ \frac{dP}{dt}= kP >0$$

en el tiempo $t=t_0$ y la población está creciendo (como era de esperarse). 

Conforme $t$ crece, $P(t)$ se vuelve mayor, por lo que $ \frac{dP}{dt} $ aumenta. A su vez, $P(t)$ crece aún más rápidamente.

Es decir, la velocidad de crecimiento crece en relación directa con la población.

El valor de $P(t)$ en  $t=0 $ se llama una **condición inicial **.

<font color =red>Si comenzamos con una condición inicial diferente obtenemos una función $P(t)$ distinta,</font> 


```python

```

# Solucion analitica de la ecuación diferencial

Si, por otra parte, conocemos el valor exacto de $P_0$ de $P(0)$ y queremos predecir el valor de $P(10)$ o $P(100)$, entonces necesitamos información más precisa sobre la función $P(t)$. El par de ecuaciónes

## $$ \frac{dP}{dt}= kP \ \ \ \ P(0) = P_0$$

se llama **problema de valor inicial. ** Y una ** solución** al problema de valor inicial es una función $P(t)$ que satisface ambas ecuaciones. Es decir

## $$ \frac{dP}{dt}= kP \ \ y \ \ P(0) = P_0$$

En consecuencia, para solucionar esta ecuación diferencial debemos hallar una función $P(t)$ cuya derivada sea el producto de $k$ con $P(t)$.

# Separación de variables

## $$ \frac{dP}{P}= kdt$$

## $$\int \frac{dP}{P}= \int kdt$$

## $$ \ln P + c_1= k t +c_2$$
## $$ e^{\ln P} = e^{k t +c}$$
## $$ P(t) = e^{k t} e^{c}$$
## $$ P(t) = ce^{k t} $$

<font color=red>**Existe un número infinito de soluciones** para la ecuación diferencial, **una para cada valor de $c$.** </font>

Para determinar cuál de ésas es la correcta para la situación considerada, usamos la condición inicial dada.

## $$P_0= P(0) = c \cdot e^{k\cdot 0} = c \cdot e^0 = c \cdot 1 = c $$

En consecuencia, debemos escoger $c= P_0$, por lo que una solución del problema del valor inicial es 

## $$ P(t) = P_0e^{k t} $$

La función $P(t)$ se llama solución al problema del valor inicial así como **solución particular** de la ecuación diferencial.

El conjunto de funciones $P(t) = ce^{kt}$ se llama ** solucion general** de la ecuación general, porque podemos usarla para encontrar la respuesta particular correspondiente a cualquier problema de valor inicial.

# Ejemplo: Método de Euler

$$ \frac{dP}{dt}= 2P \ \ \ \ P(0) = 10 $$

$$ \frac{P_{k+1}-P_{k}}{h}= 2P_k $$


$$ P_{k+1}=P_{k} + h2P_k $$

Es evidente que la ordenada $y_1$ calculada de esta manera no es igual a $F(x1)$ pues existe un pequeño error. No obstante, el valor $y_1$ sive para aproximar $F'(x)$ en el punto $P=(x_1,y_1)$ y repeteri el procedimiento anterior a fin de generar la sucesión de aproximaciones siguiente


$$ y_1 = y_0 + h f(x_0,y_0)$$
$$ y_2 = y_1 + h f(x_1,y_1)$$
$$\vdots $$
$$ y_{i+1} = y_i + h f(x_i,y_i)$$
$$\vdots $$
$$ y_{n} = y_{n-1} + h f(x_{n-1},y_{n-1})$$


```python

```

# Ejemplo:  Solución analitica de la ecuación diferencial

$$ \frac{dP}{dt}= 2P \ \ \ \ P(0) = 10, \ \ \ k=2$$


## Separación de variables

$$ \frac{dP}{P}= 2dt$$

$$\int \frac{dP}{P}= \int 2dt$$

$$ \ln P + c_1= 2 t +c_2$$
$$ e^{\ln P} = e^{2 t +c}$$
$$ P(t) = e^{2 t} e^{c}$$
$$ P(t) = ce^{2 t} $$



$$P_0= P(0) = c \cdot e^{2\cdot 0} = c \cdot e^0 = c \cdot 1 = c $$

<font color=red>Solución: </font>
$$ P(t)=  10e^{2t}$$


- Ecuación diferencial
$$ \frac{dP}{dt}= 2P \ \ \ \ P(0) = 10 $$
- Método de Euler
$$ P_{k+1}=P_{k} + h2P_k $$
- Solución particular de la ecuación diferencial.
$$ P(t) = P_0e^{2 t} $$

# Ejemplo


```python
%matplotlib inline 
import matplotlib.pylab as plt
import numpy as np
import matplotlib as mpl
label_size = 14
mpl.rcParams['xtick.labelsize'] = label_size 
mpl.rcParams['ytick.labelsize'] = label_size
```


```python

```


```python
P = lambda t: 10*np.exp(2*t) 
    
def dP (Pk,h):
    Pk1= Pk + h*2*Pk
    return Pk1

def dP_dt (P,t):
    dP= 2*P
    return dP

tfin =3
h=.1

#Solucion exacta
t = np.arange(0,tfin+.1,.1)
plt.plot(t,P(t),"b")

#Utilizando metodo de euler    
Pk=np.arange(0,tfin+h,h)
Pk[0]=10

for k in np.arange(0,Pk.size-1) :
    Pk[k+1] = dP(Pk[k],h)

tk=np.arange(0,tfin+h,h)
plt.plot(tk,Pk,"g")    
    
#Utilizando odeint
t = np.arange(0,3,.1);
x = odeint(dP_dt,10,t)
plt.plot(t,x,'r');
```


```python

```

# La población de Estados Unidos

Consideremos las cifras de los censos de Estados Unidos desde 1790.

- Medimos el tiempo en años y la población $P(t)$ en millones de personas.

- Hacemos que $t=0$ sea el año $1790$, por lo que la condición inicial es $P(0) = 3.9$

## $$ \frac{dP}{dt}= kP \ \ \ \ P(0) = 3.9$$

## tiene $ P(t) = 3.9 e^{k t} $ como solución.

Pero no podemos usar este modelo para hacer predicciones porque **no conocemos el valor de $k$**.

Sin embargo, sabemos que la población en el año $1800$ era de $5.3$ millones y podemos usar este valor para determinar $k$. Si hacemos

## $$5.3 = P(10) = 3.9 e^{k\cdot 10 }$$

tenemos entonces

## $$  e^{k\cdot 10}  = \frac{5.3}{3.9} $$

## $$  10 k   = \ln \frac{5.3}{3.9} $$

## $$ k   \approx 0.03067. $$

## $$ P(t) = 3.9 e^{ 0.03067 \cdot t} $$


```python
año = np.arange(1790,2000,10)
t = np.arange(0,210,10)
real = np.array([3.9,5.3,7.2,9.6,12,17,23,31,38,50,62,75,91,105,122,131,151,179,203,226,249])

Pt = 3.9*np.exp(0.03067*t)

plt.plot(t,real)
plt.plot(t,Pt)

_=plt.xticks(t, año, rotation='vertical')
plt.margins(0.05)
plt.xlabel("t")
plt.ylabel("P(t)")
```



Nuestro modelo es bastante bueno siempre que la población sea relativamente pequeña.

Sin embargo, con el paso del tiempo el modelo predice que la población continuará creciendo sin límite, y obviamente esto no sucede en el mundo real.

En consecuencia, si queremos un modelo que sea exacto sobre una escala grande de tiempo, debemos tomar en cuenta el hecho de que la poblaciones existen en una cantidad finita de espacio y con recursos limitados. 


```python

```

# Ejercicio: La población de Estados Unidos



## $$ \frac{dP}{dt}= kP \ \ \ \ P(0) = 3.9$$

## $$ k   \approx 0.03067. $$


```python
def df (P,k,h):
    
    return P + h*k*P

h=.5
Pk1=np.arange(0,200+h,h)
t = np.arange(0,Pk1.size-1)
k=0.03067
Pk1[0]= 3.9

for i in  t:
    Pk1[i+1] = df(Pk1[i],k,h)

x= int (1/h)


t = np.arange(0,200+h,h)
plt.plot(t,Pk1,"r")



año = np.arange(1790,2000,10)
t = np.arange(0,210,10)
real = np.array([3.9,5.3,7.2,9.6,12,17,23,31,38,50,62,75,91,105,122,131,151,179,203,226,249])

Pt = 3.9*np.exp(0.03067*t)

plt.plot(t,real,"g")
plt.plot(t,Pt,"b")


_=plt.xticks(t, año, rotation='vertical')
plt.margins(0.05)
plt.xlabel("t")
plt.ylabel("P(t)")
```


```python

```


```python

```

# Modelo logístico de la población

Para ajustar el modelo de crecimiento exponencial de la población que tome en cuenta un entorno y recursos limitados, agregamos la hipótesis:

## Si la población es pequeña, la razón de crecimiento de la población es proporcional a su tamaño. 

## Si la población es demasiado grande para ser soportada por su entorno y recursos. la  población disminuirá. Es decir, la razón de crecimiento es negativa.

Para este modelo usamos de nuevo

- $t$ = tiempo (variable inependiente).
- $P$ = población (variable dependiente)
- $k$ = constante de proporcionalidad (parámetro) entre la tasa de crecimiento de la población y el tamaño de ésta.

Sin embargo, nuestra hipótesis acerca de recursos limitados introduce otra cantidad, el tamaño de la población que corresponde a ser **"demasiado grande"**. Esta cantidad es un segundo parámetro, denotado por $N$, que llamamos la **"capacidad de soporte"** del entorno. 

En términos de la capacidad de soporte, estamos suponiendo que $P(t)$ crece si $P(t)<N$. No obstante, si $P(t)>N$, suponemos que $P(t)$ está decreciendo.

Usando esta notación, podemos reescribir nuestras hipótesis como:

- $\frac{dP}{dt} \approx kP$ si $P$ es pequeña (primera hipótesis).

- Si $P>N,\  \frac{dP}{dt} <0 $ (segunda hipótesis).

Queremos también que el modelo sea **"algebraicamente simple"** o por lo menos tan simple como sea posible, por lo que tratamos de modificar el modelo exponencial lo menos posible. 

Por ejemplo, podríamos intentar una expresión de la forma

$$\frac{dP}{dt} = k (algo) P $$

Queremos que el factor **"algo"** sea cercano a $1$ si $P$ es pequeña, pero si $P>N$, queremos que **"algo"** sea negativo. La expresión más simple que tienen estas propiedades es la función

$$ (algo) = \left(1 - \frac{P}{N} \right)$$

Note que esta expresión es igual a $1$ si $P=0$ y es negativa si $P>N$. Nuestro modelo es entonces

$$\frac{dP}{dt} = k \left(1 - \frac{P}{N} \right) P $$

Éste se llama el **modelo logístico de la población ** con velocidad de crecimiento $k$ y capacidad $N$ de soporte.

Se trata de otra ecuación diferencial de primer orden. Se dice que esta ecuación es **no lineal ** por que su lado derecho no es una función lineal de $P$ como lo era en el modelo de crecimiento exponencial.


```python
def df (P,k,h,N):
    return P + h*k*(1-P/N)*P

h=.01
N=203

Pk1=np.arange(0,200+h,h)
t = np.arange(0,Pk1.size-1)

k=0.03067

Pk1[0]= 3.9

for i in  t:
    Pk1[i+1] = df(Pk1[i],k,h,N)

x= int (100/h)
print (Pk1[x],t[x])

t = np.arange(0,200+h,h)
plt.plot(t,Pk1,"r",label="Modelo")



año = np.arange(1790,2000,10)
t = np.arange(0,210,10)
real = np.array([3.9,5.3,7.2,9.6,12,17,23,31,38,50,62,75,91,105,122,131,151,179,203,226,249])

Pt = 3.9*np.exp(0.03067*t)

plt.plot(t,real,"g",label="Real")

#plt.plot(t,Pt,"b")


_=plt.xticks(t, año, rotation='vertical')
plt.margins(0.05)
plt.xlabel("t")
plt.ylabel("P(t)")

plt.legend(loc = 'best', prop={'size': 14})    
```

# Ecuaciones de Volterra: modelo presa depredador


Se trata de un sistema de dos ecuaciones diferenciales de primer orden, acopladas, autónomas y no lineales:

$$ \displaystyle \frac{dx}{dt} = \alpha x - \beta x y $$

$$ \displaystyle \frac{dy}{dt} = -\gamma y + \delta y x $$

donde x es el número de presas (cebras en nuestro caso) e y es el número de depredadores (leones). Los parámetros son constantes positivas que representan:

- $\alpha$: tasa de crecimiento de las presas.
- $\beta$: éxito en la caza del depredador.
- $\gamma$: tasa de decrecimiento de los depredadores.
- $\delta$: éxito en la caza y cuánto alimenta cazar una presa al depredador.

Aunque este modelo puede parecer relativamente simple, ha sido la base de algunos interesantes estudios ecológicos. En particular, Volterra y D'Ancona usaron con éxito el modelo para explicar el incremento en la población de tiburones en el mar Mediterráneo durante la Primera Guerra Mundial, cuando la pesca de las especies **"presa"** decreció. El modelo puede también usarse como base para el estudio de los efectos de los pesticidas en la población de insectos depredadores e insectos presas.


```python
from scipy.integrate import odeint

def df_dt(x, t, a, b, c, d):
    """Función del sistema en forma canónica"""
    dx = a * x[0] - b * x[0] * x[1]
    dy = - c * x[1] + d * x[0] * x[1]
    return np.array([dx, dy])

# Parámetros
a = 0.1
b = 0.02
c = 0.3
d = 0.01

# Condiciones iniciales
x0 = 40   # Presas
y0 = 9    # Depredadores
conds_iniciales = np.array([x0, y0])

# Condiciones para integración
tf = 200
N = 800
t = np.linspace(0, tf, N)

solucion = odeint(df_dt, conds_iniciales, t, args=(a, b, c, d))
#plt.plot(t, solucion[:, 0], label='presa')
#plt.plot(t, solucion[:, 1], label='depredador')
plt.plot(solucion[:, 0], solucion[:, 1])
plt.xlabel('Presas', fontsize = 14) 
plt.ylabel('Depredadores', fontsize = 14) 
```

## Ecuaciones de Volterra: modelo presa depredador

$$
\begin{align}
\dot{x}_1 &=-x_1 +x_1x_2,\\
\dot{x}_2 &=x_2 -x_1x_2
\end{align}
$$


```python
def Voltera (x,t):
    
    return [-x[0] + x[0]*x[1],x[1]-x[0]*x[1]]

t = np.arange(0,200,.1)
x = odeint(Voltera,[1.2,1.2],t)
plt.plot(x[:,0],x[:,1])
```


```python
t = np.arange(0,3,.1)

x1 = np.arange(-2,2,.2)
x2=x1

N= len(x1)

for xr in x1:
    for xc in x2:
        x0 = [xr,xc]
        
        x = odeint(Voltera,x0,t)
        plt.plot(x[:,0],x[:,1])
        plt.hold(True)
        
plt.axis([-5,5,-5,5])
```

# Teorema 1

La ecuación diferencial lineal de orden $n$

$$x^{(n)} + a_1x^{(n-1)}+ a_2x^{(n-2)} + \ldots + a_{n-1}x' + a_nx =f(t) $$

puede expresarse como un sistema de $n$ ecuaciónes lineales de primer orden.

*Demostración* Defina $x_1=x ,x_2=x',x_3=x'', \ldots, x_n=x^{n-1}$ Entonces tenemos


$$ \dot{x}_1 = x_2$$
$$ \dot{x}_2 = x_3$$
$$ \vdots$$
$$ \dot{x}_{n-1} = x_n$$
$$ \dot{x}_n = -a_nx_1 -a_{n-1}x_2 - \ldots - a_1x_n +f$$

# Ejemplo

$$t^3x''' + 4t^2x'' - 8tx' + 8x =0  $$

$$x(0)=3,x'(0)=-6,x''(0)=14 $$.


Definiendo $x_1=x,x_2=x',x_3=x''$, obtenemos el sistema

$$\dot{x}_1 = x_2 $$
$$\dot{x}_2 = x_3 $$
$$\dot{x}_3 =  -\frac{8}{t^3}x_1 +\frac{8}{t^2}x_2 - \frac{4}{t}x_3  $$


```python
def df_dt2(x, t):
    """Función del sistema en forma canónica"""
    dx1 = x[1];
    dx2 = x[2]
    dx3 = - 8/t**3 *x[0] + 8/t**2 *x[1] - 4/t *x[2]
    return np.array([dx1, dx2 , dx3])

# Condiciones iniciales
x0 = 3   
x1 = -6
x2 = 14
conds_iniciales = np.array([x0,x1,x2])

# Condiciones para integración
tf = 10
N = 800
t = np.linspace(0.1, tf, N)

solucion = odeint(df_dt2, conds_iniciales, t)
plt.plot(t, solucion[:, 0])
plt.plot(t, solucion[:, 1])
plt.plot(t, solucion[:, 2])
#plt.plot(solucion[:, 0], solucion[:, 1])
```



# Un Modelo de ahorro

Supongamos que depositamos $ \$5000 $ en una cuenta de ahorro con interés incrementado a una tasa de $ 5\% $ compuesto en forma continua. Si $ A(t)$ denota la cantidad de dinero en la cuenta en el tiempo $ t$, 


entonces la ecuación diferencial para $A$ es

$$ \frac{dA}{dt} = 0.05A $$

Como vimos en la sección previa, la solución general para esta ecuación es la función exponencial

$$A(t) = ce^{0.05t} $$

donde $c = A(0).$ Entonces

$$ A(t) = 5000e^{0.05t}$$
es nuestra solución particular.

Suponiendo que la tasa de interés nunca cambian, después de 10 años tendremos 

$$ A(10) = 5000 e^{0.5} \approx  8244  $$

Pesos en esta cuenta.

# Simula la solución analitica de la ecuación diferencial


```python
def fuc(t):
    return 5000*np.exp(0.05*t)


t = np.arange(0,11,1)

A=fuc(t)

plt.plot(t,A)
```


```python
fuc(10)
```




    8243.6063535006415



# Simula usando el Método de Euler

$$ \frac{dA}{dt} = 0.05A $$

$$ A(0) = 5000$$

$$ \frac{A_{k+1} - A_k}{T} = 0.05A $$
$$ A_{k+1} =A_k + T (0.05A) $$


```python
T=1;
A = np.arange(0,11,T)
K = np.arange(0,A.size-1,1)
A[0] = 5000

for k in K:
    A[k+1] = A[k] + T* (0.05*A[k])

t = np.arange(0,11,T)
plt.plot(t,A)
```


```python
A[int (10/T) ]
```




    8138



# Simula usando la funcion odeint de scipy

# $ \frac{dx}{dt} $

# $ \frac{dA}{dt} $



```python
from scipy.integrate import odeint


def dA_dt(A,t):
    dA =(0.05*A)
    return dA
    
A0 = 5000   
T=1
t = np.arange(0,11,T)
solucion = odeint(dA_dt,A0,t) 

plt.plot(t,solucion)
print (solucion[int (10/T)])
```

# Compara los tres resultados para A(10) = ?


```python

```

20a75ec735333fd99cb410f4c83cfc488b093565
ipynb
Jupyter Notebook
Intro/Clase8/.ipynb_checkpoints/EQDifferencial-checkpoint.ipynb
RoboticaIndustrial/Calando
a539a039404c5700e2d7d34e5b93c2f8be884f26
MIT
Intro/Clase8/.ipynb_checkpoints/EQDifferencial-checkpoint.ipynb
RoboticaIndustrial/Calando
a539a039404c5700e2d7d34e5b93c2f8be884f26
MIT
Intro/Clase8/.ipynb_checkpoints/EQDifferencial-checkpoint.ipynb
RoboticaIndustrial/Calando
a539a039404c5700e2d7d34e5b93c2f8be884f26
MIT
Qwen/Qwen-72B
1. YES
2. YES