MateWiki - Contribuciones del usuario [es]

Difusión de un contaminante. Grupo 4

2014-05-19T22:38:49Z

Adrián:

{{ TrabajoED | Difusión de una sustancia contaminante. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina '''''difusión''''', proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

[[Archivo:Barra2.png |marco|centro|Configuración inicial de las concentraciones de contaminante en el tubo]]

== Modelización del problema ==

Consideramos un tubo de longitud <math>L</math> de un cierto material, con sección transversal <math>A</math> constante. La orientamos en la dirección del eje <math>x</math>, de <math>x = 0</math> a <math>x = L</math>. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de <math>x</math> y <math>t</math>. Designamos por <math>u(x, t)</math> la concentración de contaminante en la sección de la varilla que dista <math>x ≥ 0</math> del extremos <math>x = 0</math> cuando ha pasado un tiempo <math>t ≥ 0</math>.

Tomemos un trozo de tubo entre las secciones <math>x</math> y <math>x + ∆x</math>, que designaremos
por <math>[x, x + ∆x]</math>. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por <math>φ(x, t)</math> el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje <math>X</math>. Si <math>φ (x, t) > 0</math> pensamos que el flujo va hacia la derecha, si <math>φ (x, t) <0</math> va hacia la izquierda. El flujo de contaminante en <math>[x, x + ∆x]</math> está dado por <math>φ (x, t)A − φ (x + ∆x, t)A</math>.

Designamos por <math>f(x, t)</math> a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto <math>f(x)=0</math>.

La ''ley de Fick'' (análoga a la de Fourier para la temperatura) dice:
<math>F=-D \frac{∂u}{∂x}</math>

Fick demostró empíricamente que:
<math> φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}</math>

De aquí se deduce la anteriormente citada ley cuando suponemos que <math>∆x</math> tiende a <math>0</math>.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene <math>ux(x,t)>0</math>,
y <math>φ (x, t)</math> será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene <math>ux(x,t)<0</math>,
y <math>φ (x, t)</math> será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en <math>\cfrac{m^2}{s}</math>, que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ''ley de Fick'', la concentración en la sección transversal del tubo que dista x del extremo <math>x = 0</math> cuando ha pasado un tiempo <math>t</math>, <math>u(x, t)</math>, satisface:
<math>u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>x ∈(0,L),t>0</math>

es decir, la ecuación en derivadas parciales:
<math>u_t= \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>u_t-Du_{xx} = f(x,t)</math>

Como hemos señalado anteriormente <math>f(x,t)=0</math> y <math>D=1</math>:
<math>u_t-u_{xx}=0</math>

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

== Planteamiento del sistema ==
El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Supondremos D constante e igual a 1 y L con valor 5.

== Resolución del sistema ==
Método de diferencias finitas con ∆x=0.1 suponiendo que en el instante inicial se verifica:

:<math>
u(x,0)=\left\{\begin{matrix}\\0 & x≤3\\3 & x>3\end{matrix}\right.
</math>

=== Método del trapecio ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h/4;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

}}

[[Archivo:Trapeciog4.jpg]]
=== Método de Euler explícito ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler explícito
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=uu+dt*(-K*uu+F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

}}

[[Archivo:Explicitog4.jpg ]]

=== Método de Euler implícito ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler implícito
A=eye(N+1)+dt*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\(uu+dt*F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
}}

[[Archivo:Eulerimpg4.jpg]]

=== Método de Euler modificado ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler modificado
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
k1=-K*uu;
k2=-K*(uu+dt*k1);
uu=uu+dt/2*(k1+k2);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
}}
[[Archivo:Eulermodg4.jpg]]

===Comparativa de métodos===

En términos de precisión cabe destacar que es más adecuado el método de Euler implítico frente a Euler Modificado, o Runge-Kutta de orden 2.
El método que presenta peores resultados es Euler explítico dado que es muy inestable para sistemas.

Se muestra a continuación gráficas explicativas entre los tres métodos:

[[Archivo:Error_Euler_implícito_y_modificado_(ventana_grande).jpg |marco|centro|Error entre Euler implícito y Euler modificado]]
[[Archivo:Imageng4.jpg |marco|centro|Error entre Euler explícito y Euler implícito]]
[[Archivo:Error_Euler_explícito_y_modificado_(ventana_grande).jpg ||marco|centro|Error entre Euler explícito y Euler modificado]]

== Conservación de la masa total de contaminante ==

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación <math>u_t-u_{xx}=0</math> , basta con resolver la siguiente integral:

<math>\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0</math>

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración <math>u(x, t)</math> , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.
{{matlab|codigo=

A=U';
integral=trapz(A);
concentracion=U(:,26);
subplot(1,2,1)
plot(t,integral)
xlabel('Tiempo')
ylabel('Masa total de contaminante')
title('Evolución de la cantidad de contaminante en el sistema')
subplot(1,2,2)
plot(t,concentracion)
xlabel('Tiempo')
ylabel('Concentración en el punto medio del tubo')
title('Evolución de la concentración en el punto medio del tubo')

}}

[[Archivo:Integralg4.jpg]]
Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

== Resultados para tiempos grandes ==

Despreciamos <math>u_t</math>:

:<math>
\left\{\begin{matrix}\\-u_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Usamos una T=100 y hacemos los cálculos con el método del trapecio:

{{matlab|codigo=
clear all
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h/4;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
%Solución estacionaria
a=U(length(t),:);
b=U(1,:);
c=U((length(t)-1)/100+1,:);
d=U((length(t)-1)*2/100+1,:);
e=U((length(t)-1)*10/100+1,:);
hold on
plot(x,a,'b','linewidth',2);
plot(x,b,'-y','linewidth',1);
plot(x,c,'-g','linewidth',1);
plot(x,d,'-k','linewidth',1);
plot(x,e,'-r','linewidth',1);
legend('100','0','1','2','10');
xlabel('Posición')
ylabel('Concentracion')
hold off
}}

Y en la gráfica comprobamos las distintas distancias.

[[Archivo:Grafg4.jpg]]

== Cálculo del tiempo necesario para alcanzar el estado estacionario con un error del 5% ==

{{matlab|codigo=
clear all
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1;
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Estado estacionario de la concentración
est=U(length(t),:);
clear A F K L N T U dt h n t u0 uu x
% Repetimos el cálculo y lo paramos al alcanzar un error máximo del 5%
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
tiempo=(b-1)*dt
clear A F K L N T U a b dt h n t u0 uu x
% Cálculo del tiempo necesario para alcanzar el estado estacionario con un
% error del 5% para un intervalo espacial 10 veces más pequeño
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1/10;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(301,1);3*ones(200,1)];
% Discretización temporal
dt=(h*10)^2/2; % Mantenemos el mismo intervalo temporal
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1;
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Estado estacionario de la concentración
est=U(length(t),:);
clear A F K L N T U a b dt h n t u0 uu x
% Comprobamos si varía el tiempo al hacer h 10 veces más pequeño
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1/10;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(301,1);3*ones(200,1)];
% Discretización temporal
dt=(h*10)^2/2; % Mantenemos el mismo intervalo temporal
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
nuevo_tiempo=(b-1)*dt
}}
El tiempo que tarda en alcanzar la concentración el estado estacionario con un error del 5% es 8.68 segundos con h=0.1. Si cambiamos la longitud de paso a 0.01, el nuevo tiempo es practicamente invariable, 8.645 segundos.

== Caso de colocar un limpiador en un extremo ==

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right.
</math>

=== Método de diferencias finitas ===

:<math>
\left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right.
</math>

=== Resolución numérica ===

Al cambiar las condiciones de contorno, el valor estacionario de la concentración <math> u(x,t)</math> es 0 es todos los puntos del tubo, lo que se traduce en que el contaminante desaparece por completo cuando el tiempo es muy grande.

Será preciso calcular el error relativo en cada punto del tubo, obtener el valor absoluto y escoger el máximo valor.

<math>error \ relativo=\cfrac{valor \ real-valor \ estacionario}{valor \ estacionario}=5 </math>%

