solution
\begin{align*} \dot{r} & =r\left ( 1-r\right ) \\ \dot{\theta } & =1 \end{align*}
The Poincare section (line) is taken as the \(\theta =0\) axis (this means the x-axis in the Phase plane). The Poincare section can be any other line which is orthogonal to the orbits, but this line is the easiest to use because \(\theta =0\) on it.
Let \(\phi \left ( r_{0}\right ) \) be the Poincare map. This means given an initial conditions point \(r=r_{0}\) on the Poincare section then \(\phi \left ( r_{0}\right ) \) will return the location of the solution \(r\left ( t\right ) \) on this line after the solution \(r\left ( t\right ) \) has run for one period \(T\) of time. i.e. \(r\left ( t+T\right ) \), where the period is \(T\). The following diagram helps illustrate this.
The first thing to determine when finding the Poincare map \(\phi \left ( r_{0}\right ) \) is to solve the ODE. In this case, only the \(r\left ( t\right ) \) ODE needs to be solved because there is no dependency on \(\theta \left ( t\right ) \) in the \(r\left ( t\right ) \) ODE. The following solves the ODE \(\dot{r}=r\left ( 1-r\right ) \).
\begin{align*} \frac{dr}{dt} & =r\left ( 1-r\right ) \\ \frac{dr}{r\left ( 1-r\right ) } & =dt\\ \left ( \frac{1}{r}-\frac{1}{r-1}\right ) dr & =dt \end{align*}
Integrating both sides gives\begin{align*} \int \left ( \frac{1}{r}-\frac{1}{r-1}\right ) dr & =\int dt\\ \ln \left \vert r\right \vert -\ln \left \vert r-1\right \vert & =t+C_{1}\\ \ln \left \vert \frac{r}{r-1}\right \vert & =t+C_{1}\\ \frac{r}{r-1} & =C_{2}e^{t} \end{align*}
At initial conditions, \(t=0,r=r_{0}\) the above becomes\[ \frac{r_{0}}{r_{0}-1}=C_{2}\] Hence the solution is\begin{align*} \frac{r}{r-1} & =\frac{r_{0}}{r_{0}-1}e^{t}\\ r & =\left ( r-1\right ) \left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) \\ r & =r\left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) -\left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) \\ r-r\left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) & =-\left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) \\ r\left ( 1-\frac{r_{0}}{r_{0}-1}e^{t}\right ) & =-\left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) \end{align*}
Simplifying gives\begin{align*} r\left ( t\right ) & =\frac{-\left ( \frac{r_{0}}{r_{0}-1}e^{t}\right ) }{1-\frac{r_{0}}{r_{0}-1}e^{t}}\\ & =\frac{-1}{\frac{1}{\frac{r_{0}}{r_{0}-1}e^{t}}-1}\\ & =\frac{-1}{\frac{\left ( r_{0}-1\right ) e^{-t}}{r_{0}}-1}\\ & =\frac{-r_{0}}{\left ( r_{0}-1\right ) e^{-t}-r_{0}}\\ & =\frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-t}} \end{align*}
At \(t=T\) the above solution becomes\[ r\left ( T\right ) =\frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-T}}\] Therefore the Poincare map is \[ \phi \left ( r_{0}\right ) =\frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-T}}\] Taking the period as \(2\pi \), then the above becomes\begin{equation} \phi \left ( r_{0}\right ) =\frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-2\pi }} \tag{1} \end{equation} What the above function does, is that given a specific value of \(r_{0}\) on the Poincare section (the x-axis), it returns where the solution location will be on this line after one period.
For an example, taking initial conditions \(r_{0}=1.5,\theta =0\) then (1) gives\begin{align*} \phi \left ( 1.5\right ) & =\frac{1.5}{1.5-\left ( 1.5-1\right ) e^{-2\pi }}\\ & =1.00062 \end{align*}
This shows the solution \(r\left ( t\right ) \) after one period, has reached the Poincare section again at a point which is closer to the limit cycle at \(r=1\). Initially \(r\left ( 0\right ) =1.5\) and after \(2\pi \) seconds, \(r\left ( 2\pi \right ) =1.00062\).
Let \(\phi \left ( r_{n}\right ) \) be Poincare map at iteration \(n.\)\begin{align*} \phi \left ( r_{0}\right ) & =\frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-2\pi }}\\ \phi \left ( r_{1}\right ) & =\frac{\left ( \frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-2\pi }}\right ) }{\left ( \frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-2\pi }}\right ) -\left ( \left ( \frac{r_{0}}{r_{0}-\left ( r_{0}-1\right ) e^{-2\pi }}\right ) -1\right ) e^{-2\pi }}=\frac{r_{0}}{r_{0}+e^{-4\pi }-r_{0}e^{-4\pi }}\\ \phi \left ( r_{2}\right ) & =\frac{\left ( \frac{r_{0}}{r_{0}+e^{-4\pi }-r_{0}e^{-4\pi }}\right ) }{\left ( \frac{r_{0}}{r_{0}+e^{-4\pi }-r_{0}e^{-4\pi }}\right ) -\left ( \left ( \frac{r_{0}}{r_{0}+e^{-4\pi }-r_{0}e^{-4\pi }}\right ) -1\right ) e^{-2\pi }}=\frac{r_{0}}{r_{0}+e^{-6\pi }-r_{0}e^{-6\pi }}\\ \phi \left ( r_{3}\right ) & =\frac{\frac{r_{0}}{r_{0}+e^{-6\pi }-r_{0}e^{-6\pi }}}{\left ( \frac{r_{0}}{r_{0}+e^{-6\pi }-r_{0}e^{-6\pi }}\right ) -\left ( \left ( \frac{r_{0}}{r_{0}+e^{-6\pi }-r_{0}e^{-6\pi }}\right ) -1\right ) e^{-2\pi }}=\frac{r_{0}}{r_{0}+e^{-8\pi }-r_{0}e^{-8\pi }}\\ & \vdots \\ \phi \left ( r_{n}\right ) & =\frac{r_{0}}{r_{0}+e^{-\left ( 2+2n\right ) \pi }-r_{0}e^{-\left ( 2+2n\right ) \pi }} \end{align*}
As \(n\rightarrow \infty \) the above gives\[ \lim _{n\rightarrow \infty }\phi \left ( r_{n}\right ) \rightarrow \frac{r_{0}}{r_{0}}=1 \] Therefore \(\phi \left ( r_{0}\right ) \) is stable. The limit cycle at \(r=1\) is stable because any solution that starts always from \(r=1\) will eventually reach \(r=1\). The following also is a plot of \(\phi \left ( r_{0}\right ) \)
solution
The SIRS model is \begin{align} \dot{S} & =-\beta SI+\mu R\tag{1}\\ \dot{I} & =\beta SI-\upsilon I\nonumber \\ \dot{R} & =\upsilon I-\mu R\nonumber \end{align}
Where \(S=S(t)\) is the population of susceptible individuals, \(I=I(t)\), the infected population, and \(R=R(t)\) the recovered population. This diagram shows the model where now some of the recovered population can become susceptible again and later become infected. The parameter \(\mu \) indicates how much of the recovered population could become susceptible again.
The units of \(S,R,I\) are population measured in person. These values can not be negative since they are population amount. The units of \(\mu \) is \(\frac{1}{\text{time}}\) where time can be day or week or any other unit of time. The units of \(\upsilon \) is also \(\frac{1}{\text{time}}\). The units of \(\beta \) is \(\frac{1}{\left ( \text{time}\right ) \left ( \text{person}\right ) }\)
Since \(\frac{d}{dt}\left ( S+I+R\right ) =0\) then \(S+I+R=\tau \), where \(\tau \) is constant (the total size of the population). Therefore \(R=\tau -I-S\). This implies that, since \(\tau \) is known, then finding \(I\) and \(S\) allows finding \(R\) automatically from this algebraic equation. Replacing \(R\) by \(\tau -I-S\) in the first and the second equations in (1) eliminates \(R\) and gives (the restricted system) as\begin{align} \dot{S} & =-\beta SI+\mu \left ( \tau -I-S\right ) \tag{2}\\ \dot{I} & =\beta SI-\upsilon I\nonumber \end{align}
This is the new system of equations. The third ODE is not needed. Once \(S\left ( t\right ) ,I\left ( t\right ) \,\) are solved for, then \(R\left ( t\right ) \) is found from \(R=\tau -I-S\) without the need to solve a third ODE.
From (2), the critical points are solutions to restricted system\begin{align} -\beta SI+\mu \left ( \tau -I-S\right ) & =0\tag{2}\\ \beta SI-\upsilon I & =0 \tag{4} \end{align}
From (4) \(I\left ( \beta S-\upsilon \right ) =0\). This gives solutions \(I=0,S=\frac{\upsilon }{\beta }\). When \(I=0\) then (2) gives \(\mu \left ( \tau -S\right ) =0\) or \(S=\tau \) since \(\mu \) is not zero. The first critical point is therefore \(\left \{ I=0,S=\tau \right \} \).
Now, when \(S=\frac{\upsilon }{\beta }\) Eq. (2) gives \(-\beta \frac{\upsilon }{\beta }I+\mu \left ( \tau -I-\frac{\upsilon }{\beta }\right ) =0\) or \(-\beta \frac{\upsilon }{\beta }I+\mu \tau -\mu I-\mu \frac{\upsilon }{\beta }=0\), hence \begin{align*} I\left ( -\beta \frac{\upsilon }{\beta }-\mu \right ) & =-\mu \tau +\mu \frac{\upsilon }{\beta }\\ I & =-\frac{\mu \left ( \frac{\upsilon }{\beta }-\tau \right ) }{\beta \frac{\upsilon }{\beta }+\mu }\\ & =\frac{\mu \left ( \tau -\frac{\upsilon }{\beta }\right ) }{\beta \frac{\upsilon }{\beta }+\mu } \end{align*}
Therefore the second critical point is \(\left \{ I=\frac{\mu \left ( \tau -\frac{\upsilon }{\beta }\right ) }{\beta \frac{\upsilon }{\beta }+\mu },S=\frac{\upsilon }{\beta }\right \} \) or \(\left \{ I=\frac{\tau \mu \beta -\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta }\right \} \). These critical points are for the reduced system.
For each one of these critical points found above, the critical value for \(R\) can now be found from the relation \(R=\tau -I-S\). For the first critical point \(\left \{ I=0,S=\tau \right \} \) then \(R=\tau -0-\tau =0\). Therefore this critical point becomes \(\left \{ I=0,S=\tau ,R=0\right \} \).
For the second critical point found above then \begin{align*} R & =\tau -I-S\\ & =\tau -\frac{\tau \mu \beta -\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) }-\frac{\upsilon }{\beta }\\ & =\upsilon \frac{\beta \tau -\upsilon }{\beta \left ( \mu +\upsilon \right ) } \end{align*}
Hence this critical point becomes \(\left \{ I=\frac{\tau \mu \beta -\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \frac{\beta \tau -\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right \} \). Therefore the two critical points for the full system are\begin{align*} & \left \{ I=0,S=\tau ,R=0\right \} \\ & \left \{ I=\frac{\tau \mu \beta -\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \frac{\beta \tau -\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right \} \end{align*}
case (a) Under the case that \(\tau >\frac{\upsilon }{\beta }\) then the second critical point above becomes\[ I=\frac{\tau \mu \beta }{\beta \left ( \upsilon +\mu \right ) }-\frac{\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \left ( \frac{\tau }{\left ( \mu +\upsilon \right ) }-\frac{\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right ) \] Replacing \(\tau \) by \(\frac{\upsilon }{\beta }+\delta \) where \(\delta >0\) then the above becomes\begin{align*} I & =\frac{\upsilon }{\beta }\frac{\mu }{\left ( \upsilon +\mu \right ) }+\delta \frac{\mu }{\left ( \upsilon +\mu \right ) }-\frac{\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \left ( \frac{\upsilon }{\beta }\frac{1}{\left ( \mu +\upsilon \right ) }+\frac{\delta }{\mu +\upsilon }-\frac{\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right ) \\ I & =\delta \frac{\mu }{\left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\frac{\upsilon \delta }{\mu +\upsilon } \end{align*}
Since all \(\mu ,\upsilon ,\beta ,\delta >0\) then the above is a valid critical point. Hence under this case the critical points are\begin{align*} & \left \{ I=0,S=\tau ,R=0\right \} \\ & \left \{ I=\delta \frac{\mu }{\left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\frac{\upsilon \delta }{\mu +\upsilon }\right \} \end{align*}
Where \(\tau =\frac{\upsilon }{\beta }+\delta \) where \(\delta >0\,.\)
case (b) Under the case that \(\tau <\frac{\upsilon }{\beta }\) the second critical point above becomes\[ I=\frac{\tau \mu \beta }{\beta \left ( \upsilon +\mu \right ) }-\frac{\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \left ( \frac{\tau }{\left ( \mu +\upsilon \right ) }-\frac{\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right ) \] Replacing \(\tau \) by \(\frac{\upsilon }{\beta }-\delta \) where \(\delta >0\) the above becomes\begin{align*} I & =\frac{\upsilon }{\beta }\frac{\mu }{\left ( \upsilon +\mu \right ) }-\delta \frac{\mu }{\left ( \upsilon +\mu \right ) }-\frac{\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \left ( \frac{\upsilon }{\beta }\frac{1}{\left ( \mu +\upsilon \right ) }-\frac{\delta }{\mu +\upsilon }-\frac{\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right ) \\ I & =-\delta \frac{\mu }{\left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=-\frac{\upsilon \delta }{\mu +\upsilon } \end{align*}
Since all \(\mu ,\upsilon ,\beta ,\delta >0\) then the above is a not a valid critical point, because it says \(I,R\) are negative and this is not possible, as these are population size. Hence under this case the critical point is only\[ \left \{ I=0,S=\tau ,R=0\right \} \] case \(\tau =\frac{\upsilon }{\beta }\)
Under this case then the second critical point above becomes\[ I=\frac{\tau \mu \beta }{\beta \left ( \upsilon +\mu \right ) }-\frac{\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \left ( \frac{\tau }{\left ( \mu +\upsilon \right ) }-\frac{\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right ) \] Replacing \(\tau \) by \(\frac{\upsilon }{\beta }\) then the above becomes\begin{align*} I & =\frac{\upsilon }{\beta }\frac{\mu }{\left ( \upsilon +\mu \right ) }-\frac{\upsilon \mu }{\beta \left ( \upsilon +\mu \right ) },S=\frac{\upsilon }{\beta },R=\upsilon \left ( \frac{\upsilon }{\beta }\frac{1}{\left ( \mu +\upsilon \right ) }-\frac{\upsilon }{\beta \left ( \mu +\upsilon \right ) }\right ) \\ I & =0,S=\frac{\upsilon }{\beta },R=0 \end{align*}
But this is the same as the first critical point. Hence only one critical point exist under this condition. This is summary of result.
| case | critical points |
| case (a) \(\tau \geq \frac{\upsilon }{\beta }\) | \(\left \{ I=0,S=\tau ,R=0\right \} ,\left \{ I=\frac{\mu }{\upsilon +\mu }\delta ,S=\frac{\upsilon }{\beta },R=\frac{\upsilon }{\mu +\upsilon }\delta \right \} \) where \(\delta =\tau -\frac{\upsilon }{\beta }\) |
| case (b) \(\tau \leq \frac{\upsilon }{\beta }\) | \(\left \{ I=0,S=\tau ,R=0\right \} \) |
The full system (1) is (I assumed we need to find eigenvalues for full system. The question was not clear on this).\begin{align} \dot{S} & =-\beta SI+\mu R\tag{1}\\ \dot{I} & =\beta SI-\upsilon I\nonumber \\ \dot{R} & =\upsilon I-\mu R\nonumber \end{align}
The Jacobian matrix is \[ J=\begin{pmatrix} \frac{\partial \dot{S}}{\partial S} & \frac{\partial \dot{S}}{\partial I} & \frac{\partial \dot{S}}{\partial R}\\ \frac{\partial \dot{I}}{\partial S} & \frac{\partial \dot{I}}{\partial I} & \frac{\partial \dot{I}}{\partial R}\\ \frac{\partial \dot{R}}{\partial S} & \frac{\partial \dot{R}}{\partial I} & \frac{\partial \dot{R}}{\partial R}\end{pmatrix} =\begin{pmatrix} -\beta I & -\beta S & \mu \\ \beta I & \beta S-\upsilon & 0\\ 0 & \upsilon & -\mu \end{pmatrix} \]
Case a At the first critical point \(\left \{ I=0,S=\tau ,R=0\right \} \) the above becomes\[ A=\begin{pmatrix} 0 & -\beta \tau & \mu \\ 0 & \beta \tau -\upsilon & 0\\ 0 & \upsilon & -\mu \end{pmatrix} \] The eigenvalues are from \(\left \vert A-\lambda I\right \vert =0\) or\begin{align*} \begin{vmatrix} -\lambda & -\beta \tau & \mu \\ 0 & \beta \tau -\upsilon -\lambda & 0\\ 0 & \upsilon & -\mu -\lambda \end{vmatrix} & =0\\ -\lambda \begin{vmatrix} \beta \tau -\upsilon -\lambda & 0\\ \upsilon & -\mu -\lambda \end{vmatrix} +\beta \tau \begin{vmatrix} 0 & 0\\ 0 & -\mu -\lambda \end{vmatrix} +\mu \begin{vmatrix} 0 & \beta \tau -\upsilon -\lambda \\ 0 & \upsilon \end{vmatrix} & =0\\ -\lambda \left ( \beta \tau -\upsilon -\lambda \right ) \left ( -\mu -\lambda \right ) & =0\\ -\lambda \left ( \lambda +\mu \right ) \left ( \lambda +\upsilon -\beta \tau \right ) & =0 \end{align*}
Hence \(\lambda =0,\lambda =-\mu ,\lambda =\beta \tau -\upsilon \) are the eigenvalues.
At the second critical point \(\left \{ I=\frac{\mu }{\upsilon +\mu }\delta ,S=\frac{\upsilon }{\beta },R=\frac{\upsilon }{\mu +\upsilon }\delta \right \} \) the Jacobian matrix becomes\begin{align*} A & =\begin{pmatrix} -\beta I & -\beta S & \mu \\ \beta I & \beta S-\upsilon & 0\\ 0 & \upsilon & -\mu \end{pmatrix} \\ & =\begin{pmatrix} -\beta \frac{\mu }{\upsilon +\mu }\delta & -\beta \frac{\upsilon }{\beta } & \mu \\ \beta \frac{\mu }{\upsilon +\mu }\delta & \beta \frac{\upsilon }{\beta }-\upsilon & 0\\ 0 & \upsilon & -\mu \end{pmatrix} \\ & =\begin{pmatrix} -\beta \frac{\mu }{\upsilon +\mu }\delta & -\upsilon & \mu \\ \beta \frac{\mu }{\upsilon +\mu }\delta & 0 & 0\\ 0 & \upsilon & -\mu \end{pmatrix} \end{align*}
The eigenvalues are from \(\left \vert A-\lambda I\right \vert =0\) or\begin{align*} \begin{vmatrix} -\beta \frac{\mu }{\upsilon +\mu }\delta -\lambda & -\upsilon & \mu \\ \beta \frac{\mu }{\upsilon +\mu }\delta & -\lambda & 0\\ 0 & \upsilon & -\mu -\lambda \end{vmatrix} & =0\\ \left ( -\beta \frac{\mu }{\upsilon +\mu }\delta -\lambda \right ) \begin{vmatrix} -\lambda & 0\\ \upsilon & -\mu -\lambda \end{vmatrix} +\upsilon \begin{vmatrix} \frac{\tau \mu \beta -\upsilon \mu }{\left ( \upsilon +\mu \right ) } & 0\\ 0 & -\mu -\lambda \end{vmatrix} +\mu \begin{vmatrix} \beta \frac{\mu }{\upsilon +\mu }\delta & -\lambda \\ 0 & \upsilon \end{vmatrix} & =0\\ \left ( -\beta \frac{\mu }{\upsilon +\mu }\delta -\lambda \right ) \left ( -\lambda \right ) \left ( -\mu -\lambda \right ) +\upsilon \left ( \frac{\tau \mu \beta -\upsilon \mu }{\left ( \upsilon +\mu \right ) }\right ) \left ( -\mu -\lambda \right ) +\mu \upsilon \left ( \beta \frac{\mu }{\upsilon +\mu }\delta \right ) & =0\\ \lambda ^{3}+\lambda ^{2}\left ( \beta \frac{\mu }{\upsilon +\mu }\delta +\mu \right ) +\lambda \left ( \beta \frac{\mu \upsilon }{\upsilon +\mu }\delta +\beta \frac{\mu ^{2}}{\upsilon +\mu }\delta \right ) & =0\\ \lambda \left ( \lambda ^{2}+\lambda \left ( \beta \frac{\mu }{\upsilon +\mu }\delta +\mu \right ) +\left ( \beta \frac{\mu \upsilon }{\upsilon +\mu }\delta +\beta \frac{\mu ^{2}}{\upsilon +\mu }\delta \right ) \right ) & =0 \end{align*}
Hence one eigenvalue is \(\lambda _{1}=0\) and the other two are the solution to \[ \lambda ^{2}+\lambda \left ( \beta \delta \frac{\mu }{\upsilon +\mu }+\mu \right ) +\frac{\beta \delta \mu }{\upsilon +\mu }\left ( \mu +\upsilon \right ) =0 \] Which is quadratic. Using \(\lambda =\frac{-b}{2a}\pm \frac{1}{2a}\sqrt{b^{2}-4ac}\) then the second and third eigenvalues are\[ \lambda _{2,3}=-\frac{\beta \delta \frac{\mu }{\upsilon +\mu }+\mu }{2}\pm \frac{1}{2}\sqrt{\left ( \beta \delta \frac{\mu }{\upsilon +\mu }+\mu \right ) ^{2}-4\left ( \frac{\beta \delta \mu }{\upsilon +\mu }\left ( \mu +\upsilon \right ) \right ) }\]
Case b Only one critical point \(\left \{ I=0,S=\tau ,R=0\right \} \). This was found above as \(\lambda =0,\lambda =-\mu ,\lambda =\beta \tau -\upsilon \).
\(\allowbreak \)Summary of result
In the following, for case a, \(\tau -\frac{\upsilon }{\beta }=\delta \) where \(\delta >0\)
| case | result | |||||||||||||||
| case (a) \(\tau \geq \frac{\upsilon }{\beta }\) |
|
|||||||||||||||
| case (b) \(\tau \leq \frac{\upsilon }{\beta }\) |
|
|||||||||||||||
For case (a), the critical points are found to be \(\left \{ I=0,S=\tau ,R=0\right \} \) and \(\left \{ I=\frac{\mu }{\upsilon +\mu }\delta ,S=\frac{\upsilon }{\beta },R=\frac{\upsilon }{\mu +\upsilon }\delta \right \} \) where \(\tau -\frac{\upsilon }{\beta }=\delta \) and \(\delta >0\).
For the first critical point, it says that there are no infected population and no recovered population. Hence we expect this to be stable, since no infection will occur, and the population \(\tau \) is all of type \(S\left ( t\right ) \) and will remain so for all time.
For the second critical point \(\left \{ I=\frac{\mu }{\upsilon +\mu }\delta ,S=\frac{\upsilon }{\beta },R=\frac{\upsilon }{\upsilon +\mu }\delta \right \} \), it implies that amount of people who become infected each day is offset by the amount of people who are recovered but become again susceptible. And the amount of people who recover each day is offset by the same amount of people who were recovered but become susceptible. This keeps the size of infected people each day the same not changing, as well as the size of susceptible and recovered population the same. This will only happen when \(\tau \geq \frac{\upsilon }{\beta }\). A simulation is made to show this. This below shows, using \(\tau =1000\), the second critical point when \(\frac{\upsilon }{\beta }=707\), showing the current size of \(S\left ( t\right ) ,I\left ( t\right ) ,R\left ( t\right ) \) which as be seen, remain the same for all time.
This shows what happens when the infection rate \(\beta \) is relatively high at \(\beta =0.0013\)
This shows what happens when the infection rate \(\beta \) has been decreased down to \(\beta =0.0003\). Notice the flattening of the red curve (infected population).
Considering now only the restricted system found in part 1, which is\begin{align*} \dot{S} & =-\beta SI+\mu \left ( \tau -I-S\right ) \\ \dot{I} & =\beta SI-\upsilon I \end{align*}
The condition \(S+I\leq \tau \) is always true since total size of population is fixed \(\tau =S+I+R\).
For condition \(I,S\geq 0\), it means that there are some infected population and some susceptible population present. This is of interest, since it implies dynamics of infection is still in place and taking effect.
Going back now to the full model\begin{align*} \dot{S} & =-\beta SI+\mu R\\ \dot{I} & =\beta SI-\upsilon I\\ \dot{R} & =\upsilon I-\mu R \end{align*}
Considering the line \(I\left ( t\right ) \). To check if it is invariant of not. Then assuming initial conditions starts on this line only, which means \(S\left ( 0\right ) =0,R\left ( 0\right ) =0,I\left ( 0\right ) =\tau \). When this is the case, then we see that from third equation that \(R\left ( t\right ) \) will change with time, since it depends on \(I\left ( t\right ) \) which is not zero. This in turn causes \(S\left ( t\right ) \) to change from zero. Therefore the line \(I\left ( t\right ) \) is not invariant. In other words, if \(I\left ( 0\right ) =\tau \), then solution will not remain on this line for all time and eventually both \(R\left ( t\right ) \) and \(S\left ( t\right ) \) will be non-zero.
To verify this, I run the simulation program with \(I\left ( 0\right ) =\tau =1000\) and it shows that solution will not remain on axis \(I\left ( t\right ) \). This is because some infected population will recover, and some of those who recovered will become susceptible again.
Now let us consider the \(S\left ( t\right ) \) axis. This means \(S\left ( 0\right ) =\tau \) and now we check if the solution remains on \(S\left ( t\right ) \) axis or not. In this case we would expect \(S\) axis to be invariant, since there are no infected population. This means \(I\left ( 0\right ) =0,R\left ( 0\right ) =0\). From (1) this shows that \(\dot{R}\left ( 0\right ) =0\), hence \(R\left ( t\right ) \) is constant and remain zero for all time. And since \(I\left ( 0\right ) =0\) then second equation in (1) also shows that \(\dot{I}\left ( 0\right ) =0\), hence \(I\left ( t\right ) \) is constant and remain zero for all time. This means first equation of (1) gives \(\dot{S}\left ( t\right ) =0\) and \(S\left ( t\right ) \) is constant for all time which is \(\tau \). So the solution remains on the \(S\) axis. Hence \(S\) axis is invariant.
Considering now only the restricted system found in part 1, which is\begin{align*} \dot{S} & =-\beta SI+\mu \left ( \tau -I-S\right ) \\ \dot{I} & =\beta SI-\upsilon I \end{align*}
Starting in region \(\Delta \), the solution will remain in \(\Delta \), because total population is fixed and there is no death. Hence starting with any combination of initial conditions \(S\left ( 0\right ) ,I\left ( 0\right ) \) the solution will always be in this region. Hence positive invariant.
solution
\(\omega \) limit set, is the set of all points that are the limit of all positive orbits \(\gamma ^{+}\left ( x\right ) \). In other words, given a specific orbit \(\gamma ^{+}\left ( x\right ) \) that starts at some initial conditions point \(x_{0}\) and if as \(t\rightarrow \infty \) this orbit terminates at point \(p\) then \(p\) is in the \(\omega \) limit set of such orbit.
Similarly, the \(\alpha \) limit set, is the set of all points that are the limit of the negative orbits \(\gamma ^{-}\left ( x\right ) \). These are orbits where as \(t\rightarrow -\infty \) the orbit would have originated from the point \(p\). Then \(p\) is in the \(\alpha \) limit set of such orbit.
To find the \(\omega \) limit set, we need to find the points where solutions terminate at them eventually (attractive or saddle points) and the \(\alpha \) limit set are points that solutions do not terminate at them, but move away from them (negative attractions).
\begin{align*} r^{\prime } & =r\left ( 1-r^{2}\right ) \\ \theta ^{\prime } & =1 \end{align*}
The critical points are \(r=0,r=-1,r=1\). But since this is polar coordinates and \(r\) is the radius, then \(r\) can not be negative, Hence only \(r=0,r=1\) are the critical points (with \(\theta \) any value). This implies that all points on circle of radius \(1\) and the point \(r=0\). To determine \(\omega \) limit set and \(\alpha \) limit set\(\,\) we have to check the eigenvalues. The Jacobian matrix is \[ J=\begin{pmatrix} \frac{\partial r^{\prime }}{\partial r} & \frac{\partial r^{\prime }}{\partial \theta }\\ \frac{\partial \theta ^{\prime }}{\partial r} & \frac{\partial \theta ^{\prime }}{\partial \theta }\end{pmatrix} =\begin{pmatrix} 1-3r^{2} & 0\\ 0 & 0 \end{pmatrix} \] At \(r=0\)\[ \left \vert J-\lambda I\right \vert =\begin{vmatrix} 1-\lambda & 0\\ 0 & -\lambda \end{vmatrix} =0 \] Therefore \(\left ( 1-\lambda \right ) \left ( -\lambda \right ) =0\) or \(\lambda =0,\lambda =1\). Hence the point \(r=0\) is unstable. Therefore the \(\alpha \) limit set is the origin \(\left \{ r=0,\theta \right \} \) for all orbits that start with initial conditions \(r\left ( 0\right ) <1\). There is no \(\alpha \) limit set for orbits that start outside \(r>1\) since there is no critical point outside the limit cycle.
At \(r=1\) the determinant becomes \[ \left \vert J-\lambda I\right \vert =\begin{vmatrix} \left ( 1-3\right ) -\lambda & 0\\ 0 & -\lambda \end{vmatrix} =0 \] Therefore \(\left ( -2-\lambda \right ) \left ( -\lambda \right ) =0\). Hence \(\lambda =0,\lambda =-2\). Hence \(r=1\) is stable. What this means is that \(\omega \) limit set is all points on the unit circle corresponding to all orbits that start anywhere, including on the limit circle itself, (except orbits that have initial conditions \(r=0\)). The following diagram illustrates this result
\begin{align*} r^{\prime } & =r\left ( r^{2}-3r+2\right ) \\ \theta ^{\prime } & =1 \end{align*}
The critical points are \(r=0,\) all points on circle of radius \(r=2\) and all points on circle of radius \(r=1\).
To determine \(\omega \) limit set and \(\alpha \) limit set\(\,\) we have to check the eigenvalues at each set of critial points fround above. The Jacobian matrix is \[ J=\begin{pmatrix} \frac{\partial r^{\prime }}{\partial r} & \frac{\partial r^{\prime }}{\partial \theta }\\ \frac{\partial \theta ^{\prime }}{\partial r} & \frac{\partial \theta ^{\prime }}{\partial \theta }\end{pmatrix} =\begin{pmatrix} 3r^{2}-6r+2 & 0\\ 0 & 0 \end{pmatrix} \] At \(r=0\)\[ \left \vert J-\lambda I\right \vert =\begin{vmatrix} 2-\lambda & 0\\ 0 & -\lambda \end{vmatrix} =0 \] Therefore \(\left ( 2-\lambda \right ) \left ( -\lambda \right ) =0\) or \(\lambda =0,\lambda =2\). Hence the point \(r=0\) is unstable.
At \(r=1\)\begin{align*} \left \vert J-\lambda I\right \vert & =\begin{vmatrix} \left ( 3-6+2\right ) -\lambda & 0\\ 0 & -\lambda \end{vmatrix} \\ & =\begin{vmatrix} -1-\lambda & 0\\ 0 & -\lambda \end{vmatrix} =0 \end{align*}
Therefore \(\left ( -1-\lambda \right ) \left ( -\lambda \right ) =0\). Hence \(\lambda =0,\lambda =-1\). Hence all points on circle \(r=1\) are stable critical points.
At \(r=2\),\begin{align*} \left \vert J-\lambda I\right \vert & =\begin{vmatrix} \left ( 12-12+2\right ) -\lambda & 0\\ 0 & -\lambda \end{vmatrix} \\ & =\begin{vmatrix} 2-\lambda & 0\\ 0 & -\lambda \end{vmatrix} =0 \end{align*}
Therefore \(\left ( 2-\lambda \right ) \left ( -\lambda \right ) =0\). Hence \(\lambda =0,\lambda =2\). Hence \(r=2\) is not stable. The following diagram illustrates the result of what was found above showing the \(\omega \) limit set and the \(\alpha \) limit sets found.
