Resolver ecuaciones presa-depredador

Abrir script en vivo

En este ejemplo se muestra cómo resolver una ecuación diferencial en la que se representa un modelo presa-depredador utilizando métodos de integración Runge-Kutta de tamaño de paso variable. El método ode23 utiliza un par de fórmulas de segundo y tercer orden para obtener una precisión intermedia y el método ode45 utiliza un par de cuarto y quinto orden para obtener una mayor precisión.

Veamos el par de ecuaciones diferenciales ordinarias de primer orden conocidas como ecuaciones Lotka-Volterra o modelo presa-depredador:

$\begin{array}{l} \frac{dx}{dt} = x - α xy \\ \frac{dy}{dt} = - y + β xy . \end{array}$

Las variables $x$ y $y$ miden los tamaños de las poblaciones de presas y depredadores, respectivamente. El término cruzado cuadrático representa las interacciones entre las especies. La población de presas aumenta cuando no hay depredadores y la población de depredadores disminuye cuando hay escasez de presas.

Ecuaciones de código

Para simular el sistema, cree una función que devuelva un vector columna de derivadas de estado, dados los valores de estado y tiempo. Las dos variables $x$ y $y$ se pueden representar en MATLAB® como los dos primeros valores de un vector y. Igualmente, las derivadas son los dos primeros valores de un vector yp. La función debe admitir valores de t e y y devolver los valores producidos por la ecuación en yp.

yp(1) = (1 - alpha*y(2))*y(1)

yp(2) = (-1 + beta*y(1))*y(2)

En este ejemplo, las ecuaciones se encuentran en un archivo llamado lotka.m. Este archivo utiliza valores de parámetro de $α = 0.01$ y $β = 0.02$ .

type lotka

function yp = lotka(t,y)
%LOTKA  Lotka-Volterra predator-prey model.

%   Copyright 1984-2014 The MathWorks, Inc.

yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

Simular el sistema

Cree un objeto ode para definir el problema y las condiciones iniciales. Utilice una condición inicial de $x (0) = y (0) = 20$ para que las poblaciones de depredadores y presas sean iguales inicialmente. Especifique que se debe usar el solver ode23. Después, utilice el método solve para simular el sistema en el intervalo de tiempo $0 < t < 15$ .

y0 = [20; 20];   
F = ode(ODEFcn=@lotka,InitialValue=y0,Solver="ode23")

F = 
  ode with properties:

   Problem definition
               ODEFcn: @lotka
          InitialTime: 0
         InitialValue: [2×1 double]
         EquationType: standard

   Solver properties
    AbsoluteTolerance: 1.0000e-06
    RelativeTolerance: 1.0000e-03
               Solver: ode23

  Show all properties

t0 = 0;
tfinal = 15;
S = solve(F,t0,tfinal);
t = S.Time;
y = S.Solution;

Representar los resultados

Represente las poblaciones resultantes con respecto al tiempo.

plot(t,y,"-o")
title("Predator/Prey Populations Over Time")
xlabel("t")
ylabel("Population")
legend("Prey","Predators",Location="north")

Figure contains an axes object. The axes object with title Predator/Prey Populations Over Time, xlabel t, ylabel Population contains 2 objects of type line. These objects represent Prey, Predators.

Ahora, represente las poblaciones una en comparación con la otra. La gráfica de plano de fases resultante muestra la relación cíclica entre las poblaciones de forma muy clara.

plot(y(1,:),y(2,:))
title("Phase Plane Plot")
xlabel("Prey Population")
ylabel("Predator Population")

Figure contains an axes object. The axes object with title Phase Plane Plot, xlabel Prey Population, ylabel Predator Population contains an object of type line.

Comparar los resultados de diferentes solvers

Resuelva el sistema una segunda vez utilizando el solver ode45 en lugar de ode23. El solver ode45 tarda más en cada paso, pero también realiza pasos más largos. Sin embargo, el valor de salida de ode45 es fluido, porque de forma predeterminada el solver utiliza una fórmula de extensión continua para producir el valor de salida en cuatro puntos temporales igualmente espaciados en el intervalo de cada paso realizado (puede ajustar el número de puntos con el argumento nombre-valor Refine del método solve). Represente las dos soluciones para hacer una comparación.

F.Solver = "ode45";
S2 = solve(F,t0,tfinal);
t2 = S2.Time;
y2 = S2.Solution;

plot(y(1,:),y(2,:),"o",y2(1,:),y2(2,:),".");
title("Phase Plane Plot")
legend("ode23","ode45")

Figure contains an axes object. The axes object with title Phase Plane Plot contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent ode23, ode45.

Los resultados muestran que resolver ecuaciones diferenciales con diferentes métodos numéricos puede producir respuestas ligeramente diferentes.

Consulte también

ode45 | ode23