El tiempo que se tarda en alcanzar el valor del 5% es dependiente del valor que consideramos como estacionario, que en nuestro caso será la concentración correspondiente a T=100 segundos desde que se inicia el ensayo (elegimos este tiempo ya que es suficientemente pequeño el valor de la concentración <math> u(x,t)</math>. El resultado que se obtiene de tiempo para alcanzar dicho error del 5% es 99.49 segundos, quedando reflejado a continuación:

{{matlab|codigo=
% Difusión de una sustancia contaminante con limpiador en el extremo de la
% izquierda
% u_t-u_xx=0, x en (0,L)
% u(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
% Valor estacionario de la concentración en el tubo
est=U(length(t),:);
clear A F K L N T U dt h n t u0 uu x xx tt
% Repetimos el cálculo y lo paramos al alcanzar un error máximo del 5%
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
tiempo=(b-1)*dt

}}

== Método de Fourier ==

Para la resolución del sistema mediante el método de Fourier:

<math>
\left\{\begin{matrix}\\u_t-u_{xx}=0 \\u_x (0,t)=0 \\u_x (L,t)=0 \\u(x,t)=u_0\end{matrix}\right.
</math>

Buscaremos soluciones del tipo <math>u(x,t)=X(x)T(t)</math>:

<math>XT^{'}-X^{''} T=0</math>:
<math>\cfrac{T'}{T}=\cfrac{X^{''}}{X}=-μ</math>

El problema de autovalores asociado, será:

<math>
\left\{\begin{matrix}\\X^{''}-μX=0 \\X'(0)=0 \\X^{'} (L)=0\end{matrix}\right.
</math>:

<math>μ_k=(\cfrac{kπ}{L})^2</math>:

<math>φ_k=cos⁡(\cfrac{kπ}{L} x)</math>

Por otro lado, la ecuación diferencial que debe cumplir <math>T(t)</math> es:

<math> T^{'}_k (t)-μ_k T_k (t)=0</math>

Resolviéndola obtenemos:

<math> T_k (t)=C_k e^{-μ_k t}</math>

Entonces, la solución aproximada mediante el método de Fourier es la siguiente:

<math>u(x,t)= \sum_{k=0}^K C_k e^{-μt} φ_k </math>

A continuación se muestra el código de matlab que resuelve el problema mediante Fourier usando 1, 3, 5, 10 y 20 términos de la serie.

{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% Método de Fourier
% ______________________________________
%
clear all
% Datos del problema
L=5; % Longitud del tubo
a=0; b=5; % Extremos
T=5; % Tiempo final
% Numero de autofunciones de la serie = M
M=20;
% Mallado
% En x --------> Paso espacial = h
h=(b-a)/50;
x=a:h:b;
% En t --------> Paso temporal = j
j=h/4;
t=0:j:T;
[Mx,Mt]=meshgrid(x,t);
f=[zeros(1,31),3*ones(1,20)]; % Valor inicial
U=0;
for k=0:M
p=cos(k*pi*x/L); % Autofunciones
Ck=trapz(x,f.*p)/trapz(x,p.^2);
U=U+Ck*(exp(-(k^2*pi^2/L^2)*Mt).*cos(k*pi*Mx/L));
end
surf(Mx,Mt,U)
xlabel('Espacio');
ylabel('Tiempo');
}}

Para obtener los diferentes resultados basta con introducir en M el número de términos de la serie de Fourier que vayamos a usar. Para M=20 se obtiene un resultado similar a los obtenidos anteriormente con otros métodos

[[Archivo:Fourier020.jpg ]]

y si lo calculamos para M=50 se aprecia que la gráfica se suaviza y se asemeja aun más.

[[Archivo:Fourier050.jpg ]]

Las aproximaciones en t=0.5 para 1, 3, 5, 10 y 20 términos de Fourier serán

[[Archivo:comparacionfourier.jpg ]]

== Cambio en las condiciones de frontera ==
El problema pasa a ser
:<math>
\left\{\begin{matrix}\\u_t-u_{xx}=0\\u(5,t)=10sent\\u_x(0,t)=3\\u(x,0)=u_0\end{matrix}\right.
</math>
Este sistema tiene una condición de contorno tipo Dirichlet en el extremo de la varilla <math>x=L</math>.
Esto quiere decir que este extremo está en contacto con una fuente de contaminante siendo la concentración en el instante t en dicho extremo <math>u=10sent</math>

La otra condición es de tipo Neumann en cero <math>-Du_x(0,t)=φ(t)</math>

En este caso <math>φ(t)=3>0</math>, además <math>D=1>0</math>, entonces debe verificarse <math>u_x(0,t)<0</math>. Cogiendo <math>∆x>0</math> y pequeño obtenemos

<math>u(0+∆x,t)-u(0,t)<0</math>

<math>u(0+∆x,t)<u(0,t)</math>

lo que indica que el flujo de contaminante irá hacia la derecha,interpretado de forma que por el extremo x=0 está entrando una cantidad de contaminante en el instante t igual a φ=3.

{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=3
% u(L,t)=f(t)
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=10;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h; % Nodos interiores
% Aproximación de -u_xx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,2)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N,1);
F(1)=-3*2/h;
% Vector u0
u0=[zeros(31,1);3*ones(19,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Función f(t)
f=10*sin(t);
% Método de Euler explícito
U(1,:)=[u0' f(1)];
uu=u0;
for n=1:length(t)-1
F(N)=f(n)/h^2;
uu=uu+dt*(-K*uu+F);
U(n+1,:)=[uu' f(n+1)];
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
}}
[[Archivo:ultimog4.jpg]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Difusión de un contaminante. Grupo 4

2014-05-19T22:30:12Z

Adrián:

{{ TrabajoED | Difusión de una sustancia contaminante. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina '''''difusión''''', proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

[[Archivo:Barra2.png |marco|centro|Configuración inicial de las concentraciones de contaminante en el tubo]]

== Modelización del problema ==

Consideramos un tubo de longitud <math>L</math> de un cierto material, con sección transversal <math>A</math> constante. La orientamos en la dirección del eje <math>x</math>, de <math>x = 0</math> a <math>x = L</math>. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de <math>x</math> y <math>t</math>. Designamos por <math>u(x, t)</math> la concentración de contaminante en la sección de la varilla que dista <math>x ≥ 0</math> del extremos <math>x = 0</math> cuando ha pasado un tiempo <math>t ≥ 0</math>.

Tomemos un trozo de tubo entre las secciones <math>x</math> y <math>x + ∆x</math>, que designaremos
por <math>[x, x + ∆x]</math>. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por <math>φ(x, t)</math> el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje <math>X</math>. Si <math>φ (x, t) > 0</math> pensamos que el flujo va hacia la derecha, si <math>φ (x, t) <0</math> va hacia la izquierda. El flujo de contaminante en <math>[x, x + ∆x]</math> está dado por <math>φ (x, t)A − φ (x + ∆x, t)A</math>.

Designamos por <math>f(x, t)</math> a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto <math>f(x)=0</math>.

La ''ley de Fick'' (análoga a la de Fourier para la temperatura) dice:
<math>F=-D \frac{∂u}{∂x}</math>

Fick demostró empíricamente que:
<math> φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}</math>

De aquí se deduce la anteriormente citada ley cuando suponemos que <math>∆x</math> tiende a <math>0</math>.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene <math>ux(x,t)>0</math>,
y <math>φ (x, t)</math> será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene <math>ux(x,t)<0</math>,
y <math>φ (x, t)</math> será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en <math>\cfrac{m^2}{s}</math>, que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ''ley de Fick'', la concentración en la sección transversal del tubo que dista x del extremo <math>x = 0</math> cuando ha pasado un tiempo <math>t</math>, <math>u(x, t)</math>, satisface:
<math>u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>x ∈(0,L),t>0</math>

es decir, la ecuación en derivadas parciales:
<math>u_t= \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>u_t-Du_{xx} = f(x,t)</math>

Como hemos señalado anteriormente <math>f(x,t)=0</math> y <math>D=1</math>:
<math>u_t-u_{xx}=0</math>

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

== Planteamiento del sistema ==
El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Supondremos D constante e igual a 1 y L con valor 5.

== Resolución del sistema ==
Método de diferencias finitas con ∆x=0.1 suponiendo que en el instante inicial se verifica:

:<math>
u(x,0)=\left\{\begin{matrix}\\0 & x≤3\\3 & x>3\end{matrix}\right.
</math>

=== Método del trapecio ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h/4;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

}}

[[Archivo:Trapeciog4.jpg]]
=== Método de Euler explícito ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler explícito
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=uu+dt*(-K*uu+F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)

}}

[[Archivo:Explicitog4.jpg ]]

=== Método de Euler implícito ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler implícito
A=eye(N+1)+dt*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\(uu+dt*F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
}}

[[Archivo:Eulerimpg4.jpg]]

=== Método de Euler modificado ===
{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=5;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método de Euler modificado
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
k1=-K*uu;
k2=-K*(uu+dt*k1);
uu=uu+dt/2*(k1+k2);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
}}
[[Archivo:Eulermodg4.jpg]]

====Comparativa de métodos====

En términos de precisión cabe destacar que es más adecuado el método de Euler implítico frente a Euler Modificado, o Runge-Kutta de orden 2.
El método que presenta peores resultados es Euler explítico dado que es muy inestable para sistemas.

Se muestra a continuación gráficas explicativas entre los tres métodos:

[[Archivo:Error_Euler_implícito_y_modificado_(ventana_grande).jpg |marco|centro|Error entre Euler implícito y Euler modificado]]
[[Archivo:Imageng4.jpg |marco|centro|Error entre Euler explícito y Euler implícito]]
[[Archivo:Error_Euler_explícito_y_modificado_(ventana_grande).jpg ||marco|centro|Error entre Euler explícito y Euler modificado]]

== Conservación de la masa total de contaminante ==

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación <math>u_t-u_{xx}=0</math> , basta con resolver la siguiente integral:

<math>\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0</math>

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración <math>u(x, t)</math> , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.
{{matlab|codigo=

A=U';
integral=trapz(A);
concentracion=U(:,26);
subplot(1,2,1)
plot(t,integral)
xlabel('Tiempo')
ylabel('Masa total de contaminante')
title('Evolución de la cantidad de contaminante en el sistema')
subplot(1,2,2)
plot(t,concentracion)
xlabel('Tiempo')
ylabel('Concentración en el punto medio del tubo')
title('Evolución de la concentración en el punto medio del tubo')

}}

[[Archivo:Integralg4.jpg]]
Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

== Resultados para tiempos grandes ==

Despreciamos <math>u_t</math>:

:<math>
\left\{\begin{matrix}\\-u_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Usamos una T=100 y hacemos los cálculos con el método del trapecio:

{{matlab|codigo=
clear all
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h/4;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
%Solución estacionaria
a=U(length(t),:);
b=U(1,:);
c=U((length(t)-1)/100+1,:);
d=U((length(t)-1)*2/100+1,:);
e=U((length(t)-1)*10/100+1,:);
hold on
plot(x,a,'b','linewidth',2);
plot(x,b,'-y','linewidth',1);
plot(x,c,'-g','linewidth',1);
plot(x,d,'-k','linewidth',1);
plot(x,e,'-r','linewidth',1);
legend('100','0','1','2','10');
xlabel('Posición')
ylabel('Concentracion')
hold off
}}

Y en la gráfica comprobamos las distintas distancias.

[[Archivo:Grafg4.jpg]]

== Cálculo del tiempo necesario para alcanzar el estado estacionario con un error del 5% ==

{{matlab|codigo=
clear all
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1;
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Estado estacionario de la concentración
est=U(length(t),:);
clear A F K L N T U dt h n t u0 uu x
% Repetimos el cálculo y lo paramos al alcanzar un error máximo del 5%
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
tiempo=(b-1)*dt
clear A F K L N T U a b dt h n t u0 uu x
% Cálculo del tiempo necesario para alcanzar el estado estacionario con un
% error del 5% para un intervalo espacial 10 veces más pequeño
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1/10;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(301,1);3*ones(200,1)];
% Discretización temporal
dt=(h*10)^2/2; % Mantenemos el mismo intervalo temporal
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1;
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Estado estacionario de la concentración
est=U(length(t),:);
clear A F K L N T U a b dt h n t u0 uu x
% Comprobamos si varía el tiempo al hacer h 10 veces más pequeño
% Datos del problema
L=5;
T=50;
% Discretización espacial
h=0.1/10;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(1,2)=-2;
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(301,1);3*ones(200,1)];
% Discretización temporal
dt=(h*10)^2/2; % Mantenemos el mismo intervalo temporal
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
nuevo_tiempo=(b-1)*dt
}}
El tiempo que tarda en alcanzar la concentración el estado estacionario con un error del 5% es 8.68 segundos con h=0.1. Si cambiamos la longitud de paso a 0.01, el nuevo tiempo es practicamente invariable, 8.645 segundos.

== Caso de colocar un limpiador en un extremo ==

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right.
</math>

=== Método de diferencias finitas ===

:<math>
\left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right.
</math>

=== Método numérico ===

Al cambiar las condiciones de contorno, el valor estacionario de la concentración <math> u(x,t)</math> es 0 es todos los puntos del tubo, lo que se traduce en que el contaminante desaparece por completo cuando el tiempo es muy grande.

Será preciso calcular el error relativo en cada punto del tubo, obtener el valor absoluto y escoger el máximo valor.

<math>error \ relativo=\cfrac{valor \ real-valor \ estacionario}{valor \ estacionario}=5 </math>%

El tiempo que se tarda en alcanzar el valor del 5% es dependiente del valor que consideramos como estacionario, que en nuestro caso será la concentración correspondiente a T=100 segundos desde que se inicia el ensayo (elegimos este tiempo ya que es suficientemente pequeño el valor de la concentración <math> u(x,t)</math>. El resultado que se obtiene de tiempo para alcanzar dicho error del 5% es 99.49 segundos, quedando reflejado a continuación:

{{matlab|codigo=
% Difusión de una sustancia contaminante con limpiador en el extremo de la
% izquierda
% u_t-u_xx=0, x en (0,L)
% u(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
for n=1:length(t)-1
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
% Valor estacionario de la concentración en el tubo
est=U(length(t),:);
clear A F K L N T U dt h n t u0 uu x xx tt
% Repetimos el cálculo y lo paramos al alcanzar un error máximo del 5%
% Datos del problema
L=5;
T=100;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
% Aproximación de -u_xx
K=2*diag(ones(1,N+1))-diag(ones(1,N),1)-diag(ones(1,N),-1);
K(N+1,N)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N+1,1);
% Vector u0
u0=[zeros(31,1);3*ones(20,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Método del trapecio
A=eye(N+1)+dt/2*K;
U(1,:)=u0';
uu=u0;
n=1;
while max((abs(uu'-est))./est)>0.05
uu=A\((eye(N+1)-dt/2*K)*uu+dt*F);
U(n+1,:)=uu';
n=n+1;
end
% Tiempo invertido
a=size(U);
b=a(1);
tiempo=(b-1)*dt

}}

== Método de Fourier ==

Para la resolución del sistema mediante el método de Fourier:

<math>
\left\{\begin{matrix}\\u_t-u_{xx}=0 \\u_x (0,t)=0 \\u_x (L,t)=0 \\u(x,t)=u_0\end{matrix}\right.
</math>

Buscaremos soluciones del tipo <math>u(x,t)=X(x)T(t)</math>:

<math>XT^{'}-X^{''} T=0</math>:
<math>\cfrac{T'}{T}=\cfrac{X^{''}}{X}=-μ</math>

El problema de autovalores asociado, será:

<math>
\left\{\begin{matrix}\\X^{''}-μX=0 \\X'(0)=0 \\X^{'} (L)=0\end{matrix}\right.
</math>:

<math>μ_k=(\cfrac{kπ}{L})^2</math>:

<math>φ_k=cos⁡(\cfrac{kπ}{L} x)</math>

Por otro lado, la ecuación diferencial que debe cumplir <math>T(t)</math> es:

<math> T^{'}_k (t)-μ_k T_k (t)=0</math>

Resolviéndola obtenemos:

<math> T_k (t)=C_k e^{-μ_k t}</math>

Entonces, la solución aproximada mediante el método de Fourier es la siguiente:

<math>u(x,t)= \sum_{k=0}^K C_k e^{-μt} φ_k </math>

A continuación se muestra el código de matlab que resuelve el problema mediante Fourier usando 1, 3, 5, 10 y 20 términos de la serie.

{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=0
% u_x(L,t)=0
% u(0,t)=u0(x)
% Método de Fourier
% ______________________________________
%
clear all
% Datos del problema
L=5; % Longitud del tubo
a=0; b=5; % Extremos
T=5; % Tiempo final
% Numero de autofunciones de la serie = M
M=20;
% Mallado
% En x --------> Paso espacial = h
h=(b-a)/50;
x=a:h:b;
% En t --------> Paso temporal = j
j=h/4;
t=0:j:T;
[Mx,Mt]=meshgrid(x,t);
f=[zeros(1,31),3*ones(1,20)]; % Valor inicial
U=0;
for k=0:M
p=cos(k*pi*x/L); % Autofunciones
Ck=trapz(x,f.*p)/trapz(x,p.^2);
U=U+Ck*(exp(-(k^2*pi^2/L^2)*Mt).*cos(k*pi*Mx/L));
end
surf(Mx,Mt,U)
xlabel('Espacio');
ylabel('Tiempo');
}}

Para obtener los diferentes resultados basta con introducir en M el número de términos de la serie de Fourier que vayamos a usar. Para M=20 se obtiene un resultado similar a los obtenidos anteriormente con otros métodos

[[Archivo:Fourier20g4.jpg ]]

y si lo calculamos para M=50 se aprecia que la gráfica se suaviza y se asemeja aun más.

[[Archivo:Fourier50.jpg ]]
Los valores negativos que se obtienen son defecto del método.

Las aproximaciones en t=0.5 para 1, 3, 5, 10 y 20 términos de Fourier serán

[[Archivo:comparacionfourier.jpg ]]

== Cambio en las condiciones de frontera ==
El problema pasa a ser
:<math>
\left\{\begin{matrix}\\u_t-u_{xx}=0\\u(5,t)=10sent\\u_x(0,t)=3\\u(x,0)=u_0\end{matrix}\right.
</math>
Este sistema tiene una condición de contorno tipo Dirichlet en el extremo de la varilla <math>x=L</math>.
Esto quiere decir que este extremo está en contacto con una fuente de contaminante siendo la concentración en el instante t en dicho extremo <math>u=10sent</math>

La otra condición es de tipo Neumann en cero <math>-Du_x(0,t)=φ(t)</math>

En este caso <math>φ(t)=3>0</math>, además <math>D=1>0</math>, entonces debe verificarse <math>u_x(0,t)<0</math>. Cogiendo <math>∆x>0</math> y pequeño obtenemos

<math>u(0+∆x,t)-u(0,t)<0</math>

<math>u(0+∆x,t)<u(0,t)</math>

lo que indica que el flujo de contaminante irá hacia la derecha,interpretado de forma que por el extremo x=0 está entrando una cantidad de contaminante en el instante t igual a φ=3.

{{matlab|codigo=
% Difusión de una sustancia contaminante
% u_t-u_xx=0, x en (0,L)
% u_x(0,t)=3
% u(L,t)=f(t)
% u(0,t)=u0(x)
% ______________________________________
%
clear all
% Datos del problema
L=5;
T=10;
% Discretización espacial
h=0.1;
N=L/h;
% Vector de puntos en el espacio
x=0:h:L;
xint=0:h:L-h; % Nodos interiores
% Aproximación de -u_xx
K=2*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
K(1,2)=-2;
K=(1/h^2)*K;
% Vector F
F=zeros(N,1);
F(1)=-3*2/h;
% Vector u0
u0=[zeros(31,1);3*ones(19,1)];
% Discretización temporal
dt=h^2/2;
t=0:dt:T; %Vector de tiempos
% Función f(t)
f=10*sin(t);
% Método de Euler explícito
U(1,:)=[u0' f(1)];
uu=u0;
for n=1:length(t)-1
F(N)=f(n)/h^2;
uu=uu+dt*(-K*uu+F);
U(n+1,:)=[uu' f(n+1)];
end
% Dibujamos
[xx,tt]= meshgrid(x,t);
surf(xx,tt,U)
}}
[[Archivo:ultimog4.jpg]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-05-18T15:45:46Z

Adrián:

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema100.jpg]]

En la gráfica obtenida se muestra la evolución en el tiempo de las poblaciones de presas y predadores. Se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores. De igual forma, se observa que una de las poblaciones de presas, que comienza teniendo una mayor población, tiende a disminuir, mientras que la que empezó con menor población aumenta. Los máximos y los mínimos que se dan en la población de predadores se mantienen constantes en el tiempo.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo:EvolucionPresas1Presas2100.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador. Su forma, sin embargo, no nos permite decir si estamos ante un ecosistema estable.

En la tercera gráfica se aprecia la evolución de la población de un tipo de presas en función de la población del otro tipo de presas.

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo:EvolucionPresa1Presa2300.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

La evolución de las presas en el tiempo continúa la tendencia que se aprecia en los primeros 100 años. Tenderá a extinguirse un tipo de presas debido a que el depredador se alimentará más de ellas, provocando esto que el otro tipo de presas aumente progresivamente su población.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Para esta anchura de paso obtenemos la siguiente tabla de máximos y mínimos asociada:

[[Archivo:TablaEuler12.JPG]]

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Se presenta en la primera grafica de nuevo la evolución en el tiempo de los dos tipos de presas y predador. Se aprecia lo indicado en la gráfica obtenida con Euler modificado: aumentan los dos tipos de presas, teniendo en cuenta que uno de ellos dispone de muchos más ejemplares que el otro, los predadores comienzan entonces a aumentar dado que disponen de más alimento para sobrevivir y esto provoca que disminuyan las presas. Se observa que éstos alcanzan un máximo donde las presas disminuyen hasta su mínimo y por tanto apenas existe alimento. Es en ese momento cuando empiezan a competir en ellos y disminuye su población.

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior, ya que permite que pueda aumentarse la longitud de paso sin acumular tanto error.

En este caso obtenemos la siguiente tabla de máximos y mínimos:

[[Archivo:TablaEuler01.JPG]]

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matriz de parámetros
A=[A1 0 0;0 B1 0;0 0 -C1];
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*[-A2;0;C2]+xx(2)*xx(3)*[0;-B2;C3]);
k2=(A*(xx+(h/2)*k1)+(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(3))*[-A2;0;C2]+(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(3))*[0;-B2;C3]);
k3=(A*(xx+(h/2)*k2)+(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(3))*[-A2;0;C2]+(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(3))*[0;-B2;C3]);
k4=(A*(xx+h*k3)+(xx(1)+h*k3(1))*(xx(3)+h*k3(3))*[-A2;0;C2]+(xx(2)+h*k3(2))*(xx(3)+h*k3(3))*[0;-B2;C3]);
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK2.jpg]]

Es con este método con el que mejor se aprecia la evolución de las especies de presas (llegando a extinguirse una de ellas, la presa tipo 2) y de la especie predadora.

[[Archivo:TablaRK.JPG]]

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Difusión de un contaminante. Grupo 4

2014-05-18T15:45:24Z

Adrián:

{{ Beta }}

{{ TrabajoED | Difusión de una sustancia contaminante. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina '''''difusión''''', proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

== Modelización del problema ==

Consideramos un tubo de longitud <math>L</math> de un cierto material, con sección transversal <math>A</math> constante. La orientamos en la dirección del eje <math>x</math>, de <math>x = 0</math> a <math>x = L</math>. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de <math>x</math> y <math>t</math>. Designamos por <math>u(x, t)</math> la concentración de contaminante en la sección de la varilla que dista <math>x ≥ 0</math> del extremos <math>x = 0</math> cuando ha pasado un tiempo <math>t ≥ 0</math>.

Tomemos un trozo de tubo entre las secciones <math>x</math> y <math>x + ∆x</math>, que designaremos
por <math>[x, x + ∆x]</math>. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por <math>φ(x, t)</math> el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje <math>X</math>. Si <math>φ (x, t) > 0</math> pensamos que el flujo va hacia la derecha, si <math>φ (x, t) <0</math> va hacia la izquierda. El flujo de contaminante en <math>[x, x + ∆x]</math> está dado por <math>φ (x, t)A − φ (x + ∆x, t)A</math>.

Designamos por <math>f(x, t)</math> a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto <math>f(x)=0</math>.

La ''ley de Fick'' (análoga a la de Fourier para la temperatura) dice:
<math>F=-D \frac{∂u}{∂x}</math>

Fick demostró empíricamente que:
<math> φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}</math>

De aquí se deduce la anteriormente citada ley cuando suponemos que <math>∆x</math> tiende a <math>0</math>.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene <math>ux(x,t)>0</math>,
y <math>φ (x, t)</math> será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene <math>ux(x,t)<0</math>,
y <math>φ (x, t)</math> será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en <math>\cfrac{m^2}{s}</math>, que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ''ley de Fick'', la concentración en la sección transversal del tubo que dista x del extremo <math>x = 0</math> cuando ha pasado un tiempo <math>t</math>, <math>u(x, t)</math>, satisface:
<math>u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>x ∈(0,L),t>0</math>

es decir, la ecuación en derivadas parciales:
<math>u_t= \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>u_t-Du_{xx} = f(x,t)</math>

Como hemos señalado anteriormente <math>f(x,t)=0</math> y <math>D=1</math>:
<math>u_t-u_{xx}=0</math>

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

== Planteamiento del sistema ==
El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Supondremos D constante e igual a 1 y L con valor 5.

== Resolución del sistema ==
Método de diferencias finitas con ∆x=0.1 suponiendo que en el instante inicial se verifica:

:<math>
u(x,0)=\left\{\begin{matrix}\\0 & x≤3\\3 & x>3\end{matrix}\right.
</math>

=== Método del trapecio ===
{{matlab|codigo=

}}

=== Método de Euler explícito ===
{{matlab|codigo=

}}

=== Método de Euler implícito ===
{{matlab|codigo=

}}

=== Método de Euler modificado ===
{{matlab|codigo=

}}

== Conservación de la masa total de contaminante ==

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación <math>u_t-u_{xx}=0</math> , basta con resolver la siguiente integral:

<math>\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0</math>

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración <math>u(x, t)</math> , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.
{{matlab|codigo=

}}

Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

== Resultados para tiempos grandes ==

Despreciamos <math>u_t</math>:

:<math>
\left\{\begin{matrix}\\-u_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=0\end{matrix}\right.
</math>

== ==

== Caso de colocar un limpiador en un extremo ==

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right.
</math>

=== Método de diferencias finitas ===

:<math>
\left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right.
</math>

=== Método numérico ===

{{matlab|codigo=

}}

Al cambiar las condiciones de contorno se intuye que para tiempos muy grandes la concentración del tubo tenderá a 0,por tanto se puede decir que el valor estacionario de la concentración en el tubo es 100 segundos, cuando se comprueba que el valor de u es suficientemente pequeño, casi 0. Se consigue un valor estacionario de la concentración con un error del 5% transcurrido un tiempo de 99,49 segundos.

[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Difusión de un contaminante. Grupo 4

2014-05-18T13:08:48Z

Adrián:

{{ Beta }}

{{ TrabajoED | Difusión de una sustancia contaminante. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina '''''difusión''''', proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

== Modelización del problema ==

Consideramos un tubo de longitud <math>L</math> de un cierto material, con sección transversal <math>A</math> constante. La orientamos en la dirección del eje <math>x</math>, de <math>x = 0</math> a <math>x = L</math>. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de <math>x</math> y <math>t</math>. Designamos por <math>u(x, t)</math> la concentración de contaminante en la sección de la varilla que dista <math>x ≥ 0</math> del extremos <math>x = 0</math> cuando ha pasado un tiempo <math>t ≥ 0</math>.

Tomemos un trozo de tubo entre las secciones <math>x</math> y <math>x + ∆x</math>, que designaremos
por <math>[x, x + ∆x]</math>. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por <math>φ(x, t)</math> el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje <math>X</math>. Si <math>φ (x, t) > 0</math> pensamos que el flujo va hacia la derecha, si <math>φ (x, t) <0</math> va hacia la izquierda. El flujo de contaminante en <math>[x, x + ∆x]</math> está dado por <math>φ (x, t)A − φ (x + ∆x, t)A</math>.

Designamos por <math>f(x, t)</math> a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto <math>f(x)=0</math>.

La ''ley de Fick'' (análoga a la de Fourier para la temperatura) dice:
<math>F=-D \frac{∂u}{∂x}</math>

Fick demostró empíricamente que:
<math> φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}</math>

De aquí se deduce la anteriormente citada ley cuando suponemos que <math>∆x</math> tiende a <math>0</math>.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene <math>ux(x,t)>0</math>,
y <math>φ (x, t)</math> será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene <math>ux(x,t)<0</math>,
y <math>φ (x, t)</math> será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en <math>\cfrac{m^2}{s}</math>, que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ''ley de Fick'', la concentración en la sección transversal del tubo que dista x del extremo <math>x = 0</math> cuando ha pasado un tiempo <math>t</math>, <math>u(x, t)</math>, satisface:
<math>u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>x ∈(0,L),t>0</math>

es decir, la ecuación en derivadas parciales:
<math>u_t= \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>u_t-Du_{xx} = f(x,t)</math>

Como hemos señalado anteriormente <math>f(x,t)=0</math> y <math>D=1</math>:
<math>u_t-u_{xx}=0</math>

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

== Planteamiento del sistema ==
El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Supondremos D constante e igual a 1 y L con valor 5.

== Resolución del sistema ==
Método de diferencias finitas con ∆x=0.1 suponiendo que en el instante inicial se verifica:

:<math>
u(x,0)=\left\{\begin{matrix}\\0 & x≤3\\3 & x>3\end{matrix}\right.
</math>

=== Método del trapecio ===
{{matlab|codigo=

}}

=== Método de Euler explícito ===
{{matlab|codigo=

}}

=== Método de Euler implícito ===
{{matlab|codigo=

}}

=== Método de Euler modificado ===
{{matlab|codigo=

}}

== Conservación de la masa total de contaminante ==

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación <math>u_t-u_{xx}=0</math> , basta con resolver la siguiente integral:

<math>\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0</math>

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración <math>u(x, t)</math> , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.
{{matlab|codigo=

}}

Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

== Resultados para tiempos grandes ==

Despreciamos <math>u_t</math>:

:<math>
\left\{\begin{matrix}\\-u_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=0\end{matrix}\right.
</math>

== ==

== Caso de colocar un limpiador en un extremo ==

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right.
</math>

=== Método de diferencias finitas ===

:<math>
\left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right.
</math>

=== Método numérico ===

{{matlab|codigo=

}}

Por tanto se puede decir que el valor estacionario de la concentración en el tubo es 000000. Se consigue un valor estacionario de la concentración con un error del 5% transcurrido un tiempo de 000000.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Difusión de un contaminante. Grupo 4

2014-05-18T13:01:10Z

Adrián:

{{ Beta }}

{{ TrabajoED | Difusión de una sustancia contaminante. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina '''''difusión''''', proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

== Modelización del problema ==

Consideramos un tubo de longitud <math>L</math> de un cierto material, con sección transversal <math>A</math> constante. La orientamos en la dirección del eje <math>x</math>, de <math>x = 0</math> a <math>x = L</math>. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de <math>x</math> y <math>t</math>. Designamos por <math>u(x, t)</math> la concentración de contaminante en la sección de la varilla que dista <math>x ≥ 0</math> del extremos <math>x = 0</math> cuando ha pasado un tiempo <math>t ≥ 0</math>.

Tomemos un trozo de tubo entre las secciones <math>x</math> y <math>x + ∆x</math>, que designaremos
por <math>[x, x + ∆x]</math>. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por <math>φ(x, t)</math> el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje <math>X</math>. Si <math>φ (x, t) > 0</math> pensamos que el flujo va hacia la derecha, <math>si φ (x, t) <0</math> va hacia la izquierda. El flujo de contaminante en <math>[x, x + ∆x]</math> está dado por <math>φ (x, t)A − φ (x + ∆x, t)A</math>.

Designamos por <math>f(x, t)</math> a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto <math>f(x)=0</math>.

La ''ley de Fick'' (análoga a la de Fourier para la temperatura) dice:
<math>F=-D \frac{∂u}{∂x}</math>

Fick demostró empíricamente que:
<math> φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}</math>

De aquí se deduce la anteriormente citada ley cuando suponemos que <math>∆x</math> tiende a <math>0</math>.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene <math>ux(x,t)>0</math>,
y <math>φ (x, t)</math> será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene <math>ux(x,t)<0</math>,
y <math>φ (x, t)</math> será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en <math>\cfrac{m^2}{s}</math>, que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ''ley de Fick'', la concentración en la sección transversal del tubo que dista x del extremo <math>x = 0</math> cuando ha pasado un tiempo <math>t</math>, <math>u(x, t)</math>, satisface:
<math>u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>x ∈(0,L),t>0</math>

es decir, la ecuación en derivadas parciales:
<math>u_t= \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>u_t-Du_{xx} = f(x,t)</math>

Como hemos señalado anteriormente <math>f(x,t)=0</math> y <math>D=1</math>:
<math>u_t-u_{xx}=0</math>

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

== Planteamiento del sistema ==
El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Supondremos D constante e igual a 1 y L con valor 5.

== Resolución del sistema ==
Método de diferencias finitas con ∆x=0.1 suponiendo que en el instante inicial se verifica:

:<math>
u(x,0)=\left\{\begin{matrix}\\0 & x≤3\\3 & x>3\end{matrix}\right.
</math>

=== Método del trapecio ===
{{matlab|codigo=

}}

=== Método de Euler explícito ===
{{matlab|codigo=

}}

=== Método de Euler implícito ===
{{matlab|codigo=

}}

=== Método de Euler modificado ===
{{matlab|codigo=

}}

== Conservación de la masa total de contaminante ==

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación <math>u_t-u_{xx}=0</math> , basta con resolver la siguiente integral:

<math>\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0</math>

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración <math>u(x, t)</math> , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.
{{matlab|codigo=

}}

Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

== Resultados para tiempos grandes ==

Despreciamos <math>u_t</math>:

:<math>
\left\{\begin{matrix}\\-u_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=0\end{matrix}\right.
</math>

== ==

== Caso de colocar un limpiador en un extremo ==

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right.
</math>

=== Método de diferencias finitas ===

:<math>
\left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right.
</math>

=== Método numérico ===

{{matlab|codigo=

}}

Por tanto se puede decir que el valor estacionario de la concentración en el tubo es 000000. Se consigue un valor estacionario de la concentración con un error del 5% transcurrido un tiempo de 000000.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Difusión de un contaminante. Grupo 4

2014-05-18T13:00:45Z

Adrián:

{{ Beta }}

{{ TrabajoED | Difusión de una sustancia contaminante (Grupo 4) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

El proceso irreversible por el cual un grupo de partículas se distribuye de manera uniforme en un medio se denomina '''''difusión''''', proceso estadísticamente predecible en conjunto siendo totalmente aleatorio el movimiento de cada partícula aislada. La difusión de moléculas de una sustancia disuelta en otra será nuestro objeto de estudio, para lo cual disponemos de un tubo que será detallado posteriormente.

== Modelización del problema ==

Consideramos un tubo de longitud <math>L</math> de un cierto material, con sección transversal <math>A</math> constante. La orientamos en la dirección del eje <math>x</math>, de <math>x = 0</math> a <math>x = L</math>. Consideramos una solución compuesta por dos sustancias, de las cuales una de ellas es un contaminante. La sustancia se propaga por el interior de dicho tubo.

Suponemos que el tubo es delgado y que su superficie lateral es aislante, de forma que el contaminante no sale de él, ni puede entrar desde el exterior. Podemos entonces pensar que la concentración es constante a lo largo de cada sección transversal, y ver el tubo como un objeto unidimensional. La concentración va a depender entonces de <math>x</math> y <math>t</math>. Designamos por <math>u(x, t)</math> la concentración de contaminante en la sección de la varilla que dista <math>x ≥ 0</math> del extremos <math>x = 0</math> cuando ha pasado un tiempo <math>t ≥ 0</math>.

Tomemos un trozo de tubo entre las secciones <math>x</math> y <math>x + ∆x</math>, que designaremos
por <math>[x, x + ∆x]</math>. Suponemos que ∆x es una cantidad muy pequeña.

Designamos por <math>φ(x, t)</math> el flujo de contaminante, es decir, la cantidad de contaminante que fluye por unidad de tiempo y unidad de área. Como la superficie lateral del tubo está aislada, solamente habría flujo en la dirección del eje <math>X</math>. Si <math>φ (x, t) > 0</math> pensamos que el flujo va hacia la derecha, <math>si φ (x, t) <0</math> va hacia la izquierda. El flujo de contaminante en <math>[x, x + ∆x]</math> está dado por <math>φ (x, t)A − φ (x + ∆x, t)A</math>.

Designamos por <math>f(x, t)</math> a la concentración de soluto generada por posibles fuentes o sumideros, por unidad de volumen y unidad de tiempo. En el caso de este problema no existen ni fuentes ni sumideros, por lo tanto <math>f(x)=0</math>.

La ''ley de Fick'' (análoga a la de Fourier para la temperatura) dice:
<math>F=-D \frac{∂u}{∂x}</math>

Fick demostró empíricamente que:
<math> φ(x,t)≅-D \cfrac{(u(x+∆x,t)-u(x,t))}{∆x}</math>

De aquí se deduce la anteriormente citada ley cuando suponemos que <math>∆x</math> tiende a <math>0</math>.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es mayor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) > 0, y si ∆x es pequeño, se tiene <math>ux(x,t)>0</math>,
y <math>φ (x, t)</math> será negativo y el flujo irá hacia la izquierda.

Si suponemos que ∆x > 0 y la concentración en un tiempo t es menor en x + ∆x que en x, entonces u(x + ∆x) − u(x, t) < 0, y si ∆x es pequeño, se tiene <math>ux(x,t)<0</math>,
y <math>φ (x, t)</math> será positivo y el flujo de calor irá hacia la derecha.

La constante D es el coeficiente de difusión medido en <math>\cfrac{m^2}{s}</math>, que depende de las propiedades físicas de los compuestos. Este coeficiente D puede depender de x si las condiciones físicas varían en las distintas partes de la sección del tubo. También, en principio, podría depender del tiempo. Supondremos que no, pues en caso contrario el problema en términos matemáticos se complica bastante. La supondremos constante e igual a 1.

Utilizando la ''ley de Fick'', la concentración en la sección transversal del tubo que dista x del extremo <math>x = 0</math> cuando ha pasado un tiempo <math>t</math>, <math>u(x, t)</math>, satisface:
<math>u_t (x,t) = \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>x ∈(0,L),t>0</math>

es decir, la ecuación en derivadas parciales:
<math>u_t= \cfrac{∂}{∂x} Du_x+f(x,t)</math>:
<math>u_t-Du_{xx} = f(x,t)</math>

Como hemos señalado anteriormente <math>f(x,t)=0</math> y <math>D=1</math>:
<math>u_t-u_{xx}=0</math>

Obtenemos de esta forma una ecuación de difusión similar a la del calor.

== Planteamiento del sistema ==
El sistema de ecuaciones que debe satisfacer u(x,t) para que el problema esté bien propuesto es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=u_0\end{matrix}\right.
</math>

Supondremos D constante e igual a 1 y L con valor 5.

== Resolución del sistema ==
Método de diferencias finitas con ∆x=0.1 suponiendo que en el instante inicial se verifica:

:<math>
u(x,0)=\left\{\begin{matrix}\\0 & x≤3\\3 & x>3\end{matrix}\right.
</math>

=== Método del trapecio ===
{{matlab|codigo=

}}

=== Método de Euler explícito ===
{{matlab|codigo=

}}

=== Método de Euler implícito ===
{{matlab|codigo=

}}

=== Método de Euler modificado ===
{{matlab|codigo=

}}

== Conservación de la masa total de contaminante ==

Para comprobar que la masa se conserva a lo largo del tiempo, dado que estamos ante la ecuación <math>u_t-u_{xx}=0</math> , basta con resolver la siguiente integral:

<math>\frac d {dt} \int_0^5 u(x,t)\,dx= \int_0^5 u_{xx}(x,t)\,dx= u_x(x,t)|_0^5=0</math>

que tendrá solución trivial debido a las condiciones frontera.

La primera integral,que hace referencia a la concentración <math>u(x, t)</math> , da como resultado una constante, que ha sido resuelta numéricamente por el método del trapecio.
{{matlab|codigo=

}}

Queda ilustrado en las gráficas perfectamente la evolución de la cantidad de contaminante de sistema, frente al comportamiento de la concentración en el punto medio del tubo a lo largo del tiempo.

== Resultados para tiempos grandes ==

Despreciamos <math>u_t</math>:

:<math>
\left\{\begin{matrix}\\-u_{xx}=0\\u_x(0,t)=0\\u_x(L,t)=0\\u(x,0)=0\end{matrix}\right.
</math>

== ==

== Caso de colocar un limpiador en un extremo ==

Colocar un limpiador en el extremo izquierdo se traduce físicamente en que la concentración en x=0 es ahora nula. Por lo tanto el problema a resolver ahora es el siguiente:

:<math>
\left\{\begin{matrix}\\u_t-Du_{xx}=0\\u(0,t)=0 & t≥0\\u_x(L,t)=0 & t≥0\\u(x,0)=u_0 & xє(0,L)\end{matrix}\right.
</math>

=== Método de diferencias finitas ===

:<math>
\left\{\begin{matrix}\\u'_n(t) + \frac {-u_n + 2u_n(t) + u_{n+1}(t)}{h^2}\\\frac {u_{n+1}(t)-u_{n+1}(t)}{2h}\\u_0(t)=0\\u_n(0)=u_0\end{matrix}\right.
</math>

=== Método numérico ===

{{matlab|codigo=

}}

Por tanto se puede decir que el valor estacionario de la concentración en el tubo es 000000. Se consigue un valor estacionario de la concentración con un error del 5% transcurrido un tiempo de 000000.

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-04T23:39:30Z

Adrián:

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema100.jpg]]

En la gráfica obtenida se muestra la evolución en el tiempo de las poblaciones de presas y predadores. Se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores. De igual forma, se observa que una de las poblaciones de presas, que comienza teniendo una mayor población, tiende a disminuir, mientras que la que empezó con menor población aumenta. Los máximos y los mínimos que se dan en la población de predadores se mantienen constantes en el tiempo.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo:EvolucionPresas1Presas2100.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador. Su forma, sin embargo, no nos permite decir si estamos ante un ecosistema estable.

En la tercera gráfica se aprecia la evolución de la población de un tipo de presas en función de la población del otro tipo de presas.

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo:EvolucionPresa1Presa2300.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

La evolución de las presas en el tiempo continúa la tendencia que se aprecia en los primeros 100 años. Tenderá a extinguirse un tipo de presas debido a que el depredador se alimentará más de ellas, provocando esto que el otro tipo de presas aumente progresivamente su población.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Para esta anchura de paso obtenemos la siguiente tabla de máximos y mínimos asociada:

[[Archivo:TablaEuler12.JPG]]

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Se presenta en la primera grafica de nuevo la evolución en el tiempo de los dos tipos de presas y predador. Se aprecia lo indicado en la gráfica obtenida con Euler modificado: aumentan los dos tipos de presas, teniendo en cuenta que uno de ellos dispone de muchos más ejemplares que el otro, los predadores comienzan entonces a aumentar dado que disponen de más alimento para sobrevivir y esto provoca que disminuyan las presas. Se observa que éstos alcanzan un máximo donde las presas disminuyen hasta su mínimo y por tanto apenas existe alimento. Es en ese momento cuando empiezan a competir en ellos y disminuye su población.

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior, ya que permite que pueda aumentarse la longitud de paso sin acumular tanto error.

En este caso obtenemos la siguiente tabla de máximos y mínimos:

[[Archivo:TablaEuler01.JPG]]

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matriz de parámetros
A=[A1 0 0;0 B1 0;0 0 -C1];
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*[-A2;0;C2]+xx(2)*xx(3)*[0;-B2;C3]);
k2=(A*(xx+(h/2)*k1)+(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(3))*[-A2;0;C2]+(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(3))*[0;-B2;C3]);
k3=(A*(xx+(h/2)*k2)+(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(3))*[-A2;0;C2]+(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(3))*[0;-B2;C3]);
k4=(A*(xx+h*k3)+(xx(1)+h*k3(1))*(xx(3)+h*k3(3))*[-A2;0;C2]+(xx(2)+h*k3(2))*(xx(3)+h*k3(3))*[0;-B2;C3]);
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK2.jpg]]

Es con este método con el que mejor se aprecia la evolución de las especies de presas (llegando a extinguirse una de ellas, la presa tipo 2) y de la especie predadora.

[[Archivo:TablaRK.JPG]]

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-04T23:37:00Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema100.jpg]]

En la gráfica obtenida se muestra la evolución en el tiempo de las poblaciones de presas y predadores. Se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores. De igual forma, se observa que una de las poblaciones de presas, que comienza teniendo una mayor población, tiende a disminuir, mientras que la que empezó con menor población aumenta. Los máximos y los mínimos que se dan en la población de predadores se mantienen constantes en el tiempo.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo:EvolucionPresas1Presas2100.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador. Su forma, sin embargo, no nos permite decir si estamos ante un ecosistema estable.

En la tercera gráfica se aprecia la evolución de la población de un tipo de presas en función de la población del otro tipo de presas.

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo:EvolucionPresa1Presa2300.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

La evolución de las presas en el tiempo continúa la tendencia que se aprecia en los primeros 100 años. Tenderá a extinguirse un tipo de presas debido a que el depredador se alimentará más de ellas, provocando esto que el otro tipo de presas aumente progresivamente su población.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Para esta anchura de paso obtenemos la siguiente tabla de máximos y mínimos asociada:

[[Archivo:TablaEuler12.JPG]]

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Se presenta en la primera grafica de nuevo la evolución en el tiempo de los dos tipos de presas y predador. Se aprecia lo indicado en la gráfica obtenida con Euler modificado: aumentan los dos tipos de presas, teniendo en cuenta que uno de ellos dispone de muchos más ejemplares que el otro, los predadores comienzan entonces a aumentar dado que disponen de más alimento para sobrevivir y esto provoca que disminuyan las presas. Se observa que éstos alcanzan un máximo donde las presas disminuyen hasta su mínimo y por tanto apenas existe alimento. Es en ese momento cuando empiezan a competir en ellos y disminuye su población.

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior, ya que permite que pueda aumentarse la longitud de paso sin acumular tanto error.

En este caso obtenemos la siguiente tabla de máximos y mínimos:

[[Archivo:TablaEuler01.JPG]]

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matriz de parámetros
A=[A1 0 0;0 B1 0;0 0 -C1];
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*[-A2;0;C2]+xx(2)*xx(3)*[0;-B2;C3]);
k2=(A*(xx+(h/2)*k1)+(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(3))*[-A2;0;C2]+(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(3))*[0;-B2;C3]);
k3=(A*(xx+(h/2)*k2)+(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(3))*[-A2;0;C2]+(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(3))*[0;-B2;C3]);
k4=(A*(xx+h*k3)+(xx(1)+h*k3(1))*(xx(3)+h*k3(3))*[-A2;0;C2]+(xx(2)+h*k3(2))*(xx(3)+h*k3(3))*[0;-B2;C3]);
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK2.jpg]]

Es con este método con el que mejor se aprecia la evolución de las especies de presas (llegando a extinguirse una de ellas, la presa tipo 2) y de la especie predadora.

[[Archivo:TablaRK.JPG]]

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:TablaEuler12.JPG

2014-03-04T23:36:50Z

Adrián:

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-04T23:26:42Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema100.jpg]]

En la gráfica obtenida se muestra la evolución en el tiempo de las poblaciones de presas y predadores. Se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores. De igual forma, se observa que una de las poblaciones de presas, que comienza teniendo una mayor población, tiende a disminuir, mientras que la que empezó con menor población aumenta. Los máximos y los mínimos que se dan en la población de predadores se mantienen constantes en el tiempo.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo:EvolucionPresas1Presas2100.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador. Su forma, sin embargo, no nos permite decir si estamos ante un ecosistema estable.

En la tercera gráfica se aprecia la evolución de la población de un tipo de presas en función de la población del otro tipo de presas.

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo:EvolucionPresa1Presa2300.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

La evolución de las presas en el tiempo continúa la tendencia que se aprecia en los primeros 100 años. Tenderá a extinguirse un tipo de presas debido a que el depredador se alimentará más de ellas, provocando esto que el otro tipo de presas aumente progresivamente su población.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Para esta anchura de paso obtenemos la siguiente tabla de máximos y mínimos asociada:

[[Archivo:TablaEuler1.JPG]]

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Se presenta en la primera grafica de nuevo la evolución en el tiempo de los dos tipos de presas y predador. Se aprecia lo indicado en la gráfica obtenida con Euler modificado: aumentan los dos tipos de presas, teniendo en cuenta que uno de ellos dispone de muchos más ejemplares que el otro, los predadores comienzan entonces a aumentar dado que disponen de más alimento para sobrevivir y esto provoca que disminuyan las presas. Se observa que éstos alcanzan un máximo donde las presas disminuyen hasta su mínimo y por tanto apenas existe alimento. Es en ese momento cuando empiezan a competir en ellos y disminuye su población.

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior, ya que permite que pueda aumentarse la longitud de paso sin acumular tanto error.

En este caso obtenemos la siguiente tabla de máximos y mínimos:

[[Archivo:TablaEuler01.JPG]]

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matriz de parámetros
A=[A1 0 0;0 B1 0;0 0 -C1];
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*[-A2;0;C2]+xx(2)*xx(3)*[0;-B2;C3]);
k2=(A*(xx+(h/2)*k1)+(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(3))*[-A2;0;C2]+(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(3))*[0;-B2;C3]);
k3=(A*(xx+(h/2)*k2)+(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(3))*[-A2;0;C2]+(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(3))*[0;-B2;C3]);
k4=(A*(xx+h*k3)+(xx(1)+h*k3(1))*(xx(3)+h*k3(3))*[-A2;0;C2]+(xx(2)+h*k3(2))*(xx(3)+h*k3(3))*[0;-B2;C3]);
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK2.jpg]]

Es con este método con el que mejor se aprecia la evolución de las especies de presas (llegando a extinguirse una de ellas, la presa tipo 2) y de la especie predadora.

[[Archivo:TablaRK.JPG]]

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-04T23:09:59Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema100.jpg]]

En la gráfica se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo:EvolucionPresas1Presas2100.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador.

<big><big><big><big><big><big><big><big><big><big>FALTA EXPLICAR LA TERCERA!!!</big></big></big></big></big></big></big></big></big></big>

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo:EvolucionPresa1Presa2300.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

[[Archivo:TablaEuler1.JPG]]

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior.

[[Archivo:TablaEuler01.JPG]]

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matriz de parámetros
A=[A1 0 0;0 B1 0;0 0 -C1];
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*[-A2;0;C2]+xx(2)*xx(3)*[0;-B2;C3]);
k2=(A*(xx+(h/2)*k1)+(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(3))*[-A2;0;C2]+(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(3))*[0;-B2;C3]);
k3=(A*(xx+(h/2)*k2)+(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(3))*[-A2;0;C2]+(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(3))*[0;-B2;C3]);
k4=(A*(xx+h*k3)+(xx(1)+h*k3(1))*(xx(3)+h*k3(3))*[-A2;0;C2]+(xx(2)+h*k3(2))*(xx(3)+h*k3(3))*[0;-B2;C3]);
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK2.jpg]]

La situación que nos muestra esta gráfica es diferente a las anteriores. Se aprecia que la población <math>x_2</math> tiende a extinguirse. Esto tiene dos consecuencias, el aumento de los predadores y la estabilizacion de la presa <math>x_1</math> que hasta el momento de la desaparicion de <math>x_2</math> había estado creciendo

[[Archivo:TablaRK.JPG]]

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:TablaRK.JPG

2014-03-04T23:09:46Z

Adrián:

Archivo:TablaEuler1.JPG

2014-03-04T23:06:27Z

Adrián:

Archivo:TablaEuler01.JPG

2014-03-04T23:04:00Z

Adrián:

Archivo:EvolucionEcosistemaRK2.jpg

2014-03-04T22:56:59Z

Adrián:

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-04T22:50:21Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema100.jpg]]

En la gráfica se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores.

[[Archivo:EvolucionPresas1Predadores100.jpg]] [[Archivo:EvolucionPresas2Predadores100.jpg]] [[Archivo:EvolucionPresas1Presas2100.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador.

<big><big><big><big><big><big><big><big><big><big>FALTA EXPLICAR LA TERCERA!!!</big></big></big></big></big></big></big></big></big></big>

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Matrices de coeficientes del sistema
A=[A1 0 0;0 B1 0;0 0 -C1];
B=[-A2 0 C2]';
C=[0 -B2 C3]';
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
xx=x0;
for n=1:N
k1=(A*xx+xx(1)*xx(3)*B+xx(2)*xx(3)*C);
k2=(A*(xx+h*k1)+(xx(1)+h*k1(1))*(xx(3)+h*k1(3))*B+(xx(2)+h*k1(2))*(xx(3)+h*k1(3))*C);
xx=xx+h*(k1+k2)/2;
x1(n+1)=xx(1);
x2(n+1)=xx(2);
x3(n+1)=xx(3);
end
hold on
plot(t,x1,'b')
plot(t,x2,'g')
plot(t,x3,'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x1,x3,'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x2,x3,'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x1,x2,'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema300.jpg]] [[Archivo:EvolucionPresa1Predador300.jpg]] [[Archivo:EvolucionPresa2Predador300.jpg]] [[Archivo:EvolucionPresa1Presa2300.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior.

<big><big><big><big><big><big><big><big><big><big>FALTA METER LAS TABLAS!!!</big></big></big></big></big></big></big></big></big></big>

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
k3=[A1*(xx(1)+(h/2)*k2(1))-A2*(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(1)),B1*(xx(2)+(h/2)*k2(2))-B2*(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(2)),-C1*(xx(3)+(h/2)*k2(3))+C2*(xx(1)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))+C3*(xx(2)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))]';
k4=[A1*(xx(1)+h*k3(1))-A2*(xx(1)+h*k3(1))*(xx(3)+h*k3(1)),B1*(xx(2)+h*k3(2))-B2*(xx(2)+h*k3(2))*(xx(3)+h*k3(2)),-C1*(xx(3)+h*k3(3))+C2*(xx(1)+h*k3(3))*(xx(3)+h*k3(3))+C3*(xx(2)+h*k3(3))*(xx(3)+h*k3(3))]';
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK.jpg]]

La situación que nos muestra esta gráfica es diferente a las anteriores. Se aprecia que la población <math>x_2</math> tiende a extinguirse. Esto tiene dos consecuencias, el aumento de los predadores y la estabilizacion de la presa <math>x_1</math> que hasta el momento de la desaparicion de <math>x_2</math> había estado creciendo

<big><big><big><big><big><big><big><big><big><big>FALTA METER LOS MÁXIMOS Y MÍNIMOS!!!</big></big></big></big></big></big></big></big></big></big>

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:EvolucionPresa1Presa2300.jpg

2014-03-04T22:50:08Z

Adrián:

Archivo:EvolucionPresa2Predador300.jpg

2014-03-04T22:49:32Z

Adrián:

Archivo:EvolucionPresa1Predador300.jpg

2014-03-04T22:48:55Z

Adrián:

Archivo:EvolucionEcosistema300.jpg

2014-03-04T22:48:07Z

Adrián:

Archivo:EvolucionPresas1Presas2100.jpg

2014-03-04T22:40:21Z

Adrián:

Archivo:EvolucionPresas2Predadores100.jpg

2014-03-04T22:39:26Z

Adrián:

Archivo:EvolucionPresas1Predadores100.jpg

2014-03-04T22:38:32Z

Adrián:

Archivo:EvolucionEcosistema100.jpg

2014-03-04T22:37:38Z

Adrián:

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-02T13:44:12Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;

x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
xx=xx+(h/2)*(k1+k2);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema.jpg]]

En la gráfica se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores.

[[Archivo:EvolucionPresas1Predadores.jpg]] [[Archivo:EvolucionPresas2Predadores.jpg]] [[Archivo:EvolucionPresas1Presas2.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador.

<big><big><big><big><big><big><big><big><big><big>FALTA EXPLICAR LA TERCERA!!!</big></big></big></big></big></big></big></big></big></big>

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
xx=xx+(h/2)*(k1+k2);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema2.jpg]] [[Archivo:EvolucionPresa1Predador2.jpg]] [[Archivo:EvolucionPresa2Predador2.jpg]] [[Archivo:EvolucionPresa1Presa22.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior.

<big><big><big><big><big><big><big><big><big><big>FALTA METER LAS TABLAS!!!</big></big></big></big></big></big></big></big></big></big>

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
k3=[A1*(xx(1)+(h/2)*k2(1))-A2*(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(1)),B1*(xx(2)+(h/2)*k2(2))-B2*(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(2)),-C1*(xx(3)+(h/2)*k2(3))+C2*(xx(1)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))+C3*(xx(2)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))]';
k4=[A1*(xx(1)+h*k3(1))-A2*(xx(1)+h*k3(1))*(xx(3)+h*k3(1)),B1*(xx(2)+h*k3(2))-B2*(xx(2)+h*k3(2))*(xx(3)+h*k3(2)),-C1*(xx(3)+h*k3(3))+C2*(xx(1)+h*k3(3))*(xx(3)+h*k3(3))+C3*(xx(2)+h*k3(3))*(xx(3)+h*k3(3))]';
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK.jpg]]

