Predator-prey system¶
We encode the predator–prey equations via a function.
using FundamentalsNumericalComputation
function predprey(u,p,t)
α,β = p; y,z = u; # rename for convenience
s = (y*z) / (1+β*y) # appears in both equations
return [ y*(1-α*y) - s, -z + s ]
end
predprey (generic function with 1 method)
Note that the function accepts three inputs, u, p, and t, even though there is no explicit dependence on t. The second input is used to pass parameters that don’t change throughout a single instance of the problem.
To solve the IVP we must also provide the initial condition, which is a 2-vector here, and the interval for the independent variable.
u0 = [1,0.01]
tspan = (0.,80.)
α,β = 0.1,0.25
ivp = ODEProblem(predprey,u0,tspan,[α,β])
sol = solve(ivp,Tsit5());
plot(sol,label=["prey" "predator"],title="Predator-prey solution")
We can find the discrete values used to compute the interpolated solution. The sol.u value is a vector of vectors, which can make manipulating the values a bit tricky. Here we convert the solution values to a matrix with two columns (one for each component).
@show size(sol.u);
@show (sol.t[20],sol.u[20]);
size(sol.u) = (133,)
(sol.t[20], sol.u[20]) =
(8.38669870359271, [0.027739196929331474, 0.6994179006495149])
u = [ sol.u[i][j] for i=1:length(sol.t), j=1:2 ]
scatter!(sol.t,u,m=3,label="")
When there are just two components, it’s common to plot the solution in the phase plane, i.e., with \(u_1\) and \(u_2\) along the axes and time as a parameterization of the curve.
plot(sol,vars=(1,2),label="",
xlabel="y",ylabel="z",title="Predator-prey phase plane")
From this plot we can deduce that the solution approaches a periodic one, which in the phase plane is reprepresented by a closed loop.