MateWiki - Contribuciones del usuario [es]

Ecuación del calor en una placa en forma de anillo

2014-05-22T09:49:23Z

Toni Villa: /* Método del trapecio en el tiempo */

{{ TrabajoED | Ecuación del calor en una placa en forma de anillo. Grupo 6B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Elisabet Sánchez López, Ana Santos Martín, Isabel Sastre Furones, Ángel Antonio Villa Figueroa }}

== Introducción ==
Tenemos una placa plana en forma de anillo. Vamos a trabajar en coordenadas polares en el plano y sabemos que la placa está comprendida entre los radios ρ=1 y ρ=6 y su temperatura inicial viene definida por la función u(ρ,0):

<math>\left\{\begin{matrix}100(ρ-1) ……… ρ ∈ [1, 2] \\100 ……… ρ ∈ [2, 5]\\90(6-ρ)+10 ……… ρ ∈ [5, 6]\end{matrix}\right.</math>

Si colocamos en los extremos interior ρ=1 y exterior ρ=6 objetos a una temperatura constante de 0 y 10 grados respectivamente, es decir, las condiciones de frontera serán:

<math>\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.</math>

<math> u_t-\Delta u= u_t- \frac{\partial^2 u }{\partial ρ^2} -\frac{1 }{ρ} \frac{\partial u }{\partial ρ}-\frac{1 }{ρ^2} \frac{\partial^2 u }{\partial θ^2}= u_t- \frac{\partial^2 u }{\partial ρ^2} -\frac{1 }{ρ} \frac{\partial u }{\partial ρ} = 0 </math>

Así la ecuación queda de la siguiente manera:

<math> u_t-\frac{1 }{ρ} u_ρ - u_(ρρ) = 0 </math>

Suponiendo que la temperatura u de la placa únicamente depende de la coordenada radial y del tiempo, es decir, u=u(ρ,t) y sabiendo que satisface la ecuación del calor <math> u_t-\Delta u=0</math> , vamos a plantear el sistema de ecuaciones que satisface u(ρ,t) en coordenadas polares:

== Resolución por el método de diferencias finitas ==
=== Método del trapecio en el tiempo===

Tomamos los pasos de discretización del espacio <math> \Delta ρ=0.1</math> y del tiempo <math> \Delta t=(\Delta ρ/4)</math> en t ∈ [0; 10]

{{matlab|codigo=
%Datos
a=1; b=6; h=0.1;
ya=0; yb=10;
N=(b-a)/h; dt=h/4;
r=[a:h:b]'; rr=[a+h:h:b-h]';
%Definimos K1
K1=diag(ones(1,N-1)*2)+diag(ones(1,N-2)*(-1),-1)+diag(ones(1,N-2)*(-1),1);
K1=K1*(1/h^2);
%Definimos K2
K2=diag(ones(N-2,1)*(-1),1)+diag(ones(N-2,1),-1);
K2=(1/h^2)*K2;
M=diag(1./rr);
K2=M*K2;

T0=0; TF=10;
t=T0

%Definimos la malla
[R,T]=meshgrid(r,t);
%Defino F:
F=zeros(1,N-1)
F(end)=F(end)+10/(h^2)+10/(2*5.9*h); %(donde (b-h)=5.9)
%Condición inicial
for i=1:length(rr)
if rr(i)<2
U(i)=100*(rr(i)-1);
elseif rr(i)<5
U(i)=100;
else
U(i)=90*(6-rr(i))+10;
end
end
sol(1,:)=[0,U,10];
%Método del trapecio
for m=1:length(t)-1
Z=U'+(dt/2)*(-(K1+K2)*U'+F');
W=eye(N-1)+(dt/2)*(K1+K2);
W=W\Z;
sol(m+1,:)=[0,W',10];
U=W'
end
%Resultados gráficos
surf(R,T,sol)
}}

[[Archivo:imagentrapecio.jpg|marco|centro|Superficie temperatura en el eje radial ρ∈[1,6] en el intervalo de tiempo]]

Como se puede observar en la gráfica, la temperatura va disminuyendo hasta alcanzar valores casi constantes a lo largo del tiempo. En el centro (del eje radial) se ve un descenso más brusco, debido a que la temperatura inicial es más alta, y en los extremos se mantiene el mismo valor. Por lo tanto, existe un enfriamiento de la placa, sobre todo en la parte central, hasta unos valores constantes de temperatura: 0 en ρ=1, 10 en ρ=6 y entre 0-10 en ρ∈(1,6).

'''Temperatura en ρ=3'''

Partiendo del programa anterior, queremos ahora dibujar el comportamiento de la temperatura en los puntos para los cuales '''ρ=3''' en una gráfica 2D temperatura/tiempo. Implementamos el siguiente código Matlab y obtenemos dicha gráfica.

{{matlab|codigo=
%Cálculo del punto N que corresponde a ρ=3
N3=(3-a)/h+1;
%Dibujo
plot(sol(:,N3),t)
}}

[[Archivo:imagenP3.jpg|marco|centro|Temperatura en los puntos de la corona ρ=3]]

Con esto, se demuestra lo dicho en el anterior apartado, pero se puede ver de forma mucho más clara para ρ=3 que existe un enfriamiento en la corona, llegando a una temperatura entre los 0 y 10ºC.

[[Archivo:mayorT.jpg|marco|centro|]]
i
[[Archivo:mayorTP3.jpg|marco|centro|Temperatura en los puntos de la corona ρ=3]]
Que demuestra lo comentado anteriormente que cuando se supera el intevalo de tiempo se llegan a valores constantes.

=== Método de Euler explícito en el tiempo ===

Puesto que los datos son los mismos que los usados en el método de diferencias finitas con el método del trapecio en el tiempo, únicamente incluimos en este apartado la implementación del nuevo método utilizado.

{{matlab|codigo=
%Método de Euler explícito
for m=1:length(t)-1
Z=U'+(dt)*(-(K1+K2)*U'+F');
sol(m+1,:)=[0,Z',10];
U=Z';
end
surf(R,T,sol)
}}

[[Archivo:imageneuler.jpg|marco|centro|]]

En esta gráfica resultante se puede observar que, como el método es explícito, el sistema es inestable. Esto significa que el paso tomado es grande, por eso tiene los límites muy elevados.

=== Método de Euler implícito en el tiempo ===

Análogamente, en este apartado se especifica la aplicación del método de Euler implícito, obviando los datos (que continúan siendo los mismos).

{{matlab|codigo=
%Método de Euler implícito:
for m=1:length(t)-1
Z=U'+dt*(F');
W=eye(N-1)+dt*(K1+K2);
W=W\Z;
sol(m+1,:)=[0,W',10]
U=W';
end
}}

[[Archivo:imageneimp.jpg|marco|centro|]]

Obtenemos una superficie muy parecida a la que salía con el método del trapecio.

=== Método de Euler modificado en el tiempo ===

Al utilizar el método de Euler modificado en el tiempo, observamos que éste es inestable, al igual que el método de Euler explícito.

{{matlab|codigo=
%Método de Euler modificado:
for m=1:length(t)-1
M1=-(K1+K2)*U'+F';
M2=-(K1+K2)*(U'+M1*dt)+F';
W=U'+(dt/2)*(M1+M2);
sol(m+1,:)=[0,W',10];
U=W';
end
}}

[[Archivo:imagenemod.jpg|marco|centro|]]

== Estado estacionario ==
=== Con las mismas condiciones ===

El sistema a resolver para tiempos grandes queda de la siguiente manera:

<math>\left\{\begin{matrix} -\frac{1 }{ρ} u_ρ - u_{ρρ} = 0 \\ u(1)=0 ; u(6)=10\end{matrix}\right.</math>

El paso que se toma es de 0.1 y 0.01. Como se trata de un problema de contorno, lo resolvereos implementando el siguiente código:

{{matlab|codigo=
%Datos
a=1; b=6; ya=0; yb=10;
h=[0.1,0.01];
color='rg'
%Para obtener ambas soluciones en la misma gráfica
hold on
for i=1:2
N=(b-a)/h(i); r=[a:h(i):b]';rr=[a+h(i):h(i):b-h(i)]';
K1=1/h(i)^2*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
K2=1/(2*h(i))*( diag(ones(N-2,1),-1)-diag(ones(N-2,1),1) );
M=diag(1./rr);
K2=M*K2;
F=zeros(1,N-1);
F(end)=F(end)+10/h(i)^2+10/(2*5.9*h(i));
u=(K1+K2)\F';
uu=[ya;u;yb];
plot(r,uu,color(i));
end
hold off
}}

[[Archivo:tgrande.jpg|marco|centro|Distribución de la temperatura para tiempos grandes con ambos pasos]]

Como podemos observar, casi no se aprecia la diferencia entre hacerlo con h=0.1 y h=0.01. Además, demuestra que da igual que sean tiempos grandes porque para nuestra placa se obtendrá un valor dentro del intervalo de tiempo dado y, como se había dicho antes, se llegan a valores constantes de temperatura.

[[Archivo:tgrandezoom.jpg|marco|centro|Detalle de la imagen anterior: rojo(0.1) y verde(0.01)]]

En el zoom podemos comprobar como dicha diferencia es mínima.

=== Con las nuevas condiciones ===

Cuando colocamos en la frontera exterior de la placa una pieza aislante, no se produce pérdida de calor en dicho extremo, por lo que el flujo de temperatura en la dirección radial es nulo.

Una vez alcanzado el estado estacionario, el problema que modeliza la temperatura de la placa es:

<math>\left\{\begin{matrix} -\frac{1}{ρ} u_ρ - u_{ρρ} = 0 \\ u(1)=0 ; u_ρ(6)=0\end{matrix}\right.</math>

Para obtener el valor estacionario de la temperatura de a placa, resolvemos el problema aplicando el cambio de variable '''<math> u_ρ = v </math>'''.

Así, tenemos que resolver una ecuación de variables separables cuya solución será: '''<math> u_ρ = v = \frac{c_1}{ρ} </math>'''

Integrando, obtenemos '''<math> u(ρ)=c_1·log ρ + c_2 </math>'''.

Finalmente, obligando al cumplimiento de las condiciones, obtenemos el valor de las constantes:
<math>\left\{\begin{matrix} u(1)=0 => c_2=0 \\ u_ρ(6)=0 => c_1=0\end{matrix}\right.</math>

Por lo tanto, '''<math> u(ρ)=0 </math>'''.

Ahora comprobamos cuánto tarda la temperatura en alcanzar el estado estacionario con un error del 5%. Para ello, utilizamos un método cualquiera de los progamados anteriormente y añadimos al bucle del tiempo las siguientes líneas:

{{matlab|codigo=
for ...
...
if (max(abs(sol)))<0.05
break
end
end
}}

Así, obtenemos como resultado '''t= '''.

== Disco completo ==

Ahora consideraremos que la placa ocupa todo el disco '''ρ<6'''. La solución sólo dependerá de ρ y t, y la condición de frontera será '''<math> u(6,t)=0 </math>'''.

=== Problema de autovalores ===

Para obtener el problema de autovalores asociado '''(PA)''', partimos del problema '''(P)''', cuyas condiciones se imponen teniendo en cuenta la condición de frontera dada y que las soluciones tienen que ser difereciables en '''ρ=0'''. Debido a este último dato, sabemos que en ese punto tiene que existir tangente horizontal (pendiente nula) y obtenemos, así, la condición '''<math> u_ρ(0)=0 </math>'''.

El problema '''(P)''' inicial es:

<math>\left\{\begin{matrix} u_t -\frac{1}{ρ} u_ρ - u_{ρρ} = 0 \\ u_ρ(0)=0 ; u(6)=0\end{matrix}\right.</math>

Aplicando el método de separación de variables, ya que nuestro problema es homogéneo en las condiciones de frontera, buscamos que '''<math> u(ρ,t) </math>''' sea el producto de dos funciones que dependan respectivamente de una de las variables. Así: '''<math> u(ρ,t)=R(ρ)T(t) </math>'''.

Derivando, sustituyendo en '''(P)''' y dividiendo entre '''<math> R(ρ)T(t) </math>''', llegamos al problema de autovalores asociado '''(PA)''':

<math>\left\{\begin{matrix} R''(ρ) +\frac{R'(ρ)}{ρ} = -λR(ρ) \\ R'(0)=0 ; R(6)=0\end{matrix}\right.</math>

=== Relación con las funciones de Bessel ===

La función de Bessel de primera especie, <math> J_0(r) </math>, es la solución de la ecuación diferencial de Bessel:

<math>\left\{\begin{matrix} r^2 φ''(r) + r φ'(r) + r^2 φ(r) = 0 \\ φ(0)=1 ; φ'(0)=0\end{matrix}\right.</math>

Para demostrar que las autofunciones del problema anterior son de la forma <math> φ_k = J_0 (ρλ^{1/2}) </math>, siendo λ el autovalor correspondiente, nos basamos en las similitudes entre a ecuación diferencial de Bessel y nuestro problema de autovalores. <math> R = J_0 (ρλ^{1/2}) </math> tiene que satisfacer '''(PA)'''. Esto se comprueba sin más que derivar aplicando la regla de la cadena y sustituir en el problema de autovalores, lo que nos proporciona una ecuación diferencial de Bessel de la cual '''R''' ya es solución.

Si <math> φ_k = J_0 (ρλ^{1/2}) </math> son las autofunciones del problema que nos ocupa, siendo λ los autovalores, los ceros de la función de Bessel se relacionan con los autovalores de la siguiente forma:
<math> φ_k (6)=0 \\ J_0 (6·λ^{1/2})=0 </math>

Esto sucede si <math> 6·λ^{1/2} </math> es un cero de la función de Bessel. Por tanto, si llamamos '''<math> r_0 </math>''' a las raíces de la función de Bessel, obtenemos la relación:
<math> r_0 = 6·λ^{1/2} </math>

=== Autovalores y autofunciones ===

'''Autovalores'''

Una vez demostrado que las autofunciones de '''(P)''' son de la forma <math> φ_k = J_0 (ρλ^{1/2}) </math>, para calcular los cinco primeros autovalores utilizaremos la relación anterior e implementaremos un pequeño cálculo vectorial en Matlab.

Tomamos como condición <math> φ_k (6)=0 </math>, obtenemos la relación <math> r_0 = 6·λ^{1/2} </math> y utilizamos los valores de los cinco primeros ceros de la función de Bessel dados. Así, tenemos:

{{matlab|codigo=
%Cálculo vectorial para la obtención de los 5 primeros autovalores de (P)
r0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=r0.^2/36 %Vector de autovalores
}}

Los resultados que obtenemos son:
'''<math> λ_1 = 0.1606, λ_2=0.8464, λ_3=2.0802, λ_4=3.8622, λ_5=6.1925 </math>'''

'''Autofunciones'''

Para obtener la gráfica con las primeras cinco autofunciones, utilizaremos la función '''besselj''' de Matlab. '''J=besselj(nu,Z)''' computa la función de Bessel de orden nu (en nuestro caso, nu=0) para cada elemento del vector Z (en nuestro caso, el vector de <math> ρλ^{1/2} </math>). Programando en Matlab, obtenemos los siguientes resultados:

{{matlab|codigo=
%Datos del problema
clear all
a=0; b=6;
r0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=r0.^2/36; %Vector de autovalores
K=5; %Número de autovalores y autofunciones buscadas
%Datos de la discretización
h=0.1;
N=(b-a)/h;
%Vector de radios
p=a:h:b;
%Matriz de ceros de la función de Bessel y matriz de autofunciones
J=zeros(N,K);
z=zeros(N,K);
%Inicialización de las matrices
j0=besselj(0,0);
J(:,1)=j0;
z0=0;
z(:,1)=z0;
for k=1:K
for i=1:N
z(i,k)=(p(1,i))'*sqrt(A(k));
J(i,k)=besselj(0,z(i,k));
end
end
%Representación gráfica de las autofunciones J
plot (J)
}}

[[Archivo:AutofuncionesGraph.jpg]]

=== Resolución por el método de Fourier ===

Por último, aproximaremos la solución de la ecuación del calor utilizando cinco términos del desarrollo de Fourier cuando el dato inicial es <math> u(ρ,0)=-log \frac{|ρ+0.1|}{6.1} </math>. Ensayamos con soluciones de la forma:: <math> u(ρ,t)= T_k(t)φ_k(ρ)\\ k=1,2,3,4,5 </math>

Debe satisfacer la condición inicial del problema:: <math> u(ρ,0)= T_k(0)J_0(ρλ_k^{1/2}) =-log \frac{|ρ+0.1|}{6.1} </math>

De la serie de Fourier de <math> h(ρ)=-log \frac{|ρ+0.1|}{6.1} </math>, obtenemos la expresión de los coeficientes de Fourier que se utilizará para programar en Matlab. En ella, tendremos en cuenta que, al tratarse de funciones de Bessel, las autofuciones son ortogonales respecto al producto escalar.

El programa de Matlab para obtener la solución de la ecuación del calor por el método de Fourier es:

{{matlab|codigo=
a=0; b=6;
r0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=r0.^2/36; %Vector de autovalores
K=5; %Número de autovalores y autofunciones buscadas
%Datos de la discretización
h=0.1;
N=(b-a)/h;
%Vector de radios
p=a:h:b;
%Matriz de ceros de la función de Bessel y matriz de autofunciones
J=zeros(N,K);
z=zeros(N,K);
%Inicialización de las matrices
j0=besselj(0,0);
J(:,1)=j0;
z0=0;
z(:,1)=z0;
%Matriz de soluciones
U=0;
%Mallado en el tiempo
t0=0; tf=2;
m=0.1; %Longitud de paso
t=t0:m:tf;
[Mp,Mt]=meshgrid(p,t);
%Bucle
for k=1:K
for i=1:N
z(i,k)=(p(1,i))'*sqrt(A(k));
J(i,k)=besselj(0,z(i,k));
h=-log((abs(p(1,i)+0.1))/6.1); %Vector de condiciones iniciales
ck=trapz(p,(p').*h.*J(i,k))/trapz(p,(p').*J(i,k).*J(i,k)); %Coeficientes de Fourier
U=U+ck*exp(-A(k)*Mt).*J(i,k);
end
end
%Representación gráfica de la solución obtenida
mesh(Mp,Mt,U)
}}

La gráfica obtenida es la siguiente:

[[Archivo:FourierGraph.jpg|centro]]

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

Archivo:MayorTP3.jpg

2014-05-22T09:45:34Z

Toni Villa:

Archivo:MayorT.jpg

2014-05-22T09:45:13Z

Toni Villa:

Ecuación del calor en una placa en forma de anillo

2014-05-19T19:48:59Z

Toni Villa: /* Tiempos grandes (estado estacionario) */

Ecuación del calor en una placa en forma de anillo

2014-05-19T19:36:18Z

Toni Villa: /* Euler explícito, implícito y modificado */

Ecuación del calor en una placa en forma de anillo

2014-05-19T19:33:42Z

Toni Villa: /* Temperatura en ρ=3 */

Ecuación del calor en una placa en forma de anillo

2014-05-19T19:32:11Z

Toni Villa: /* Temperatura en ρ=3 */

Ecuación del calor en una placa en forma de anillo

2014-05-19T19:29:20Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:57:19Z

Toni Villa: /* Tiempos grandes (estado estacionario) */

Archivo:Tgrandezoom.jpg

2014-05-18T22:51:40Z

Toni Villa:

Archivo:Tgrande.jpg

2014-05-18T22:50:33Z

Toni Villa:

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:43:45Z

Toni Villa: /* Tiempos grandes (estado estacionario) */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:21:34Z

Toni Villa: /* Tiempos grandes (estado estacionario) */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:20:22Z

Toni Villa: /* Tiempos grandes (estado estacionario) */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:18:03Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:17:27Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:16:27Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:15:24Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:14:53Z

Toni Villa:

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:09:03Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T22:08:17Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T21:54:15Z

Toni Villa: /* Euler explícito, implícito y modificado */

Archivo:Imagenemod.jpg

2014-05-18T21:43:00Z

Toni Villa:

Archivo:Imageneimp.jpg

2014-05-18T21:42:40Z

Toni Villa:

Archivo:Imageneuler.jpg

2014-05-18T21:42:01Z

Toni Villa:

Ecuación del calor en una placa en forma de anillo

2014-05-18T21:16:05Z

Toni Villa: /* Euler modificado */

Ecuación del calor en una placa en forma de anillo

2014-05-18T21:09:04Z

Toni Villa: /* Temperatura en ρ=3 */

Archivo:ImagenP3.jpg

2014-05-18T21:05:43Z

Toni Villa:

Ecuación del calor en una placa en forma de anillo

2014-05-18T21:05:15Z

Toni Villa: /* Temperatura en ρ=3 */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:56:52Z

Toni Villa: /* Euler modificado */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:56:31Z

Toni Villa: /* Euler implícito */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:54:48Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:54:29Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:53:13Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:52:45Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:44:19Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Archivo:Imagentrapecio.jpg

2014-05-18T20:41:35Z

Toni Villa:

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:31:31Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:31:06Z

Toni Villa: /* Método de diferencias finitas con el método del Trapecio */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:30:09Z

Toni Villa: /* Método de diferencias finitas */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:23:44Z

Toni Villa: /* Método de diferencias finitas */

Ecuación del calor en una placa en forma de anillo

2014-05-18T20:20:58Z

Toni Villa: /* Método de diferencias finitas */

Ecuación del calor en una placa en forma de anillo

2014-05-18T19:55:52Z

Toni Villa: /* Introducción */

Ecuación del calor en una placa en forma de anillo

2014-05-18T18:16:05Z

Toni Villa: /* 2.Método de diferencias finitas */

Ecuación del calor en una placa en forma de anillo

2014-05-18T17:31:32Z

Toni Villa: /* Disco */

Ecuación del calor en una placa en forma de anillo

2014-05-18T17:12:31Z

Toni Villa: /* 2.Método de diferencias finitas */

Ecuación del calor en una placa en forma de anillo

2014-05-18T17:06:37Z

Toni Villa: /* Anillo circular */

Modelo predador-presa de A. Lotka y V. Volterra

2014-03-04T10:21:37Z

Toni Villa: /* Interpretación de resultados */

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Grupo 6B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Paula De Santos Muñoz, María del Mar García Reinaldos, Elisabet Sánchez López, Ana Santos Martín, Isabel Sastre Furones, Ángel Antonio Villa Figueroa }}
El modelo predador-presa fue propuesto de forma independiente a comienzos del siglo pasado por A. Lotka y V. Volterra. Es un modelo mátematico para la dinámica de poblaciones de especies competidoras (predadores y presas), que consiste en un sistema de ecuaciones diferenciales de primer orden no lineales correspondientes al cambio de las poblaciones de presa y predador por unidad de tiempo.
== Interpretación del problema de valor inicial ==

El problema planteado trabaja la existencia de dos especies de presas distintas que comparten una población de depredadores. Llamaremos:
* x1: especie presa 1
* x2: especie presa 2
* x3: población depredadora

Nuestro trabajo se basa en un estudio del crecimiento de las presas y del depredador en el tiempo, observando las interacciones entre ellas y teniendo en cuenta que:
# En condiciones iniciales, x1 cuenta con una población de p0, x2 con una población de q0 y finalmente x3 con una población d0.
# Las poblaciones evolucionan interactuando entre ellas, sin tener en cuenta agentes externos.
# Las presas cuenta con una cantidad suficiente de comida para ser alimentada.
# La alimentación de los depredadores depende completamente de la presa.

Con estos datos podemos definir las ecuaciones de nuestro sistema. La primera ecuación nos define la evolución de la presa en el tiempo cuando está interactuando con la especie depredadora donde:
* A1x1 es la tasa de natalidad de la presa 1, siguiendo una ley exponencial.
* A2x1x3 es la tasa de mortalidad de la presa 1 afectada por los depredadores.

:<math>\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\</math>

De forma análoga a la primera ecuación, podemos definir la segunda.

:<math>\frac{\mathrm{d} x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\</math>

En la tercera ecuación definimos la evolución de x3, definiendo la constante C1 multiplicada por (-1), lo que nos hace ver
la extinción de los depredadores ante la inexistencia de presas y el crecimiento de los mismos gracias a la alimentación aportada por las presas.

:<math>\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\</math>

Con estás tres ecuaciones queda definido el sistema de ecuaciones de Lotka y Volterra a resolver como un problema de valor inicial.

:<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>

== Resolución por el método de Euler ==
=== Programación del método numérico ===
El programa utilizado para resolver el problema de valor inicial mediante el método de Euler en un intervalo de tiempo '''tє[0,300]''', para una longitud de paso '''h=0.1''', constantes '''A1=0.35''', '''A2=0.6''', '''B1=0.3''', '''B2=0.5''', '''C1=0.37''', '''C2=0.04''' y '''C3=0.035''', y valores iniciales de '''p0=0.8''' millones de presas de un tipo, '''q0=2.4''' millones de presas de otro tipo y '''d0=0.2''' millones de predadores es el siguiente:
{{matlab|codigo=
%Datos del problema
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;
%Datos de la discretización
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
p0=0.8;
q0=2.4;
d0=0.2;
y0=[p0 q0 d0]';
%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 y matriz de soluciones aproximadas
t=t0:h:tf;
y=zeros(3,N+1);
%Inicialización
y(:,1)=y0;
yy=y0;
for n=1:N
yy=yy+h*(A*yy+B*yy(1)*yy(3)+C*yy(2)*yy(3));
y(:,n+1)=yy;
end
figure(1)
plot(t,y)
title('Evolución de las poblaciones de presas (azul y verde) y predador (rojo) en función del tiempo')
xlabel('Tiempo t (años)')
ylabel('Millones de ejemplares')
axis tight
figure(2)
subplot(1,2,1)
plot(y(1,:),y(3,:))
title('Población del primer tipo de presa en función de la de predadores')
xlabel('Presas tipo 1 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
subplot(1,2,2)
plot(y(2,:),y(3,:))
title('Población del segundo tipo de presa en función de la de predadores')
xlabel('Presas tipo 2 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
}}
=== Gráficas ===
==== Considerando la evolución de las especies en el tiempo ====
De la aplicación del método numérico se obtienen las siguientes gráficas, de las que sacaremos conclusiones posteriormente.

[[Archivo:GraphEuler.jpg|centro]]

En esta primera gráfica, obtenida como resultado de utilizar como longitud de paso '''h=0.1''', se muestran las tres funciones equivalentes a la evolución en el tiempo de los dos tipos de presas y del predador. Cada curva presenta unos máximos periódicos desplazados los unos respecto de los otros, lo cual se debe al desfase con el que evolucionan las distintas especies: primero aumentan las poblaciones de ambas presas (creciendo más aquella de la que inicialmente hay más ejemplares) y después comienza a aumentar la población de la especie predadora, pues puede alimentarse fácilmente. Ésta última alcanzará su primer máximo cuando las presas disminuyen, tras lo cual disminuye también el número de predadores al escasear el alimento. Este proceso se dará periódicamente a lo largo de los 300 años. Es notable la diferencia en la evolución de ambos tipos de presas: mientras que uno tiende a aumentar, el otro disminuye.

Si ahora utilizamos como longitud de paso '''h=1''', obtenemos esta segunda gráfica de la cual podemos concluir que, debido precisamente al aumento de ''h'', el método de Euler pierde su validez.

[[Archivo:GraphEuler1.jpg|centro]]

==== Considerando la interacción entre especies ====
Dado que hemos comprobado anteriormente que el método de Euler con longitud de paso h=1 no es efectivo, evaluaremos la interacción entre especies utilizando '''h=0.1''' como longitud de paso. Las gráficas resultantes son las siguientes:

[[Archivo:GraphEuler4.jpg|centro]]

=== Interpretación de resultados ===
Al aumentar la longitud paso de '''h=0.1''' a '''h=1''', lo que observamos es una gran pérdida de precisión en los resultados, ya que el valor de la aproximación calculada por el método de Euler se aleja demasiado de la solución exacta del problema de valor inicial.
Esto se debe a que este método numérico se basa en aproximaciones mediante la tangente a una determinada función. En nuestro caso, las diferencias entre los máximos y mínimos de las funciones que describen la evolución de las especies son considerables (son funciones oscilatorias) y se producen en intervalos de tiempo ''cortos'', por lo que, al utilizar una longitud de paso ''grande'' (como puede ser ''h=1''), la inexactitud de la aproximación calculada resulta inaceptable. Por el contrario, si disminuimos la longitud de paso, el valor de la aproximación se considera válido. Esto se explica porque el método de Euler alcanza mayor nivel de convergencia cuanto menor es el valor de ''h''; sin embargo, hay que tener en cuenta que, si disminuimos mucho este valor, el programa puede llegar a fallar.

== Resolución por el método de Euler modificado ==
=== Programación del método numérico ===
A continuación, se resuelve el problema de valor inicial mediante el método de Euler modificado. Como datos tenemos una longitud de paso '''h=0.1''', constantes '''A1=0.35''', '''A2=0.6''', '''B1=0.3''', '''B2=0.5''', '''C1=0.37''', '''C2=0.04''' y '''C3=0.035''', y valores iniciales de '''p0=2''' millones de presas de un tipo, '''q0=1.4''' millones de presas de otro tipo y '''d0=1''' millones de predadores. Primero tomaremos un intervalo de tiempo '''tє[0,100]'''. El programa se describe a continuación:
{{matlab|codigo=
%Datos del problema
t0=0;
tf=100;
A1=0.35; A2=0.6;
B1=0.3; B2=0.5;
C1=0.37; C2=0.04; C3=0.035;
%Datos de la discretización
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
p0=2;
q0=1.4;
d0=1;
y0=[p0 q0 d0]';
%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;
%Inicialización
y=y0;
for n=1:N
k1=(A*y+y(1)*y(3)*B+y(2)*y(3)*C);
k2=(A*(y+h*k1)+(y(1)+h*k1(1))*(y(3)+h*k1(3))*B+(y(2)+h*k1(2))*(y(3)+h*k1(3))*C);
y=y+h*(k1+k2)/2;
y1(n+1)=y(1);
y2(n+1)=y(2);
y3(n+1)=y(3);
end
%Gráfica poblaciones-tiempo
figure(1)
hold on
plot(t,y1,'b')
plot(t,y2,'g')
plot(t,y3,'r')
title('Evolución de las poblaciones de presas (azul y verde) y predador (rojo) en función del tiempo')
xlabel('Tiempo t (años)')
ylabel('Millones de ejemplares')
axis tight
hold off
%Gráficas presa-predador
figure(2)
subplot(1,2,1)
plot(y1,y3,'b')
title('Población del primer tipo de presa en función de la de predadores')
xlabel('Presas tipo 1 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
subplot(1,2,2)
plot(y2,y3,'r')
title('Población del segundo tipo de presa en función de la de predadores')
xlabel('Presas tipo 2 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
}}
Hacemos lo mismo para el segundo intervalo de tiempo '''tє[0,300]'''. El único cambio que sufre el programa es '''tf=300'''.

=== Gráficas ===
==== Considerando la evolución de las especies en el tiempo ====
De la aplicación del método numérico para un tiempo '''tє[0,100]''', tenemos la siguiente gráfica:

[[Archivo:GraphMEuler1.jpg|centro]]

De ella podemos concluir que, durante los primeros 100 años, la evolución de ambos tipos de presas sigue direcciones opuestas: mientras la curva que describe la evolución de la presa de la que hay más ejemplares inicialmente tiende a que sus máximos sean cada vez más pequeños, los máximos de la otra son cada vez mayores. Se observa que, al finalizar el período, el número de ejemplares de ambas presas es similar. Por otro lado, la evolución de la cantidad de predadores en el tiempo es similar a la obtenida al aplicar el método de Euler, aumentando el número de ejemplares tras producirse un pico en las cantidades de presas. Cabe destacar, sin embargo, que con el método de Euler modificado obtenemos gráficas más suaves, lo cual permite aumentar un poco la longitud de paso ''h'' sin acumular demasiado error.

Si ahora aumentamos nuestro intervalo de tiempo hasta '''tє[0,300]''', podemos comprobar en la gráfica siguiente que la evolución del número de presas sigue la tendencia indicada durante los primeros 100 años: la especie predadora mantiene su evolución oscilatoria con máximos constantes, las presas cuyos máximos disminuían progresivamente casi se han extinguido transcurridos los 300 años y las otras alcanzan sus cifras máximas al final del período. Esto se debe a que los predadores se centran en un determinado tipo de presa, con lo cual la otra se encuentra menos amenazada.

[[Archivo:GraphMEuler2.jpg|centro]]

==== Considerando la interacción entre especies ====
Para evaluar la interacción entre especies tendremos en cuenta las gráficas obtenidas para un tiempo '''tє[0,300]''', ya que se observa mejor la variación de las poblaciones de presas respecto de la de predadores.

[[Archivo:GraphMEuler5.jpg|centro]]

Estas gráficas nos permiten concluir que los puntos en los que el ecosistema sería estable son los orígenes (0,0), los cuales despreciamos (pues si lo interpretamos en términos de poblaciones, no hay especies interactuando), y los puntos donde se encuentran los centros de las espirales. Como podemos comprobar, ninguna de las dos curvas converge al mismo centro, así que no podemos asegurar que los ecosistemas sean estables.

== Resolución por el método trapezoidal ==
El método trapezoidal es un método implícito que, en nuestro caso, requiere la resolución de una ecuación no lineal. Al no ser posible su resolución, se descarta el método y no se tendrá en cuenta a la hora de comparar resultados.

== Resolución por el método de Runge-Kutta ==
=== Programación del método numérico ===
A continuación, se resuelve el problema de valor inicial mediante el método de Runge- Kutta de cuarto orden. Como datos tenemos una longitud de paso '''h=0.1''', constantes '''A1=0.4''', '''A2=0.7''', '''B1=0.2''', '''B2=0.4''', '''C1=0.37''', '''C2=0.04''' y '''C3=0.03''', y valores iniciales de '''p0=3.5''' millones de presas de un tipo, '''q0=0.2''' millones de presas de otro tipo y '''d0=0.4''' millones de predadores. En este caso tomaremos un intervalo de tiempo '''tє[0,500]'''. El programa de obtención del método numérico se especifica a continuación:
{{matlab|codigo=
t0=0; tN=500; %Extremos del intervalo
y0=[3.5;0.2;0.4]; %Valores de población iniciales
y=y0;
P1(1)=y(1); %Valor inicial de P1=presas tipo I
P2(1)=y(1); %Valor inicial de P2=presas tipo II
D(1)=y(1); %Valor inicial de D=depredadores
h=0.1; %Amplitud de los subintervalos
N=(tN-t0)/h; %Número de subintervalos
A=[ 0.4 0 0;0 0.2 0;0 0 -0.37]; %Parámetros
%Algoritmo
for k=1:N
k1=(A*y+y(1)*y(3)*[-0.7;0;0.04]+y(2)*y(3)*[0;-0.4;0.03]);
k2=(A*(y+(h/2)*k1)+(y(1)+(h/2)*k1(1))*(y(3)+(h/2)*k1(3))*[-0.7;0;0.04]+...
(y(2)+(h/2)*k1(2))*(y(3)+(h/2)*k1(3))*[0;-0.4;0.03]);
k3=(A*(y+(h/2)*k2)+(y(1)+(h/2)*k2(1))*(y(3)+(h/2)*k2(3))*[-0.7;0;0.04]+...
(y(2)+(h/2)*k2(2))*(y(3)+(h/2)*k2(3))*[0;-0.4;0.03]);
k4=(A*(y+h*k3)+(y(1)+h*k3(1))*(y(3)+h*k3(3))*[-0.7;0;0.04]+...
(y(2)+h*k3(2))*(y(3)+h*k3(3))*[0;-0.4;0.03]);
y=y+h/6*(k1+2*k2+2*k3+k4);
P1(k+1)=y(1);
P2(k+1)=y(2);
D(k+1)=y(3);
end

t=t0:h:tN;

%CÁLCULO DE INSTANTES MÁXIMOS Y MÍNIMOS

[maximopresaI,x1]= max(P1)
t1=t(x1)
[minimopresaI,x2]= min(P1)
t2=t(x2)
[maximopresaII,x3]= max(P2)
t3=t(x3)
[minimopresaII,x4]= min(P2)
t4=t(x4)
[maximopredadores,x5]= max(D)
t5=t(x5)
[minimopredadores,x6]= min(D)
t6=t(x6)

%GRÁFICA PRESA TIPO I/PRESA TIPO II/DEPREDADORES

hold on
figure(1)
plot(t,P1,'c')
plot(t,P2,'r')
plot(t,D,'g')
title(' Población de presas (tipos I y II) y depredadores')
xlabel('t(años)')
ylabel('población(millones de individuos)')
hold off
}}

=== Gráficas: Evolución de las especies en el tiempo ===
[[Archivo:GraficaRungeKutta.JPG|centro]]

=== Interpretación de resultados ===

En este gráfico se observan las variaciones poblacionales a lo largo del tiempo de las tres especies; la presa tipo I en color cian, la presa de tipo II en color rojo y los predadores se muestran en verde. Interpretamos que la población de predadores empieza a aumentar cuando la población de presas es máxima, es decir cuando tienen mayor cantidad de comida. Este crecimiento se mantiene hasta que el nivel de presas alcanza un punto crítico, donde deja de haber suficientes presas para sostener la población de predadores, lo que genera una disminución en el número de dichos depredadores, que a su vez da lugar a un aumento en el número de presas.

Podemos observar dos comportamientos distintos según el tipo de presa al que nos referimos:
:* En la presa de tipo II observamos a simple vista que la población inicial es muy reducida, incluso menor que la de depredadores. Esto provoca que a pesar de que este tipo de presa no sea el único del que se alimenta el depredador, esta especie termina por extinguirse, como bien se aprecia en la gráfica.
:*En la presa de tipo I, vemos que a parte de que su población inicial es muchísimo mayor que la de depredadores, su velocidad de crecimiento también es mayor, y esto asegura la supervivencia de ambas especies.
Como resultado conluimos que tanto la cantidad de depredadores como la de presas del tipo que no se han extinguido mantienen sus ciclos de crecemiento y decrecimiento por lo que ambas se mantienen constantes en lo que se refiere a máximos y mínimos.

Dichos máximos y mínimos se adjuntan en la siguiente tabla, con los correspondientes instantes en los que tienen lugar:

[[Archivo:Maximosminimos.JPG|centro]]

== Conclusiones ==
Tras analizar los resultados obtenidos aplicando los diferentes métodos numéricos, llegamos a las siguientes conclusiones:
* El método de Euler es un método de orden 1, por lo que nos ofrecerá peores resultados que los métodos de órdenes superiores, como son el método de Euler modificado (Runge-Kutta de orden 2) y el Runge-Kutta de orden 4. Esto se comprueba al comparar los distintos resultados obtenidos.
* En comparación con los otros métodos, es preferible el '''método de Runge-Kutta''', pues aporta soluciones más precisas. Esto se debe a que permite buenas aproximaciones aumentando el valor de la longitud de paso ''h'', ya que no se acumula el error del cálculo mediante tangentes que se desvían de la solución exacta.

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

Modelo predador-presa de A. Lotka y V. Volterra

2014-03-04T10:20:53Z

Toni Villa: /* Interpretación de resultados */

{{ Beta }}

{{ TrabajoED | Modelo predador-presa. Grupo 6B | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED13/14|Curso 2013-14]] | Paula De Santos Muñoz, María del Mar García Reinaldos, Elisabet Sánchez López, Ana Santos Martín, Isabel Sastre Furones, Ángel Antonio Villa Figueroa }}
El modelo predador-presa fue propuesto de forma independiente a comienzos del siglo pasado por A. Lotka y V. Volterra. Es un modelo mátematico para la dinámica de poblaciones de especies competidoras (predadores y presas), que consiste en un sistema de ecuaciones diferenciales de primer orden no lineales correspondientes al cambio de las poblaciones de presa y predador por unidad de tiempo.
== Interpretación del problema de valor inicial ==

El problema planteado trabaja la existencia de dos especies de presas distintas que comparten una población de depredadores. Llamaremos:
* x1: especie presa 1
* x2: especie presa 2
* x3: población depredadora

Nuestro trabajo se basa en un estudio del crecimiento de las presas y del depredador en el tiempo, observando las interacciones entre ellas y teniendo en cuenta que:
# En condiciones iniciales, x1 cuenta con una población de p0, x2 con una población de q0 y finalmente x3 con una población d0.
# Las poblaciones evolucionan interactuando entre ellas, sin tener en cuenta agentes externos.
# Las presas cuenta con una cantidad suficiente de comida para ser alimentada.
# La alimentación de los depredadores depende completamente de la presa.

Con estos datos podemos definir las ecuaciones de nuestro sistema. La primera ecuación nos define la evolución de la presa en el tiempo cuando está interactuando con la especie depredadora donde:
* A1x1 es la tasa de natalidad de la presa 1, siguiendo una ley exponencial.
* A2x1x3 es la tasa de mortalidad de la presa 1 afectada por los depredadores.

:<math>\frac{\mathrm{d} x_1}{\mathrm{d} t}=A_1x_1-A_2x_1x_3\\</math>

De forma análoga a la primera ecuación, podemos definir la segunda.

:<math>\frac{\mathrm{d} x_2}{\mathrm{d} t}=B_1x_2-B_2x_2x_3\\</math>

En la tercera ecuación definimos la evolución de x3, definiendo la constante C1 multiplicada por (-1), lo que nos hace ver
la extinción de los depredadores ante la inexistencia de presas y el crecimiento de los mismos gracias a la alimentación aportada por las presas.

:<math>\frac{\mathrm{d}x_3}{\mathrm{d} t}=-C_1x_3+C_2x_1x_3+C_3x_2x_3\\</math>

Con estás tres ecuaciones queda definido el sistema de ecuaciones de Lotka y Volterra a resolver como un problema de valor inicial.

:<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>

== Resolución por el método de Euler ==
=== Programación del método numérico ===
El programa utilizado para resolver el problema de valor inicial mediante el método de Euler en un intervalo de tiempo '''tє[0,300]''', para una longitud de paso '''h=0.1''', constantes '''A1=0.35''', '''A2=0.6''', '''B1=0.3''', '''B2=0.5''', '''C1=0.37''', '''C2=0.04''' y '''C3=0.035''', y valores iniciales de '''p0=0.8''' millones de presas de un tipo, '''q0=2.4''' millones de presas de otro tipo y '''d0=0.2''' millones de predadores es el siguiente:
{{matlab|codigo=
%Datos del problema
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;
%Datos de la discretización
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
p0=0.8;
q0=2.4;
d0=0.2;
y0=[p0 q0 d0]';
%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 y matriz de soluciones aproximadas
t=t0:h:tf;
y=zeros(3,N+1);
%Inicialización
y(:,1)=y0;
yy=y0;
for n=1:N
yy=yy+h*(A*yy+B*yy(1)*yy(3)+C*yy(2)*yy(3));
y(:,n+1)=yy;
end
figure(1)
plot(t,y)
title('Evolución de las poblaciones de presas (azul y verde) y predador (rojo) en función del tiempo')
xlabel('Tiempo t (años)')
ylabel('Millones de ejemplares')
axis tight
figure(2)
subplot(1,2,1)
plot(y(1,:),y(3,:))
title('Población del primer tipo de presa en función de la de predadores')
xlabel('Presas tipo 1 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
subplot(1,2,2)
plot(y(2,:),y(3,:))
title('Población del segundo tipo de presa en función de la de predadores')
xlabel('Presas tipo 2 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
}}
=== Gráficas ===
==== Considerando la evolución de las especies en el tiempo ====
De la aplicación del método numérico se obtienen las siguientes gráficas, de las que sacaremos conclusiones posteriormente.

[[Archivo:GraphEuler.jpg|centro]]

En esta primera gráfica, obtenida como resultado de utilizar como longitud de paso '''h=0.1''', se muestran las tres funciones equivalentes a la evolución en el tiempo de los dos tipos de presas y del predador. Cada curva presenta unos máximos periódicos desplazados los unos respecto de los otros, lo cual se debe al desfase con el que evolucionan las distintas especies: primero aumentan las poblaciones de ambas presas (creciendo más aquella de la que inicialmente hay más ejemplares) y después comienza a aumentar la población de la especie predadora, pues puede alimentarse fácilmente. Ésta última alcanzará su primer máximo cuando las presas disminuyen, tras lo cual disminuye también el número de predadores al escasear el alimento. Este proceso se dará periódicamente a lo largo de los 300 años. Es notable la diferencia en la evolución de ambos tipos de presas: mientras que uno tiende a aumentar, el otro disminuye.

Si ahora utilizamos como longitud de paso '''h=1''', obtenemos esta segunda gráfica de la cual podemos concluir que, debido precisamente al aumento de ''h'', el método de Euler pierde su validez.

[[Archivo:GraphEuler1.jpg|centro]]

==== Considerando la interacción entre especies ====
Dado que hemos comprobado anteriormente que el método de Euler con longitud de paso h=1 no es efectivo, evaluaremos la interacción entre especies utilizando '''h=0.1''' como longitud de paso. Las gráficas resultantes son las siguientes:

[[Archivo:GraphEuler4.jpg|centro]]

=== Interpretación de resultados ===
Al aumentar la longitud paso de '''h=0.1''' a '''h=1''', lo que observamos es una gran pérdida de precisión en los resultados, ya que el valor de la aproximación calculada por el método de Euler se aleja demasiado de la solución exacta del problema de valor inicial.
Esto se debe a que este método numérico se basa en aproximaciones mediante la tangente a una determinada función. En nuestro caso, las diferencias entre los máximos y mínimos de las funciones que describen la evolución de las especies son considerables (son funciones oscilatorias) y se producen en intervalos de tiempo ''cortos'', por lo que, al utilizar una longitud de paso ''grande'' (como puede ser ''h=1''), la inexactitud de la aproximación calculada resulta inaceptable. Por el contrario, si disminuimos la longitud de paso, el valor de la aproximación se considera válido. Esto se explica porque el método de Euler alcanza mayor nivel de convergencia cuanto menor es el valor de ''h''; sin embargo, hay que tener en cuenta que, si disminuimos mucho este valor, el programa puede llegar a fallar.

== Resolución por el método de Euler modificado ==
=== Programación del método numérico ===
A continuación, se resuelve el problema de valor inicial mediante el método de Euler modificado. Como datos tenemos una longitud de paso '''h=0.1''', constantes '''A1=0.35''', '''A2=0.6''', '''B1=0.3''', '''B2=0.5''', '''C1=0.37''', '''C2=0.04''' y '''C3=0.035''', y valores iniciales de '''p0=2''' millones de presas de un tipo, '''q0=1.4''' millones de presas de otro tipo y '''d0=1''' millones de predadores. Primero tomaremos un intervalo de tiempo '''tє[0,100]'''. El programa se describe a continuación:
{{matlab|codigo=
%Datos del problema
t0=0;
tf=100;
A1=0.35; A2=0.6;
B1=0.3; B2=0.5;
C1=0.37; C2=0.04; C3=0.035;
%Datos de la discretización
h=0.1;
N=(tf-t0)/h;
%Condiciones iniciales
p0=2;
q0=1.4;
d0=1;
y0=[p0 q0 d0]';
%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;
%Inicialización
y=y0;
for n=1:N
k1=(A*y+y(1)*y(3)*B+y(2)*y(3)*C);
k2=(A*(y+h*k1)+(y(1)+h*k1(1))*(y(3)+h*k1(3))*B+(y(2)+h*k1(2))*(y(3)+h*k1(3))*C);
y=y+h*(k1+k2)/2;
y1(n+1)=y(1);
y2(n+1)=y(2);
y3(n+1)=y(3);
end
%Gráfica poblaciones-tiempo
figure(1)
hold on
plot(t,y1,'b')
plot(t,y2,'g')
plot(t,y3,'r')
title('Evolución de las poblaciones de presas (azul y verde) y predador (rojo) en función del tiempo')
xlabel('Tiempo t (años)')
ylabel('Millones de ejemplares')
axis tight
hold off
%Gráficas presa-predador
figure(2)
subplot(1,2,1)
plot(y1,y3,'b')
title('Población del primer tipo de presa en función de la de predadores')
xlabel('Presas tipo 1 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
subplot(1,2,2)
plot(y2,y3,'r')
title('Población del segundo tipo de presa en función de la de predadores')
xlabel('Presas tipo 2 (millones de ejemplares)')
ylabel('Predadores (millones de ejemplares)')
axis tight
}}
Hacemos lo mismo para el segundo intervalo de tiempo '''tє[0,300]'''. El único cambio que sufre el programa es '''tf=300'''.

=== Gráficas ===
==== Considerando la evolución de las especies en el tiempo ====
De la aplicación del método numérico para un tiempo '''tє[0,100]''', tenemos la siguiente gráfica:

[[Archivo:GraphMEuler1.jpg|centro]]

De ella podemos concluir que, durante los primeros 100 años, la evolución de ambos tipos de presas sigue direcciones opuestas: mientras la curva que describe la evolución de la presa de la que hay más ejemplares inicialmente tiende a que sus máximos sean cada vez más pequeños, los máximos de la otra son cada vez mayores. Se observa que, al finalizar el período, el número de ejemplares de ambas presas es similar. Por otro lado, la evolución de la cantidad de predadores en el tiempo es similar a la obtenida al aplicar el método de Euler, aumentando el número de ejemplares tras producirse un pico en las cantidades de presas. Cabe destacar, sin embargo, que con el método de Euler modificado obtenemos gráficas más suaves, lo cual permite aumentar un poco la longitud de paso ''h'' sin acumular demasiado error.

Si ahora aumentamos nuestro intervalo de tiempo hasta '''tє[0,300]''', podemos comprobar en la gráfica siguiente que la evolución del número de presas sigue la tendencia indicada durante los primeros 100 años: la especie predadora mantiene su evolución oscilatoria con máximos constantes, las presas cuyos máximos disminuían progresivamente casi se han extinguido transcurridos los 300 años y las otras alcanzan sus cifras máximas al final del período. Esto se debe a que los predadores se centran en un determinado tipo de presa, con lo cual la otra se encuentra menos amenazada.

[[Archivo:GraphMEuler2.jpg|centro]]

==== Considerando la interacción entre especies ====
Para evaluar la interacción entre especies tendremos en cuenta las gráficas obtenidas para un tiempo '''tє[0,300]''', ya que se observa mejor la variación de las poblaciones de presas respecto de la de predadores.

[[Archivo:GraphMEuler5.jpg|centro]]

Estas gráficas nos permiten concluir que los puntos en los que el ecosistema sería estable son los orígenes (0,0), los cuales despreciamos (pues si lo interpretamos en términos de poblaciones, no hay especies interactuando), y los puntos donde se encuentran los centros de las espirales. Como podemos comprobar, ninguna de las dos curvas converge al mismo centro, así que no podemos asegurar que los ecosistemas sean estables.

== Resolución por el método trapezoidal ==
El método trapezoidal es un método implícito que, en nuestro caso, requiere la resolución de una ecuación no lineal. Al no ser posible su resolución, se descarta el método y no se tendrá en cuenta a la hora de comparar resultados.

== Resolución por el método de Runge-Kutta ==
=== Programación del método numérico ===
A continuación, se resuelve el problema de valor inicial mediante el método de Runge- Kutta de cuarto orden. Como datos tenemos una longitud de paso '''h=0.1''', constantes '''A1=0.4''', '''A2=0.7''', '''B1=0.2''', '''B2=0.4''', '''C1=0.37''', '''C2=0.04''' y '''C3=0.03''', y valores iniciales de '''p0=3.5''' millones de presas de un tipo, '''q0=0.2''' millones de presas de otro tipo y '''d0=0.4''' millones de predadores. En este caso tomaremos un intervalo de tiempo '''tє[0,500]'''. El programa de obtención del método numérico se especifica a continuación:
{{matlab|codigo=
t0=0; tN=500; %Extremos del intervalo
y0=[3.5;0.2;0.4]; %Valores de población iniciales
y=y0;
P1(1)=y(1); %Valor inicial de P1=presas tipo I
P2(1)=y(1); %Valor inicial de P2=presas tipo II
D(1)=y(1); %Valor inicial de D=depredadores
h=0.1; %Amplitud de los subintervalos
N=(tN-t0)/h; %Número de subintervalos
A=[ 0.4 0 0;0 0.2 0;0 0 -0.37]; %Parámetros
%Algoritmo
for k=1:N
k1=(A*y+y(1)*y(3)*[-0.7;0;0.04]+y(2)*y(3)*[0;-0.4;0.03]);
k2=(A*(y+(h/2)*k1)+(y(1)+(h/2)*k1(1))*(y(3)+(h/2)*k1(3))*[-0.7;0;0.04]+...
(y(2)+(h/2)*k1(2))*(y(3)+(h/2)*k1(3))*[0;-0.4;0.03]);
k3=(A*(y+(h/2)*k2)+(y(1)+(h/2)*k2(1))*(y(3)+(h/2)*k2(3))*[-0.7;0;0.04]+...
(y(2)+(h/2)*k2(2))*(y(3)+(h/2)*k2(3))*[0;-0.4;0.03]);
k4=(A*(y+h*k3)+(y(1)+h*k3(1))*(y(3)+h*k3(3))*[-0.7;0;0.04]+...
(y(2)+h*k3(2))*(y(3)+h*k3(3))*[0;-0.4;0.03]);
y=y+h/6*(k1+2*k2+2*k3+k4);
P1(k+1)=y(1);
P2(k+1)=y(2);
D(k+1)=y(3);
end

t=t0:h:tN;

%CÁLCULO DE INSTANTES MÁXIMOS Y MÍNIMOS

[maximopresaI,x1]= max(P1)
t1=t(x1)
[minimopresaI,x2]= min(P1)
t2=t(x2)
[maximopresaII,x3]= max(P2)
t3=t(x3)
[minimopresaII,x4]= min(P2)
t4=t(x4)
[maximopredadores,x5]= max(D)
t5=t(x5)
[minimopredadores,x6]= min(D)
t6=t(x6)

%GRÁFICA PRESA TIPO I/PRESA TIPO II/DEPREDADORES

hold on
figure(1)
plot(t,P1,'c')
plot(t,P2,'r')
plot(t,D,'g')
title(' Población de presas (tipos I y II) y depredadores')
xlabel('t(años)')
ylabel('población(millones de individuos)')
hold off
}}

=== Gráficas: Evolución de las especies en el tiempo ===
[[Archivo:GraficaRungeKutta.JPG|centro]]

=== Interpretación de resultados ===

En este gráfico se observan las variaciones poblacionales a lo largo del tiempo de las tres especies; la presa tipo I en color cian, la presa de tipo II en color rojo y los predadores se muestran en verde. Interpretamos que la población de predadores empieza a aumentar cuando la población de presas es máxima, es decir cuando tienen mayor cantidad de comida. Este crecimiento se mantiene hasta que el nivel de presas alcanza un punto crítico, donde deja de haber suficientes presas para sostener la población de predadores, lo que genera una disminución en el número de dichos depredadores, que a su vez da lugar a un aumento en el número de presas.

Podemos observar dos comportamientos distintos según el tipo de presa al que nos referimos:
:* En la presa de tipo II observamos a simple vista que la población inicial es muy reducida, incluso menor que la de depredadores. Esto provoca que a pesar de que este tipo de presa no sea el único del que se alimenta el depredador, esta especie termina por extinguirse, como bien se aprecia en la gráfica.
:*En la presa de tipo I, vemos que a parte de que su población inicial es muchísimo mayor que la de depredadores, su velocidad de crecimiento también es mayor, y esto asegura la supervivencia de ambas especies.
Como resultado conluimos que tanto la cantidad de depredadores como la de presas del tipo que no se han extinguido mantienen sus ciclos de crecemiento y decrecimiento por lo que ambas se mantienen constantes en lo que se refiere a máximos y mínimos.

Dichos máximos y mínimos se adjuntan en la siguiente tabla, con los correspondientes instantes en los que tienen lugar:

== Conclusiones ==
Tras analizar los resultados obtenidos aplicando los diferentes métodos numéricos, llegamos a las siguientes conclusiones:
* El método de Euler es un método de orden 1, por lo que nos ofrecerá peores resultados que los métodos de órdenes superiores, como son el método de Euler modificado (Runge-Kutta de orden 2) y el Runge-Kutta de orden 4. Esto se comprueba al comparar los distintos resultados obtenidos.
* En comparación con los otros métodos, es preferible el '''método de Runge-Kutta''', pues aporta soluciones más precisas. Esto se debe a que permite buenas aproximaciones aumentando el valor de la longitud de paso ''h'', ya que no se acumula el error del cálculo mediante tangentes que se desvían de la solución exacta.

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

Archivo:Maximosminimos.JPG

2014-03-04T10:20:23Z

Toni Villa: