Modelo Depredador-Presa de Lotka-Volterra (grupo 16)

1 Resolución mediante el método Runge-Kutta

A continuación se muestra el programa usado para la resolución numérica de nuestro problema mediante el método de Runge-Kutta de 4º orden para 1000 intervalos.

%Valores iniciales
t0=0; tN=100; %Extremos del intervalo a estudiar
y0=[3,1]'; %Valores de población iniciales
N=1000; h=(tN-t0)/N;%Nº y  amplitud de subintervalos
y=y0;
R(1)=y(1);%Valor inicial de R (población de conejos)
F(1)=y(2);% Valor inicial de F (población de zorros)
a=0.4; b=0.37; c=0.3; d=0.05; %Parámetros
%Algoritmo
for n=1:N
    k1=[a*y(1)-c*y(1)*y(2);-b*y(2)+d*y(2)*y(1)];
    k2=[a*(y(1)+1/2*h*k1(1))-c*(y(1)+1/2*h*k1(1))*(y(2)+1/2*h*k1(2));-b*(y(2)+1/2*h*k1(2))+d*(y(2)+1/2*h*k1(2))*(y(1)+1/2*h*k1(1))];
    k3=[a*(y(1)+1/2*h*k2(1))-c*(y(1)+1/2*h*k2(1))*(y(2)+1/2*h*k2(2));-b*(y(2)+1/2*h*k2(2))+d*(y(2)+1/2*h*k2(2))*(y(1)+1/2*h*k2(1))];
    k4=[a*(y(1)+h*k3(1))-c*(y(1)+h*k3(1))*(y(2)+h*k3(2));-b*(y(2)+h*k3(2))+d*(y(2)+h*k3(2))*(y(1)+h*k3(1))];
    y=y+h/6*(k1+2*k2+2*k3+k4);
    R(n+1)=y(1);
    F(n+1)=y(2);
end
%Dibujo
t=t0:h:tN;
figure(1)
hold on
plot(t,R,'g')%Población de conejos en función del tiempo
plot(t,F,'r')%Población de zorros en función del tiempo
title('Población de conejos(verde) y zorros(rojo) en función del tiempo');
ylabel('t(años)');
hold off
figure(2)
plot(R,F,'b')%Población de zorros en función de la de conejos
title('Población de zorros en función de la de conejos');
xlabel('R(miles de conejos)');
ylabel('F(miles de zorros)');

Obtenemos como resultado las siguientes gráficas:

• 
• 

1.1 Ventajas

En comparación con otros métodos de resolución numérica, el de Runge-Kutta tiene ciertas propiedades que pueden hacerlo preferible frente a los otros.

• Si comparamos con el método de Euler tomando diferentes amplitudes (h) de subintervalo en ambos métodosobservamos que el método de Runge-Kutta aporta soluciones más precisas.
  

Esto puede deberse a que con el método de Euler la aproximación se realiza por la recta tangente y por ello, en soluciones cíclicas u oscilatorias, como es nuestro caso, no dan muy buenos resultados a no ser que h sea muy pequeña con respecto al intervalo de estudio. Esto se explica porque el método de Euler alcanza un mayor nivel de convergencia cuanto menor es el valor de h.

Destacar al respecto que podría surgir el problema de que al escoger un valor de h demasiado pequeño el programa no responda, por lo que es necesario encontrar un valor adecuado de h que aporte una buena precisión y que pueda ser soportado.

Sin embargo, esto no ocurre con el método de Runge-Kutta, que permite buenas aproximaciones con mayores valores de h, ya que no depende de tangentes que se desvían de la solución acumulando error que puede llegar a hacer inservible la solución final. Así lo vemos en las siguientes gráficas para las que se ha tomado un valor de h=1:

Con el método de Euler:

• 
• 

Con Runge-Kutta:

• 
• 
• *El empleo del método de los trapecios en la resolución de nuestro sistema resulta inviable, ya que se trata de un sistema no lineal.