2

I'm working on my thesis on prey-predator system. I have 3 species here, 1 prey and 1 predator with 1 species having a mutualism with the prey and protecting the prey from predators. Here is the system :
$$\frac{dx}{dt}=r_1x(1-\frac{x}{\gamma})-\frac{xz}{1+\alpha N+x}+\theta_1xy$$ $$\frac{dy}{dt}=r_2y(1-\frac{y}{\gamma})+\theta_2xy-\Delta_1yz$$ $$\frac{dz}{dt}=\frac{Bz(x+N)}{1+\alpha N+x}-\Delta_2yz-Mz$$ Finding equilibria without the presence of predator is pretty easy. So far 2 equilibrium $E_0=(0,0,0)$ and $E_1=(x_1^*,y_1^*,0)$ has been found. The problem is for the third equilibria where all species live together. It is really hard to find the closed form solutions of the third equilibria. So, I stop here at these equations by $f_i(x,y,z)=0$ we get : $$r_1x(1-\frac{x}{\gamma})-\frac{xz}{1+\alpha N+x}+\theta_1xy=0$$ $$\frac{B(x+N)}{1+\alpha N+x}-\Delta_2y-M=0$$ $$\frac{r_2}{\Delta_1}(1-\frac{y}{\gamma})+\frac{\theta_2}{\Delta_1}x=z$$

From here I don't know how to proceed to analyze the stability of the third equilibrium since how can you input the equation to the jacobian matrix of the system. How shall I approach on this problem? Can I do Newton-Raphson or any numerical approximation in my problem here? I wish I could somehow make it shorter.
From above, i try to make $x$ as a parameter and solving for $y$ and $z$ that gives me polynomial of $x$ in degree of 3. I tried using software Mathematica to work on the solution of this system. Mathematica give me 3 non-trivial equilibrium (this is because the polynomial in 3rd degree) which is a very large output. I couldn't find online how to analyze the stability of a non linear system like this. Any suggestions are very much appreciated. Thank you.

$r_i$ is the rate of growth species-i. $\gamma$ is the system capacity that can hold how much of the species. $alpha$ is the additional food quality. $N$ is the quantity of the additional food. $B$ is maximum growth of predator. $\theta_i$ is the amount of profit from mutualism to the species. $\Delta_i$ is the rate of death because interaction of predator and species 2. And $M$ is the rate of death of predator.

Edit : all symbols are all constant parameters, the only variables are $x$, $y$, and $z$

Edit 2 : Adding additional info on the parameters.

Aji Wibowo
  • 45
  • 1
    Are $B, \alpha, \Delta_2, M$ constants? Maybe it would help to define all constants in the problem statement. – Moo Dec 10 '24 at 21:18
  • @Moo yes they are all constants, the only variable is $x,y$, and $z$ – Aji Wibowo Dec 11 '24 at 06:44
  • 1
    You already found $z$ in the 3rd equation. You can single out $y$ from the 2nd equation. Putting everything together in the 1st equation you can obtain $x$ in terms of constants. – obareey Dec 11 '24 at 09:11
  • @obareey yes, but it gives more question to me. I already did that as i said, if i try to take $y$ and $z$ depends on x, there will be a polynomial function $f(x)$ with degree of 3. This is not an easy task as there are so many parameters. Furthermore, if we leave $f(x)$ as it is, we arrive at my question above, can we input this polynomial $f(x)$, $y(x,t)$, and $z(x,t)$ into the jacobian matrix? I think this is not possible and there need to be another approach on how to tackle this problem – Aji Wibowo Dec 11 '24 at 15:19
  • You don't need to input the polynomial into the Jacobian. Solving the polynomial will give you the fixed points, which are constant values. When calculating the Jacobian around a fixed point, you enter the value of the selected fixed point into the Jacobian for linearization. Jacobian matrix is calculated by the partial derivatives of $x,y,z$, which should be straightforward. – obareey Dec 11 '24 at 15:38
  • @obareey the problem is the root of the polynominal are humongous. I have tried to find it with software like Mathematica and Maplesoft. Couldn't even write all that down on paper – Aji Wibowo Dec 11 '24 at 16:56

1 Answers1

1

Hint.

To use criteria as Routh-Hurwitz the system should be linearized. Let the nonlinear system be represented as

$$ \dot x = f(x,\theta) $$

at $x_0$ we have the linear approximation

$$ \dot x = f(x_0,\theta)+\nabla_x f(x_0,\theta)(x-x_0) = A(x_0,\theta)x+b(x_0,\theta) $$

applying the Laplace transform with null initial conditions we get

$$ X(s) = (I s - A(x_0,\theta))^{-1}\frac{b(x_0,\theta)}{s} $$

and the stability is determined by the roots of

$$ \det(I s - A(x_0,\theta))=0 $$

which is a third order polynomial. Applying now the Routh-Hurwitz criterion for instance, for stable roots, we will obtain a set of conditions $g_i(x_0,\theta) \ge 0$ so the problem now is to determine sets of $\{x_0,\theta\}$ which obey those conditions. This last step can be done with the help of an optimization solver.

Cesareo
  • 36,341