Predator-Prey System

4230 days ago by medlock

Linear Functional Response
var('N P alpha beta') dNdt = N * (1 - N - P) dPdt = beta * P * (N - alpha) show(dNdt) show(dPdt) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(N + P - 1\right)} N


\newcommand{\Bold}[1]{\mathbf{#1}}{\left(N - \alpha\right)} P \beta
\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(N + P - 1\right)} N


\newcommand{\Bold}[1]{\mathbf{#1}}{\left(N - \alpha\right)} P \beta
equilibria = solve([dNdt == 0, dPdt == 0], [N, P], solution_dict = True) show(equilibria) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{P:\: 0, N:\: 0\right\}, \left\{P:\: 0, N:\: 1\right\}, \left\{P:\: -\alpha + 1, N:\: \alpha\right\}\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{P:\: 0, N:\: 0\right\}, \left\{P:\: 0, N:\: 1\right\}, \left\{P:\: -\alpha + 1, N:\: \alpha\right\}\right]
J = jacobian([dNdt, dPdt], [N, P]) show(J) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr}
-2 \, N - P + 1 & -N \\
P \beta & {\left(N - \alpha\right)} \beta
\end{array}\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr}
-2 \, N - P + 1 & -N \\
P \beta & {\left(N - \alpha\right)} \beta
\end{array}\right)
for e in equilibria: show(J(e).eigenvalues()) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\alpha \beta, 1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[-{\left(\alpha - 1\right)} \beta, -1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, \alpha - \frac{1}{2} \, \sqrt{4 \, {\left(\alpha^{2} - \alpha\right)} \beta + \alpha^{2}}, -\frac{1}{2} \, \alpha + \frac{1}{2} \, \sqrt{4 \, {\left(\alpha^{2} - \alpha\right)} \beta + \alpha^{2}}\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\alpha \beta, 1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[-{\left(\alpha - 1\right)} \beta, -1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{2} \, \alpha - \frac{1}{2} \, \sqrt{4 \, {\left(\alpha^{2} - \alpha\right)} \beta + \alpha^{2}}, -\frac{1}{2} \, \alpha + \frac{1}{2} \, \sqrt{4 \, {\left(\alpha^{2} - \alpha\right)} \beta + \alpha^{2}}\right]
NMin = 0 NMax = 1.1 PMin = 0 PMax = 1.1 def phasePortrait(alpha, beta): VF = plot_vector_field((dNdt(alpha = alpha, beta = beta), dPdt(alpha = alpha, beta = beta)), (N, NMin, NMax), (P, PMin, PMax)) EQ = point2d([(e[N](alpha = alpha, beta = beta), e[P](alpha = alpha, beta = beta)) for e in equilibria], size = 25) NNull = line([(0, PMin), (0, PMax)], color = 'red') + line([(NMin, 1 - NMin), (1 - PMin, PMin)], color = 'red') PNull = line([(alpha, PMin), (alpha, PMax)], color = 'green') + line([(NMin, 0), (NMax, 0)], color = 'green') return EQ + NNull + PNull + VF 
       
def ODERHS(t, Y, params): (N, P) = Y (alpha, beta) = params return (dNdt(N = N, P = P, alpha = alpha, beta = beta), dPdt(N = N, P = P, alpha = alpha, beta = beta)) def computeSolution(alpha, beta, IC): T = ode_solver(function = ODERHS, algorithm = 'rk8pd') T.ode_solve(y_0 = IC, t_span=[0, 50], params = (alpha, beta), num_points = 200) (t, Y) = zip(*T.solution) (N, P) = zip(*Y) return (t, N, P) def numericalSolution(alpha, beta, IC): (t, Nn, Pn) = computeSolution(alpha, beta, IC) return line(zip(t, Nn), color = 'red') + line(zip(t, Pn), color = 'green') def numericalPhaseTrajectory(alpha, beta, IC): (t, Nn, Pn) = computeSolution(alpha, beta, IC) return line(zip(Nn, Pn)) ICs = [(n, 0.99) for n in (0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)] + [(0.99, p) for p in (0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)] def numericalPhaseTrajectories(alpha, beta): return phasePortrait(alpha, beta) + sum(numericalPhaseTrajectory(alpha, beta, IC) for IC in ICs) 
       
alpha = 0.8 beta = 0.5 show(phasePortrait(alpha, beta)) show(numericalSolution(alpha, beta, [random(), random()])) show(numericalPhaseTrajectories(alpha, beta)) 
       

Type II Functional Response
var('N P alpha beta gamma') dNdt = N * (1 - N / gamma - P / (1 + N)) dPdt = beta * P * (N / (1 + N) - alpha) show(dNdt) show(dPdt) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(\frac{P}{N + 1} + \frac{N}{\gamma} - 1\right)} N