La situación que nos muestra esta gráfica es diferente a las anteriores. Se aprecia que la población <math>x_2</math> tiende a extinguirse. Esto tiene dos consecuencias, el aumento de los predadores y la estabilizacion de la presa <math>x_1</math> que hasta el momento de la desaparicion de <math>x_2</math> había estado creciendo

<big><big><big><big><big><big><big><big><big><big>FALTA METER LOS MÁXIMOS Y MÍNIMOS!!!</big></big></big></big></big></big></big></big></big></big>

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-02T13:41:39Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;

x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
xx=xx+(h/2)*(k1+k2);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema.jpg]]

En la gráfica se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores.

[[Archivo:EvolucionPresas1Predadores.jpg]] [[Archivo:EvolucionPresas2Predadores.jpg]] [[Archivo:EvolucionPresas1Presas2.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador.

La tercera...

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
xx=xx+(h/2)*(k1+k2);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema2.jpg]] [[Archivo:EvolucionPresa1Predador2.jpg]] [[Archivo:EvolucionPresa2Predador2.jpg]] [[Archivo:EvolucionPresa1Presa22.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]] [[Archivo:EvolucionPresa1Predador01.jpg]] [[Archivo:EvolucionPresa2Predador01.jpg]]

Este resultado se aproxima más al que hemos obtenido con Euler modificado, pero el método de Euler es un método inestable para sistemas. Por ello es preferible usar el método anterior.

