Resolver ecuaciones presa-depredador
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:
Las variables 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 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 y .
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 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 .
y0 = [20; 20];
F = ode(ODEFcn=@lotka,InitialValue=y0,Solver="ode23")
F = ode with properties: Problem definition ODEFcn: @lotka InitialTime: 0 InitialValue: [2x1 double] 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")
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")
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")
Los resultados muestran que resolver ecuaciones diferenciales con diferentes métodos numéricos puede producir respuestas ligeramente diferentes.