Stability analysis examines how the solutions of a system behave as time progresses, particularly in response to changes in initial conditions or parameters. Following are the key steps involved in stability analysis of system Eq. (5).
Local stability
The possible equilibria of model Eq. (5) are given below as
-
1.
(S_{1} left( {frac{Q}{sigma },0,0,frac{Qbeta + Asigma }{{dsigma }}} right)),
-
2.
(S_{2} = left( {frac{{rleft( {Q + Leta } right)}}{Leta theta + rsigma },,frac{{Lleft( {rsigma – Qtheta } right)}}{Leta theta + rsigma },,0,,frac{Qrbeta + Lrbeta eta + ALeta theta – LQtheta lambda + Arsigma + Lrlambda sigma }{{dleft( {Leta theta + rsigma } right)}}} right)),
-
3.
(S_{3} = left( {C_{3} ,,0,,V_{3} ,,G_{3} } right)),
-
4.
(S_{4} = left( {C_{4} ,,N_{4} ,,V_{4} ,,G_{4} } right)),
where
(C_{3} = frac{{Qalpha + L_{1} delta_{1} left( {alpha – omega } right)}}{alpha sigma }), (V_{3} = L_{1} – frac{{L_{1} omega }}{alpha }), (G_{3} = frac{{Qalpha beta + Aalpha sigma + L_{1} beta delta_{1} left( {alpha – omega } right)}}{dalpha sigma }), (C_{4} = frac{{rleft( {left. {Qalpha + Lalpha eta + L_{1} delta_{1} left( alpha right. – omega } right)} right)}}{{alpha left( {Leta theta + rsigma } right)}}), (C_{4} = frac{{Lleft( {left. { – Qalpha theta + ralpha sigma + L_{1} delta_{1} theta left( alpha right. – omega } right)} right)}}{{alpha left( {Leta theta + rsigma } right)}}), and (C_{4} = frac{{Qalpha left( {left. {rbeta – Ltheta lambda } right) + sigma left( {Lrbeta eta + ALeta theta + Arsigma + Lrlambda sigma } right) + L_{1} delta_{1} left( {rbeta – Ltheta lambda } right)left( {alpha – omega } right.} right)}}{{alpha dleft( {Leta theta + rsigma } right)}}).
The local stability at equilibrium points (S_{1}), (S_{2}), (S_{3}) and (S_{4}) is determined through the sign of eigenvalues of the following Jacobian matrix (J) that is defined as
$$J = left[ {begin{array}{*{20}l} { – sigma } hfill & eta hfill & {delta_{1} } hfill & 0 hfill \ { – Ntheta } hfill & {frac{ – Nr}{L} + left( {1 – frac{N}{L}} right) – Ctheta } hfill & 0 hfill & 0 hfill \ 0 hfill & 0 hfill & {frac{ – Valpha }{{L_{1} }} + left( {1 – frac{V}{{L_{1} }}} right)alpha – omega } hfill & 0 hfill \ beta hfill & lambda hfill & 0 hfill & { – d} hfill \ end{array} } right].$$
(6)
Take (J_{i}) (left( {i = 1,2,3,4} right)) Jacobian matrix that is determined at equilibrium point (S_{i} left( {i = 1,2,3,4} right)).
For (S_{1}), the eigenvalues of (J_{1}) are (left( { – d,,, – sigma ,,,frac{rsigma – Qtheta }{sigma },,,alpha – omega } right)), and the system is unstable at this equilibrium point if (alpha > omega). Moreover, the eigenvalues of the (J_{2}) are
$$left( { – d,,frac{{ – b + sqrt {b^{2} – 4ac} }}{2a},frac{{ – b – sqrt {b^{2} – 4ac} }}{2a},alpha – omega } right),$$
(7)
where (a = Leta theta + rsigma), (b = Qrtheta – r^{2} sigma – Leta theta sigma – rsigma^{2}), and (c = – LQeta theta^{2} – Qrtheta sigma + Lreta theta sigma + r^{2} sigma^{2}). Thus, the equilibrium point (S_{2}) is unstable if (alpha > omega).
The characteristic equation of (J_{i} (i = 3,4)) for (S_{i} (i = 3,4)) is given by
$$x^{4} + A_{1} x^{3} + A_{2} x^{2} + A_{3} x + A_{4} = 0,$$
(8)
where
$$A_{1} = frac{1}{{LL_{1} }}left( {dLL_{1} – LL_{1} r + 2L_{1} N_{i} r + LL_{1} sigma + C_{i} LL_{1} sigma + LVleft( {V_{i} – L_{1} } right)alpha omega } right),$$
$$A_{2} = frac{1}{{LL_{1} }}left( begin{gathered} – dLL_{1} r + 2dL_{1} N_{i} r + dLL_{1} sigma + C_{i} dLL_{1} sigma – LLrsigma + 2L_{1} N_{i} rsigma + LL_{1} N_{i} eta sigma hfill \ + C_{i} LL_{1} sigma^{2} + dLV_{i} left( {V_{i} – L_{1} } right)alpha omega – LrV_{i} left( {V_{i} – L_{1} } right)alpha omega + 2N_{i} rV_{i} left( {V_{i} – L_{1} } right)alpha omega hfill \ + LV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega + C_{i} LVleft( {V_{i} – L_{1} } right)alpha sigma omega hfill \ end{gathered} right),$$
$$A_{3} = frac{1}{{LL_{1} }}left( begin{gathered} – dLL_{1} rsigma + 2dL_{1} N_{i} rsigma + dLL_{1} N_{i} eta sigma + C_{i} dLL_{1} sigma^{2} – dLrV_{i} left( {V_{i} – L_{1} } right)alpha omega hfill \ + 2dN_{i} rV_{i} left( {V_{i} – L_{1} } right)alpha omega + dLV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega + C_{i} dLV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega hfill \ – LrV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega + 2N_{i} rV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega + LN_{i} V_{i} left( {V_{i} – L_{1} } right)alpha eta sigma omega hfill \ + C_{i} LV_{i} left( {V_{i} – L_{1} } right)alpha sigma^{2} omega hfill \ end{gathered} right),$$
$$A_{4} = frac{1}{{LL_{1} }}left( begin{gathered} – dLrV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega + 2dN_{i} rV_{i} left( {V_{i} – L_{1} } right)alpha sigma omega hfill \ + dLN_{i} V_{i} left( {V_{i} – L_{1} } right)alpha eta sigma omega + C_{i} dLV_{i} left( {V_{i} – L_{1} } right)alpha sigma^{2} omega hfill \ end{gathered} right),$$
where the coefficients (A_{i} left( {i = 1,2,3,4} right)) are positive. Equation (8) has either negative or positive eigenvalues iff the following Routh–Hurwitz condition is satisfied.
$$A_{3} left( {A_{1} A_{2} – A_{3} } right) – A_{1}^{2} A_{4} > 0$$
(9)
Therefore, the system is stable at (S_{i} (i = 3,4)) if Eq. (9) holds.
Theorem 1
The system in Eq. (5) at equilibriums (S_{1}) and (S_{2}) is always unstable under the condition (alpha > omega). The system at equilibrium points (S_{i}) is locally asymptotically stable iff Eq. (9) is hold.
Global stability
Now, stability of the Liapunov function is applied to determine the global stability. The following theorem illustrates the conditions of global stability.
Theorem 2
The Eq. (5) is globally stable in region (Omega) if following conditions are hold
$$m_{2} < frac{{2delta_{1}^{2} L_{1} }}{alpha sigma },$$
(10)
$$m_{3} < min left{ {frac{eta rd}{{2theta Llambda^{2} }},frac{dsigma }{{4beta^{2} }}} right}.$$
(11)
Proof
Consider the following positive function:
$$V = frac{1}{2}left( {C – C^{ * } } right)^{2} + m_{1} left( {N – N^{ * } – N^{ * } ln frac{N}{{N^{ * } }}} right) + m_{2} left( {V – V^{ * } – V^{ * } ln frac{V}{{V^{ * } }}} right) + frac{{m_{3} }}{2}left( {G – G^{ * } } right)^{2} ,$$
(12)
where (m_{1}), (m_{2}) and (m_{3}) are positive constants. Equation (5) is globally stable if (frac{dV}{{dt}} < 0) at all equilibrium points. Therefore, the derivative of Eq. (11) is calculated as
$$frac{dV}{{dt}} = left( {C – C^{ * } } right)frac{dC}{{dt}} + m_{1} frac{{left( {N – N^{ * } } right)}}{N}frac{dN}{{dt}} + m_{2} frac{{left( {V – V^{ * } } right)}}{V}frac{dV}{{dt}} + m_{3} left( {G – G^{ * } } right)frac{dG}{{dt}}$$
(13)
Using Eqs. (1–4), we obtain
$$begin{aligned} & frac{dV}{{dt}} = left( {C – C^{ * } } right)left[ {Q + delta_{1} V + eta N – sigma C} right] + m_{1} left( {N – N^{ * } } right)left[ {rleft( {1 – frac{N}{L}} right) – theta C} right] \ & quad + ,m_{2} left( {V – V^{ * } } right)left[ {alpha left( {1 – frac{V}{{L_{1} }}} right) – omega } right] + m_{3} left( {G – G^{ * } } right)left[ {A + beta C + lambda N – dG} right]. \ end{aligned}$$
(14)
Equation (14) is rewritten when the condition for finding the equilibrium point (S^{*} left( {C^{*} ,N^{*} ,V^{*} ,G^{*} } right)) for Eq. (5) is
$$begin{aligned} frac{dV}{{dt}} & = left( {C – C^{ * } } right)left[ {Q + delta_{1} V + eta N – sigma C – left{ {Q + delta_{1} V^{ * } + eta N^{ * } – sigma C^{ * } } right}} right] \ & quad + ,m_{1} left( {N – N^{ * } } right)left[ {rleft( {1 – frac{N}{L}} right) – theta C – left{ {rleft( {1 – frac{{N^{ * } }}{L}} right) – theta C^{ * } } right}} right]. \ & quad + ,m_{2} left( {V – V^{ * } } right)left[ {alpha left( {1 – frac{V}{{L_{1} }}} right) – omega – left{ {alpha left( {1 – frac{{V^{ * } }}{{L_{1} }}} right) – omega } right}} right] \ & quad + ,m_{3} left( {G – G^{ * } } right)left[ {A + beta C – dG + lambda N – left{ {A + beta C^{ * } – dG^{ * } + lambda N^{ * } } right}} right]. \ end{aligned}$$
(15)
Rewriting Eq. (15) as
$$begin{aligned} frac{dV}{{dt}} & = left( {C – C^{ * } } right)left[ {delta_{1} left( {V – V^{ * } } right) + eta left( {N – N^{ * } } right) – sigma left( {C – C^{ * } } right)} right] \ & quad + ,m_{1} left( {N – N^{ * } } right)left[ {frac{ – r}{L}left( {N – N^{ * } } right) – theta left( {C – C^{ * } } right)} right] + m_{2} left( {V – V^{ * } } right)left[ {frac{ – alpha }{{L_{1} }}left( {V – V^{ * } } right)} right] \ & quad + ,m_{3} left( {G – G^{ * } } right)left[ {beta left( {C – C^{ * } } right) – dleft( {G – G^{ * } } right) + lambda left( {N – N^{ * } } right)} right], \ end{aligned}$$
(16)
$$begin{aligned} frac{dV}{{dt}} & = delta_{1} left( {V – V^{ * } } right)left( {C – C^{ * } } right) + left( {eta – m_{1} theta } right)left( {N – N^{ * } } right)left( {C – C^{ * } } right) – sigma left( {C – C^{ * } } right)^{2} \ & quad – frac{{m_{1} r}}{L}left( {N – N^{ * } } right)^{2} – frac{{m_{2} alpha }}{{L_{1} }}left( {V – V^{ * } } right)^{2} + m_{3} beta left( {G – G^{ * } } right)left( {C – C^{ * } } right) \ & quad + ,m_{3} lambda left( {G – G^{ * } } right)left( {N – N^{ * } } right) – m_{3} dleft( {G – G^{ * } } right)^{2} . \ end{aligned}$$
(17)
Choosing (m_{1} = frac{eta }{theta }), Eq. (17) becomes
$$begin{aligned} frac{dV}{{dt}} & = – frac{sigma }{2}left( {C – C^{ * } } right)^{2} + m_{3} beta left( {G – G^{ * } } right)left( {C – C^{ * } } right) – frac{{m_{3} d}}{2}left( {G – G^{ * } } right)^{2} \ & quad – frac{sigma }{2}left( {C – C^{ * } } right)^{2} + delta_{1} left( {V – V^{ * } } right)left( {C – C^{ * } } right) – frac{{m_{2} alpha }}{{L_{1} }}left( {V – V^{ * } } right)^{2} \ & quad – frac{eta r}{{theta L}}left( {N – N^{ * } } right)^{2} + m_{3} lambda left( {G – G^{ * } } right)left( {N – N^{ * } } right) – frac{{m_{3} d}}{2}left( {G – G^{ * } } right)^{2} . \ end{aligned}$$
(18)
Note that (a_{1} x^{2} + a_{2} xy + a_{3} y^{2}) is negatively defined if (a_{1} < 0) and (a_{2}^{2} < 4a_{1} a_{3}). Using this condition, (frac{dV}{{dt}}) is negatively defined within the region of attraction (Omega) when the provided conditions in Eqs. (10) and (11) for (m_{2}) and (m_{3}) hold.