<big><big><big><big><big><big><big><big><big><big>FALTA METER LAS TABLAS!!!</big></big></big></big></big></big></big></big></big></big>

=== Método de Runge-Kutta ===

Resolvemos mediante el método de Runge-Kutta con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 500 años y tomando <math>A_1=0.4</math>, <math>A_2=0.7</math>, <math>B_1=0.2</math>, <math>B_2=0.4</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.03</math>. Dadas las condiciones iniciales: <math>p_{0}=3.5</math> millones de presas <math>x_1</math>, <math>q_{0}=0.2</math> millones de presas <math>x_2</math> y <math>d_{0}=0.4</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Runge Kutta
clear all
%Datos
t0=0;
tf=500;
x0= [3.5 0.2 0.4]';
A1=0.4;
A2=0.7;
B1=0.2;
B2=0.4;
C1=0.37;
C2=0.04;
C3=0.03;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
k3=[A1*(xx(1)+(h/2)*k2(1))-A2*(xx(1)+(h/2)*k2(1))*(xx(3)+(h/2)*k2(1)),B1*(xx(2)+(h/2)*k2(2))-B2*(xx(2)+(h/2)*k2(2))*(xx(3)+(h/2)*k2(2)),-C1*(xx(3)+(h/2)*k2(3))+C2*(xx(1)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))+C3*(xx(2)+(h/2)*k2(3))*(xx(3)+(h/2)*k2(3))]';
k4=[A1*(xx(1)+h*k3(1))-A2*(xx(1)+h*k3(1))*(xx(3)+h*k3(1)),B1*(xx(2)+h*k3(2))-B2*(xx(2)+h*k3(2))*(xx(3)+h*k3(2)),-C1*(xx(3)+h*k3(3))+C2*(xx(1)+h*k3(3))*(xx(3)+h*k3(3))+C3*(xx(2)+h*k3(3))*(xx(3)+h*k3(3))]';
xx=xx+(h/6)*(k1+2*k2+2*k3+k4);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
}}