\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(\alpha + \frac{-N}{N + 1}\right)} P \beta
\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(\frac{P}{N + 1} + \frac{N}{\gamma} - 1\right)} N


\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(\alpha + \frac{-N}{N + 1}\right)} P \beta
equilibria = solve([dNdt == 0, dPdt == 0], [N, P], solution_dict = True) show(equilibria) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{P:\: 0, N:\: 0\right\}, \left\{P:\: 0, N:\: -1\right\}, \left\{P:\: 0, N:\: \gamma\right\}, \left\{P:\: \frac{-{\left(\alpha - 1\right)} \gamma + \alpha}{{\left(\alpha^{2} - 2 \, \alpha + 1\right)} \gamma}, N:\: \frac{-\alpha}{\alpha - 1}\right\}\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{P:\: 0, N:\: 0\right\}, \left\{P:\: 0, N:\: -1\right\}, \left\{P:\: 0, N:\: \gamma\right\}, \left\{P:\: \frac{-{\left(\alpha - 1\right)} \gamma + \alpha}{{\left(\alpha^{2} - 2 \, \alpha + 1\right)} \gamma}, N:\: \frac{-\alpha}{\alpha - 1}\right\}\right]
equilibria.pop(1) show(equilibria) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{P:\: 0, N:\: 0\right\}, \left\{P:\: 0, N:\: \gamma\right\}, \left\{P:\: \frac{-{\left(\alpha - 1\right)} \gamma + \alpha}{{\left(\alpha^{2} - 2 \, \alpha + 1\right)} \gamma}, N:\: \frac{-\alpha}{\alpha - 1}\right\}\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{P:\: 0, N:\: 0\right\}, \left\{P:\: 0, N:\: \gamma\right\}, \left\{P:\: \frac{-{\left(\alpha - 1\right)} \gamma + \alpha}{{\left(\alpha^{2} - 2 \, \alpha + 1\right)} \gamma}, N:\: \frac{-\alpha}{\alpha - 1}\right\}\right]
J = jacobian([dNdt, dPdt], [N, P]) show(J) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr}
{\left(\frac{P}{{\left(N + 1\right)}^{2}} + \frac{-1}{\gamma}\right)} N + \frac{-P}{N + 1} + \frac{-N}{\gamma} + 1 & \frac{-N}{N + 1} \\
{\left(\frac{1}{N + 1} + \frac{-N}{{\left(N + 1\right)}^{2}}\right)} P \beta & -{\left(\alpha + \frac{-N}{N + 1}\right)} \beta
\end{array}\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr}
{\left(\frac{P}{{\left(N + 1\right)}^{2}} + \frac{-1}{\gamma}\right)} N + \frac{-P}{N + 1} + \frac{-N}{\gamma} + 1 & \frac{-N}{N + 1} \\
{\left(\frac{1}{N + 1} + \frac{-N}{{\left(N + 1\right)}^{2}}\right)} P \beta & -{\left(\alpha + \frac{-N}{N + 1}\right)} \beta
\end{array}\right)
for e in equilibria: show(J(e).eigenvalues()) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\alpha \beta, 1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{-{\left(\alpha - 1\right)} \beta \gamma + \alpha \beta}{\gamma + 1}, -1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{{\left(\alpha^{2} - \alpha\right)} \gamma + \alpha - \sqrt{\alpha^{4} + {\left(\alpha^{4} - 2 \, \alpha^{3} + 4 \, {\left(\alpha^{4} - 3 \, \alpha^{3} + 3 \, \alpha^{2} - \alpha\right)} \beta + \alpha^{2}\right)} \gamma^{2} + 2 \, \alpha^{3} + 2 \, {\left(\alpha^{4} + 2 \, {\left(\alpha^{4} - 2 \, \alpha^{3} + \alpha^{2}\right)} \beta - \alpha^{2}\right)} \gamma + \alpha^{2}} + \alpha^{2}}{2 \, {\left(\alpha - 1\right)} \gamma}, \frac{{\left(\alpha^{2} - \alpha\right)} \gamma + \alpha + \sqrt{\alpha^{4} + {\left(\alpha^{4} - 2 \, \alpha^{3} + 4 \, {\left(\alpha^{4} - 3 \, \alpha^{3} + 3 \, \alpha^{2} - \alpha\right)} \beta + \alpha^{2}\right)} \gamma^{2} + 2 \, \alpha^{3} + 2 \, {\left(\alpha^{4} + 2 \, {\left(\alpha^{4} - 2 \, \alpha^{3} + \alpha^{2}\right)} \beta - \alpha^{2}\right)} \gamma + \alpha^{2}} + \alpha^{2}}{2 \, {\left(\alpha - 1\right)} \gamma}\right]
\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\alpha \beta, 1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{-{\left(\alpha - 1\right)} \beta \gamma + \alpha \beta}{\gamma + 1}, -1\right]