[[Archivo:EvolucionEcosistemaRK.jpg]]

La situación que nos muestra esta gráfica es diferente a las anteriores. Se aprecia que la población <math>x_2</math> tiende a extinguirse. Esto tiene dos consecuencias, el aumento de los predadores y la estabilizacion de la presa <math>x_1</math> que hasta el momento de la desaparicion de <math>x_2</math> había estado creciendo

<big><big><big><big><big><big><big><big><big><big>FALTA METER LOS MÁXIMOS Y MÍNIMOS!!!</big></big></big></big></big></big></big></big></big></big>

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:EvolucionEcosistemaRK.jpg

2014-03-02T13:38:33Z

Adrián:

Archivo:EvolucionPresa2Predador01.jpg

2014-03-02T13:26:20Z

Adrián:

Archivo:EvolucionPresa1Predador01.jpg

2014-03-02T13:25:36Z

Adrián:

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-02T13:24:53Z

Adrián:

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Sandra Carrillo del Cura 81, Sergio Castillo Herrero 85, Andrea García Prieto 171, Patricia González Peinado 198, Adrián Salas Calvo 385 }}

== Planteamiento e interpretación ==

Se conoce como '''modelo de Volterra-Lotka''' el modelo matemático que describe la lucha constante por la supervivencia entre dos especies que viven en un mismo hábitat siendo una de ellas el alimento de la otra. Se sabe que la especie predadora se extinguiría si no dispusiera de la presa y que a su vez las presas crecerían (exponencialmente) sin la presencia de los predadores.

Asumiendo las siguientes hipótesis:

* Las tres poblaciones son homogéneas

* Hay una cantidad suficiente de comida disponible para la alimentación de las poblaciones de las presas <math>x_1(t)</math> y <math>x_2(t)</math>, con lo cual sus respectivas tasas de natalidad siguen una ley malthusiana o exponencial:

:<math>A_1x_1</math>
:<math>B_1x_2</math>

:Siendo <math>A_1</math> y <math>B_1</math> constantes posistivas

 

* La alimentación de la población de predadores <math>x_3(t)</math> depende completamente de las presas, por lo que su tasa de natalidad depende de las iteraciones con ellas:

:<math>C_2 x_1 x_3 + C_3 x_2 x_3</math>

:Siendo <math>C_2</math> y <math>C_3</math> constantes positivas

 

* La tasa de mortalidad de las poblaciones de presas dependerán de las iteraciones presa-predador:

:<math>A_2 x_1 x_3</math>
:<math>B_2 x_2 x_3</math>

:Siendo <math>A_2</math> y <math>B_2</math> constantes positivas

 

* Ya que sin alimentación la población iría desapareciendo, la tasa de mortalidad de los predadores será:

:<math>C_1 x_3</math>

:Siendo <math>C_1</math> una constante positiva

 