\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{{\left(\alpha^{2} - \alpha\right)} \gamma + \alpha - \sqrt{\alpha^{4} + {\left(\alpha^{4} - 2 \, \alpha^{3} + 4 \, {\left(\alpha^{4} - 3 \, \alpha^{3} + 3 \, \alpha^{2} - \alpha\right)} \beta + \alpha^{2}\right)} \gamma^{2} + 2 \, \alpha^{3} + 2 \, {\left(\alpha^{4} + 2 \, {\left(\alpha^{4} - 2 \, \alpha^{3} + \alpha^{2}\right)} \beta - \alpha^{2}\right)} \gamma + \alpha^{2}} + \alpha^{2}}{2 \, {\left(\alpha - 1\right)} \gamma}, \frac{{\left(\alpha^{2} - \alpha\right)} \gamma + \alpha + \sqrt{\alpha^{4} + {\left(\alpha^{4} - 2 \, \alpha^{3} + 4 \, {\left(\alpha^{4} - 3 \, \alpha^{3} + 3 \, \alpha^{2} - \alpha\right)} \beta + \alpha^{2}\right)} \gamma^{2} + 2 \, \alpha^{3} + 2 \, {\left(\alpha^{4} + 2 \, {\left(\alpha^{4} - 2 \, \alpha^{3} + \alpha^{2}\right)} \beta - \alpha^{2}\right)} \gamma + \alpha^{2}} + \alpha^{2}}{2 \, {\left(\alpha - 1\right)} \gamma}\right]
def phasePortrait(alpha, beta, gamma): NMin = 0 NMax = 1.1 * gamma PMin = 0 PMax = 1.2 * (gamma + 1)^2 / 4 / gamma VF = plot_vector_field((dNdt(alpha = alpha, beta = beta, gamma = gamma), dPdt(alpha = alpha, beta = beta, gamma = gamma)), (N, NMin, NMax), (P, PMin, PMax)) EQ = point2d([(e[N](alpha = alpha, beta = beta, gamma = gamma), e[P](alpha = alpha, beta = beta, gamma = gamma)) for e in equilibria], size = 25) NNull = line([(0, PMin), (0, PMax)], color = 'red') \ + plot((1 + N) * (1 - N / gamma), (N, NMin, NMax), color = 'red') PNull = line([(alpha / (1 - alpha), PMin), (alpha / (1 - alpha), PMax)], color = 'green') \ + line([(NMin, 0), (NMax, 0)], color = 'green') return EQ + NNull + PNull + VF 
       
def ODERHS(t, Y, params): (N, P) = Y (alpha, beta, gamma) = params return (dNdt(N = N, P = P, alpha = alpha, beta = beta, gamma = gamma), dPdt(N = N, P = P, alpha = alpha, beta = beta, gamma = gamma)) def computeSolution(alpha, beta, gamma, IC): T = ode_solver(function = ODERHS, algorithm = 'rk8pd') T.ode_solve(y_0 = IC, t_span=[0, 100], params = (alpha, beta, gamma), num_points = 200) (t, Y) = zip(*T.solution) (N, P) = zip(*Y) return (t, N, P) def numericalSolution(alpha, beta, gamma, IC): (t, Nn, Pn) = computeSolution(alpha, beta, gamma, IC) return line(zip(t, Nn), color = 'red') + line(zip(t, Pn), color = 'green') def numericalPhaseTrajectory(alpha, beta, gamma, IC): (t, Nn, Pn) = computeSolution(alpha, beta, gamma, IC) return line(zip(Nn, Pn)) ICs = [(n, 0.99) for n in (0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)] + [(0.99, p) for p in (0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99)] def numericalPhaseTrajectories(alpha, beta, gamma): return phasePortrait(alpha, beta, gamma) \ + sum(numericalPhaseTrajectory(alpha, beta, gamma, (1.1 * gamma * N0, 1.2 * (gamma + 1)^2 / 4 / gamma * P0)) for (N0, P0) in ICs) 
       
alpha = 0.65 beta = 0.5 gamma = 5 show(phasePortrait(alpha, beta, gamma)) show(numericalSolution(alpha, beta, gamma, [gamma * random(), (gamma + 1)^2 / 4 / gamma * random()])) show(numericalPhaseTrajectories(alpha, beta, gamma))