Este modelo explica la evolución conjunta de las especies mediante el siguiente sistema de ecuaciones:

 

:<math>
\left\{\begin{matrix}\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\\frac{\mathrm{d}x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\ x_1(0)=p_{0} ; x_2(0)=q_{0} ; x_3(0)=d_{0}\end{matrix}\right.
</math>

:Siendo <math>t>0</math>

== Resolución numérica ==

=== Método de Euler modificado ===

Resolvemos mediante el método de Euler modificado con <math>h=0.1</math> en un intervalo de tiempo de 0 a 100 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=2</math> millones de presas <math>x_1</math>, <math>q_{0}=1.4</math> millones de presas <math>x_2</math> y <math>d_{0}=1</math> millón de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=100;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;

x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
xx=xx+(h/2)*(k1+k2);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema.jpg]]

En la gráfica se aprecia que cuando abunda la especie presa, la predadora tiene mucho alimento, entonces aumenta su población disminuyendo las presas. Si la población de predadores es muy numerosa comen demasiadas presas y disminuye esta población de tal forma que los predadores llegan a pasar hambre y comienza a disminuir su población. Al disminuir el número de predadores, las presas se encuentran relativamente seguras y su población vuelve a crecer. Se observa un ciclo de aumentos y disminuciones interrelacionadas de las poblaciones de las presas y de los predadores.

[[Archivo:EvolucionPresas1Predadores.jpg]] [[Archivo:EvolucionPresas2Predadores.jpg]] [[Archivo:EvolucionPresas1Presas2.jpg]]

Aquí se ha representado la evolución de una población en función de otra, a diferencia de la anterior que muestra la evolución en el tiempo.

Las dos primeras muestran la evolución de los predadores en función de cada una de las presas, obteniéndose una gráfica con forma cíclica que corrobora la relación presa-predador.

La tercera...

Tomando ahora un intervalo de 0 a 300 años, el código para matlab será:

{{matlab|codigo=
%Euler modificado
clear all
%Datos
t0=0;
tf=300;
x0= [2 1.4 1]';
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
%Datos discretización
h=0.1;
N=(tf-t0)/h;
%Vector de tiempos
t=t0:h:tf;
x=zeros(3,N+1);
%Inicialización
x(:,1)=x0;
xx=x0;
for n=1:N
k1=[A1*xx(1)-A2*xx(1)*xx(3),B1*xx(2)-B2*xx(2)*xx(3),-C1*xx(3)+C2*xx(1)*xx(3)+C3*xx(2)*xx(3)]';
k2=[A1*(xx(1)+(h/2)*k1(1))-A2*(xx(1)+(h/2)*k1(1))*(xx(3)+(h/2)*k1(1)),B1*(xx(2)+(h/2)*k1(2))-B2*(xx(2)+(h/2)*k1(2))*(xx(3)+(h/2)*k1(2)),-C1*(xx(3)+(h/2)*k1(3))+C2*(xx(1)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))+C3*(xx(2)+(h/2)*k1(3))*(xx(3)+(h/2)*k1(3))]';
xx=xx+(h/2)*(k1+k2);
x(:,n+1)=xx;
end
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'r')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'k')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(4)
plot(x(1,:),x(2,:),'c')
title('Evolución de las presas tipo 1 frente a las presas tipo 2')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Presas tipo 2 (millones de individuos')
}}

[[Archivo:EvolucionEcosistema2.jpg]] [[Archivo:EvolucionPresa1Predador2.jpg]] [[Archivo:EvolucionPresa2Predador2.jpg]] [[Archivo:EvolucionPresa1Presa22.jpg]]

En estas gráficas se aprecia de forma más clara lo indicado anteriormente.

=== Método de Euler ===

Resolvemos mediante el método de Euler con <math>h=1</math> y con <math>h=0.1</math>, en un intervalo de tiempo de 0 a 300 años y tomando <math>A_1=0.35</math>, <math>A_2=0.6</math>, <math>B_1=0.3</math>, <math>B_2=0.5</math>, <math>C_1=0.37</math>, <math>C_2=0.04</math> y <math>C_3=0.035</math>. Dadas las condiciones iniciales <math>p_{0}=0.8</math> millones de presas <math>x_1</math>, <math>q_{0}=2.4</math> millones de presas <math>x_2</math> y <math>d_{0}=0.2</math> millones de predadores, el código en Matlab es el siguiente:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistemaEuler.jpg]] [[Archivo:EvolucionPresa1PredadorEuler.jpg]] [[Archivo:EvolucionPresa2PredadorEuler.jpg]]

Hemos obtenido unos resultados muy disparados que no muestran las trayectorias con claridad. Esto es debido a que la anchura de paso h es muy grande y se produce un error importante en la aproximación.

Tomando <math>h=0.1</math> el codigo será:

{{matlab|codigo=
%Euler
clear all
%datos
t0=0;
tf=300;
A1=0.35;
A2=0.6;
B1=0.3;
B2=0.5;
C1=0.37;
C2=0.04;
C3=0.035;
x0=[0.8;2.4;0.2];
%datos discretización
h=0.1;
N=(tf-t0)/h;
t=t0:h:tf;
%vector soluciones
xx=x0;
x(:,1)=xx;
%Coeficientes de la parte lineal del sistema
A=[A1,0,0;0,B1,0;0,0,-C1];
for n=1:N
xx=xx+h*(A*xx+[-A2;0;C2]*xx(1)*xx(3)+[0;-B2;C3]*xx(2)*xx(3));
x(:,n+1)=xx;
end
figure(1)
hold on
plot(t,x(1,:),'b')
plot(t,x(2,:),'g')
plot(t,x(3,:),'r')
title('Evolución del ecosistema')
xlabel('Tiempo (años)')
ylabel('Población (millones de individuos)')
legend('Presa tipo 1','Presa tipo 2','Predador')
hold off
figure(2)
plot(x(1,:),x(3,:),'k')
title('Evolución de las presas tipo 1 frente a los predadores')
xlabel('Presas tipo 1 (millones de individuos)')
ylabel('Predadores (millones de individuos')
figure(3)
plot(x(2,:),x(3,:),'r')
title('Evolución de las presas tipo 2 frente a los predadores')
xlabel('Presas tipo 2 (millones de individuos)')
ylabel('Predadores (millones de individuos')
}}

[[Archivo:EvolucionEcosistema01.jpg]]

<math> </math>

[[Categoría:Ecuaciones Diferentiales]]
[[Categoría:ED13/14]]
[[Categoría:Trabajos 2013-14]]

Archivo:EvolucionEcosistema01.jpg

2014-03-02T13:24:39Z

Adrián:

Archivo:EvolucionPresa2PredadorEuler.jpg

2014-03-02T13:19:40Z

Adrián:

Archivo:EvolucionPresa1PredadorEuler.jpg

2014-03-02T13:13:23Z

Adrián:

Archivo:EvolucionEcosistemaEuler.jpg

2014-03-02T13:12:09Z

Adrián:

Archivo:EvolucionPresa1Presa22.jpg

2014-03-02T13:01:24Z

Adrián:

Archivo:EvolucionPresa2Predador2.jpg

2014-03-02T13:00:35Z

Adrián:

Archivo:EvolucionPresa1Predador2.jpg

2014-03-02T12:59:51Z

Adrián:

Archivo:EvolucionEcosistema2.jpg

2014-03-02T12:57:13Z

Adrián:

Archivo:EvolucionEcosistema.jpg

2014-03-02T12:48:03Z

Adrián: Este gráfico muestra la evolución en el tiempo de las presas tipo 1 (en azul), las presas tipo 2 (en verde) y los predadores (en rojo)

Este gráfico muestra la evolución en el tiempo de las presas tipo 1 (en azul), las presas tipo 2 (en verde) y los predadores (en rojo)

Archivo:EvolucionPresas1Presas2.jpg

2014-03-02T12:41:06Z

Adrián: Este gráfico muestra la evolución de las presas tipo 1 frente a las presas tipo 2

Este gráfico muestra la evolución de las presas tipo 1 frente a las presas tipo 2

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-02T12:34:22Z

Adrián:

Archivo:EvolucionPresas2Predadores.jpg

2014-03-02T12:33:02Z

Adrián: Este gráfico muestra la evolución de las presas tipo 2 frente a los predadores

Este gráfico muestra la evolución de las presas tipo 2 frente a los predadores

Archivo:EvolucionPresas1Predadores.jpg

2014-03-02T12:30:13Z

Adrián: Este gráfico muestra la evolución de las presas tipo 1 frente los predadores

Este gráfico muestra la evolución de las presas tipo 1 frente los predadores

Archivo:EvolucionPresa1Predador.jpg

2014-03-02T12:23:20Z

Adrián: Este gráfico muestra la evolución de las presas tipo 1 frente los predadores

Este gráfico muestra la evolución de las presas tipo 1 frente los predadores

Modelo Predador-Presa. Dinámica de poblaciones (Grupo 4)

2014-03-02T11:40:07Z

Adrián: Página creada con «{{ Beta }} {{ TrabajoED | Modelo predador-presa. Dinámica de poblaciones. Grupo 4 | Ecuaciones Diferenciales|:Categoría:ED13/1...»

Modelos de mezclas (grupo 5)

2014-02-02T17:23:57Z

Adrián:

{{Trabajo|Modelo de mezclas|[[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:Trabajos 2012-13|2012-13]]}}
== Introducción ==
Tenemos dos pantanos '''A''' y '''B''', con '''100 <math>Hm^3</math>''' de agua cada uno, que están unidos por una presa que deja pasar agua de A a B. El pantano A recibe '''3 <math>Hm^3/día</math>''' de agua limpia proveniente de ríos y el B '''1.5 <math>Hm^3/día</math>'''. Para mantener el nivel de los pantanos estable la presa de A a B deja pasar una media de '''3 <math>Hm^3/día</math>''' mientras que la presa al final de B desaloja '''4.5 <math>Hm^3/día</math>'''. Se produce un vertido tóxico en el pantano A que deja '''20 toneladas''' de un cierto contaminante.
Supongamos que se dan las siguientes hipótesis::

* El contaminante está disuelto de forma homogénea en el agua de los pantanos.
* Al entrar o salir agua en un pantano, esta se mezcla con el agua del pantano de forma inmediata creando una mezcla homogénea.
* La variación de contaminante en un lago es la diferencia entre el contaminante que entra y sale en el lago, es decir si denotamos por <math> x_A (t) </math> la cantidad de contaminante en el lago A, se tiene:: <math> \frac{dx_A}{dt} = velocidad de entrada – velocidad de salida</math>
[[Archivo:PantanosAyB.jpg|800px|centro|]]

== Planteamiento ==

=== Sin plan de limpieza ===

Sean <math> x_A </math> y <math> x_B </math> las cantidades de contaminante en cada uno de los lagos, la cantidad de contaminante en un instante de tiempo determinado se obtiene como diferencia entre la velocidad de entrada y salida de agua contaminada. En el lago A sólo entra agua limpia mientras que la velocidad de salida se obtiene como el producto de la concentración de contaminante en el lago A <math> \frac{X_A}{100} </math> por el caudal de agua que sale '''3 <math>^3/día </math>'''. Así mismo, en el lago B la velocidad de entrada de agua contaminada es la misma que la velocidad de salida de agua contaminada del lago A, ya que el lago A vierte el agua contaminada en B, y la velocidad de salida del agua de B es el producto de su concentración <math> \frac{X_B}{100} </math> por el caudal de agua que sale '''4.5 <math> Hm^3/día </math>'''.

Para modelizar el problema, se plantea el siguiente sistema dos ecuaciones diferenciales de primer orden:

<math>
\left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d} t}=-\frac{3X_A}{100}\\\frac{\mathrm{d} X_B}{\mathrm{d} t}=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right.
</math>

=== Con plan de limpieza ===

A continuación se propone activar un plan de limpieza basado en bombear '''1<math> Hm^3/día </math>''' de agua del pantano B al A, por lo que para mantener el volumen de agua entre los lagos constante se tiene que modificar el flujo entre ambos.

De este modo, si al lago A le llegan '''3<math>Hm^3/día</math>''' provenientes de los ríos y '''1<math>Hm^3/día</math>''' proveniente de B, el lago A deberá verter a B '''4<math>Hm^3/día </math>''' en lugar de los ''''3 <math> Hm^3/día </math>''' que estaba vertiendo antes de activar el plan de limpieza. El lago B recibe '''1.5 <math> Hm^3/día </math>''' provenientes de los ríos y, como se ha visto, '''4 <math> Hm^3/día </math>''' provenientes del lago A. Tras activar el plan de limpieza bombea '''1 <math> Hm^3/día </math>''' a A, por lo que el vertido que desaloja B aguas abajo se mantendrá constante igual a '''4.5 <math> Hm^3/día </math>''' para mantener el volumen de agua.
Así se plantea un sistema de ecuaciones diferenciales similar al anterior::
<math>
\left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d}t}=-\frac{4X_A}{100}+\frac{1X_B}{100}\\\frac{\mathrm{d}X_B}{\mathrm{d} t}=\frac{4X_A}{100}-\frac{5.5X_B}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right.
</math>

===Unión de un tercer pantano===
[[Archivo:PantanosA,ByC.jpg|800px|centro|]]
Si hubiese un tercer pantano, al que se le llamará C, unido al anterior que desalojase '''6 <math>Hm^3/día </math>''' y recibiese '''1.5 <math>Hm^3/día </math>''', la velocidad de entrada de agua sería igual a la de salida de B, es decir, '''4.5 <math> Hm^3/día </math>''' ,mientras que la de salida serían estos '''4.5 <math>Hm^3/día </math>''' más '''1.5 <math> Hm^3/día </math>''' provenientes de los ríos, dando lugar a los '''6 <math>Hm^3/día </math>''' para que se mantenga constante el volumen. Así, el sistema de ecuaciones diferenciales quedará de la siguiente forma::

<math>
\left\{\begin{matrix}\frac{\mathrm{d}X_A}{\mathrm{d}t}=\frac{-3X_A}{100}\\\frac{\mathrm{d}X_B}{\mathrm{d} t}=\frac{3X_A}{100}-\frac{4.5X_B}{100}\\\frac{\mathrm{d}X_C}{\mathrm{d} t}=\frac{4.5X_B}{100}-\frac{6X_C}{100}\\X_A(t_{0})=20, \quad X_B(t_{0})=0\end{matrix}\right.
</math>

==Resolución numérica==

A partir de este momento nos olvidamos del pantano C. Para resolver estos sistemas se utilizará un programade análisis numérico (en este caso Matlab®).

=== Método de Euler ===

Se utilizan los siguientes códigos::

'''Sin plan de limpieza:'''
{{matlab|codigo=
clear all

N=100000; %número de iteraciones
t0=0;tN=325; %tiempo inicial y final
h=(tN-t0)/N; %paso: espacio entre dos iteraciones
x0=[20 0]'; %condiciones iniciales
A=[-0.03 0;0.03 -0.045];
x=x0;
xA(1)=x(1);
xB(1)=x(2);

for n=1:N
x=x+h*A*x;
xA(n+1)=x(1);
xB(n+1)=x(2);
end
y=t0:h:tN;

%dibujamos la gráfica
hold on
plot(y,xA,'b') %contaminante en pantano A: azul
plot(y,xB,'r') %contaminante en pantano B: rojo
title('Sin plan de limpieza')
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
hold off
}}

'''Con plan de limpieza:'''
{{matlab|codigo=
clear all

N=100000; %número de iteraciones
t0=0;tN=325; %tiempo inicial y final
h=(tN-t0)/N; %paso: espacio entre dos iteraciones
x0=[20 0]'; %condiciones iniciales
B=[-0.04 0.01;0.04 -0.055];
X=x0;
Xa(1)=X(1);
Xb(1)=X(2);

for n=1:N
X=X+h*B*X;
Xa(n+1)=X(1);
Xb(n+1)=X(2);
end
y=t0:h:tN;

%dibujamos la gráfica
hold on
plot(y,Xa,'g') %contaminante en pantano A: verde
plot(y,Xb,'m') %contaminante en pantano B: magenta
title('Con plan de limpieza')
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
hold off
}}

Si queremos ver la diferencia de tiempo que tarda en desaparecer la mitad de contaminante en el pantano A activando el plan de limpieza o sin activarlo no tenemos más que observar las gráficas.

{|
|-
| [[Archivo:Eulersinplan.jpg|500px|left|]] || [[Archivo:Eulerconplan.jpg|500px|left|]]
|}

Trazando una horizontal desde 10 toneladas cortará a la gráfica sin limpieza en el tiempo <math> t_{mitad1}=23 días </math> y a la gráfica con limpieza en <math> t_{mitad2}=18 días </math>, por lo que la diferencia de tiempos será <math> t_{mitad1}- t_{mitad2}=23-18=5 días</math>.

Si hacemos lo mismo para ver cuánto tarda en desparecer la tercera parte <math> t_{tercera1}-t_{tercera2}=13-10=3 días</math>.

=== Método de Runge-Kutta ===

Aplicando el método de Runge-Kutta de cuarto orden se utilizan el siguiente código:

{{matlab|codigo=
clear all

N=10000; %número de interacciones
t0=0; tN=325; %tiempo inicial y final
x0=[20 0]'; %condiciones iniciales
h=(tN-t0)/N; %paso: espacio entre dos iteraciones

%sin plan de limpieza
A=[-0.03 0;0.03 -0.045];
x=x0;
xA(1)=x(1);
xB(1)=x(2);
for n=1:N
k1=A*x;
k2=A*(x+(1/2*k1*h));
k3=A*(x+(1/2*k2*h));
k4=A*(x+(k3*h));
x=x+h/6*(k1+2*k2+2*k3+k4);
xA(n+1)=x(1);
xB(n+1)=x(2);
end

%Con plan de limpieza
B=[-0.04 0.01;0.04 -0.055];
X=x0;
Xa(1)=X(1);
Xb(1)=X(2);
for n=1:N
k1=B*X;
k2=B*(X+(1/2*k1*h));
k3=B*(X+(1/2*k2*h));
k4=B*(X+(k3*h));
X=X+h/6*(k1+2*k2+2*k3+k4);
Xa(n+1)=X(1);
Xb(n+1)=X(2);
end

y=t0:h:tN;

%dibujamos las dos gráficas
hold on
figure(1) %sin plan de limpieza
plot(y,xA,'g'); %contaminante pantano A: verde
plot(y,xB,'m'); %contaminante pantano B: magenta
title('Sin plan de limpieza')
xlabel('Tiempo (dias)')
ylabel('Contaminante (T)')
figure(2) %con plan de limpieza
plot(y,Xa,'b') %contaminante pantano A: azul
plot(y,Xb,'r') %contaminante pantano B: rojo
title('Con plan de limpieza')
xlabel('Tiempo (dias)')
ylabel('Contaminante (T)')
hold off
}}

Ambas soluciones quedan reflejadas en las siguientes gráficas:

[[Archivo:Rungekutta.jpg|800px|centro|]]

===Comparación método Euler con método Runge-Kutta===
Para comparar ambos métodos basta con comparar sus gráficas. Realizamos el procedimiento con ambos sistemas.
La gráfica del primer sistema es la siguiente:
[[Archivo:Comparacion1.jpg|600px|centro|]]

Como puede verse la diferencia entre ambos métodos en tan pequeña que las gráficas aparecen superpuestas. Para poder apreciar la diferencia entre ellas ampliamos la grafica cogiendo dos intervalos pequeños de tiempo. Runge-Kutta aparece representado en azul y Euler en verde.
[[Archivo:Comparacion2.jpg|800px|centro|]]

Realizamos el mismo proceso con el segundo sistema y como es de esperar sucede lo mismo:
[[Archivo:Comparacion3.jpg|600px|centro|]]
[[Archivo:Comparacion4.jpg|800px|centro|]]

Una vez comparados ambos métodos se puede afirmar que la diferencia entre ellos es casi inapreciable.

==Contaminante inicial desconocido==

No conocemos la cantidad de contaminante inicial pero sabemos que tras unos días se redujo el contaminante la sólo una tonelada en A y dos en B. Si queremos estimar el contaminante que se vertió inicialmente lo podemos hacer a través del método de Euler con el siguiente código:

{{matlab|codigo=clear all

t0 = 0;
tN = -600;
N = 10000;
h = (tN-t0)/N;
x0 = [1;2];
A = [-0.03 0;0.03 -0.045];
x = x0;
xA(10001) = x(1);
xB(10001) = x(2);
for i=N+1:-1:2
x = x+h*A*x;
if x(2)<=0
xA(1:i-1) = [];
xB(1:i-1) = [];
r = length(xA);
break
else
r = N+1;
end
xA(i-1) = x(1);
xB(i-1) = x(2);
end
t = tN:-h:t0;
tr = length(t);
if tr~=r
t(1:tr-r) = [];
end
hold on
xlabel('Tiempo (días)')
ylabel('Contaminante (T)')
plot(t,xA,'b');
plot(t,xB,'r');
hold off
}}

[[Archivo:SincontaminanteB2.jpg|800px|centro|]]

Obtendremos el vertido inicial en el mismo instante en el que el contaminante en el pantano B se anule. Comopuede verse en la gráfica el contaminante en B nunca se hace cero ya que crece exponencialmente hacia la izquierda hasta un tiempo <math> t=-\infty </math>. Este resultado no es lógico, por lo que modificamos mínimamente los datos que nos dan sobre la cantidad de contaminante en ambos pantanos tras unos días para ver si podemos conseguir un resultado aproximado que sea lógico.

Cambiamos la cantidad de contaminante en B de 2 toneladas a 1.8:

[[Archivo:SincontaminanteB1_8.jpg|800px|centro|]]

Cantidad de contaminante en B=1.999 toneladas:

[[Archivo:SincontaminanteB1_999.jpg|800px|centro|]]

Como puede comprobarse según aumentamos la cantidad de contaminante en B la cantidad de contaminante inicialse va incrementando hasta que, utilizando la condición dada de 2 toneladas de contaminante en B, tiende a <math> t=+\infty </math>. Por lo tanto no se puede concretar ninguna cantidad inicial.

== Autores ==

* Raquel Gómez Vázquez, 187
* Eva Périz Rodríguez, 334
* Daniel Prieto Alonso, 341
* Adrián Salas Calvo, 385

[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:Trabajos 2012-13]]