How to Draw Phase Plane of Autonomous Od

phaseR: An R bundle for phase plane analysis of ane and two-dimensional autonomous ODE systems [1]. The phaseR packet uses stability analysis to allocate equilibrium points.

Fixed Points and Stage Portraits

Equilibrium points and stability. Equilibrium points of an autonomous ODE are defined at $ten^*=0$ such that $f(x^*)=0$. Why do we desire to find the equilibrium points? Starting time at a point $x^*$ such that $f(x^*)=0$, the organization, if unperturbed, will remain at $10^*$. Hence, these points determine the long-term behavior of a differential equation.

Phase portrait: A stage portrait is defined as the geometrical representation of the trajectories of the dynamical organization in the phase plane of the system equation. Every set of the initial status is represented by a different curve or signal in the phase aeroplane.

Linearization

Fixed Points $\mathbf{(x^{\star}, y^{\star})}$: occur for $\mathrm{\dot x = 0}$ and $\mathrm{\dot y = 0}$
Substitute $\mathrm{\dot x}$ and $\mathrm{\dot y}$ into the Jacobian Matrix

\[ \Big \mathbf{A} = \begin{bmatrix} \mathrm{\frac{\partial{\dot 10}}{\partial x}} & \mathrm{\frac{\partial \dot x}{\partial y}} \\ \mathrm{\frac{\partial \dot y}{\partial ten}} & \mathrm{\frac{\partial \dot y}{\partial y}} \end{bmatrix} \]

Evaluate matrix $\mathbf{A}$ at whatever fixed bespeak $\mathbf(x^{\star}, y^{\star})$
The Eigen value $\mathbf{\lambda}$ for the fixed point is calculated with the characteristic equation $\mathrm{det}\mathbf{(A-\lambda I)}=0$.

\[ \brainstorm{align} \large \mathbf{A}_{\pocket-size\mathrm{(x^{\star}, y^{\star})}} \mathbf{-} \brainstorm{bmatrix} \lambda & 0 \\ 0 & \lambda \stop{bmatrix} &= 0 & \\ \begin{bmatrix} a_1 - \lambda & a_2 \\ a_3 & a_4 -\lambda \cease{bmatrix} &= 0 & \\ (a_1 - \lambda) \times (a_4 -\lambda) - (a_2) \times (a_3) &= 0 \end{align} \]

The roots of the to a higher place equation can exist calculated every bit

\[\lambda = \mathrm{\frac{-b \pm \sqrt{b^2 - 4ac}}{2a}}\]

Sketch the Phase Portrait

Types of Fixed Points

a. Stable Node

The fixed point volition exist a stable node if both eigen values $\lambda_{1,2}$ are existent numbers with the aforementioned sign: $\lambda_{1,2} > 0$ or $\lambda_{one,2} < 0$.

b. Saddle Node

The fixed point volition exist a saddle betoken if i eigen value $\lambda_1$ is greater than 0, and the other eigen value $\lambda_2$ is less than 0.

c. Eye

The origin of the system will form a center if both eigen values $\lambda_{i,2}$ are pure imaginary.
The stock-still signal will be an outgoing spirals if both eigen values $\lambda_{one,2}$ are circuitous with positive real parts.
The fixed point will be oscillating in nature if the eigen values $\lambda_i$ change from real to imaginary.

d. Isolated

The origin volition be an isolated fixed point if the determinant $\nabla$ comes out to exist 0.

due east. Unstable Star

The fixed indicate will be a unstable star if there is just one eigen value $\lambda_1$ for the system.

Equilibrium points

Instance. Discover the equilibrium points to the post-obit ODE:

\[\frac{dy}{dt}=4-y^two\]

\[ \brainstorm{align} \frac{dy}{dt} &= 4-y^2 = 0 \\ 0 & = (ii-y)(2+y) \\ \Rightarrow& y = -2,2 \end{align} \]

Stability of the equilibrium points

Definition. if for every $\epsilon>0$, in that location exists $\delta>0$ such that whenever $|y(0)-y_*|<\delta$ then $|y(t)-y_*|<\epsilon$ for all $t$. Hence, a fixed betoken is stable if a system placed a small distance away from the fixed point continues to remain close to the fixed bespeak. And a fixed point is unstable if a system placed a small perturbation away from the fixed signal causes the solution to diverge.

There are 2 methods for determining the stability of a fixed signal: Stage Portrait Analysis or the Taylor Series Expansion.

Phase Portrait Analysis

Graphical Interpretation. A stage portrait plots the derivative $\dot x$ against the dependent variable $x$. We tin find the equilibrium points at locations where $f$ crosses the $\mathcal x-\text{axis}$.

We can represent the menstruum of $f$ in the phase portrait by placing arrows along the dependent variable's axis, indicating whether $f$ would be increasing or decreasing. Points where arrows on both side of the equilibrium point towards each other $\rightarrow\;\; \leftarrow$ denotes stability. And points where arrows on both side of the equilibrium point away from each other $\leftarrow\;\; \rightarrow$ denotes instability.

In the following, we plot the phase portrait for $\frac{dy}{dt}=4-y^2$. The trajectories plotted shows that solutions converge towards $y=2$, only away from $y=-ii$. Hence, the equilibrium point $y^*=2$ is stable and the equilibrium point $y^*=-two$ is unstable.

$Figure 1: Phase Portrait for $\frac{dy}{dt}=4-y^2$. The trajectories plotted shows that solutions converge towards $y=2$, but away from $y=-2$.$

Figure 1: Phase Portrait for $\frac{dy}{dt}=4-y^2$. The trajectories plotted shows that solutions converge towards $y=ii$, merely abroad from $y=-2$.

Uniqueness Theorem. The above method assumes $f$ to be continuous and differentiable. Hence, these weather guarantee just unique solutions to the democratic differential equation. Therefore, the solution curves cannot touch, wait when they converge at equilibrium points.

Taylor Series Expansion

The second method for performing stability analyses utilizes the Taylor Series expansion of $f$.

Assumptions: Suppose we are a minor altitude $\delta(0)$ away from fixed point $y_*$. Then $y(0)=y_*+\delta(0)$ and $y(t)=y_*+\delta(t)$. And then, the Taylor Series of $f$ can exist written as the post-obit, such that the $y^*$ input represents the point where we perform stability assay:

\[f(y^*+\delta)=f(y^*)+\delta\frac{\partial f}{\partial y}(y^*)+o(\delta)\]

Example 1. 1-Dimesnional ODE

Consider the one-dimensional democratic ODE: \[\frac{dy}{dt}=y(ane-y)(2-y)\]

Menstruation Field

The following plots the flow field and various trajectories, adding horizontal lines at equilibrium points:

Figure 2: The flow field and various trajectories, adding horizontal lines at equilibrium points.

Figure 2: The menstruum field and diverse trajectories, calculation horizontal lines at equilibrium points.

Fixed Points

The horizontal lines on the graph indicate that three equilibrium points take been identified at $y^*=0,i,ii$. If nosotros set $\chapeau y=0,$ nosotros can analytically solve for the three equilibrium points:

\[ \begin{align} y^* (i-y^* )(2-y^* ) &= 0 \\ y^* &= 0, 1, two \finish{marshal} \]

Stability of Fixed Points

Method 1. Phase Portrait

Plotting the phase portrait, we find that $y^*=0$ and $y^*=ii$ are unstable; and $y^*=1$ is stable

Figure 3: The flow field and various trajectories, adding horizontal lines at equilibrium points.

Method 2. Taylor Series Approach

\[ \begin{align} \frac{dy}{dt}&=y (1-y )(ii-y ) \\ &= y^iii-4y^two+2y \terminate{align} \]

Using the Taylor Series approach we have:

\[ \begin{align} \frac{d}{dy}\left.\left(\frac{dy}{dt}\right)\right|_{y=y^*} &= 3y^{*^two} - 6y^* + 2 = \begin{cases} ii, & \ y^*=0,\\ -1 ,& \ y^*=one,\\ 2, & \ y^*=ii. \end{cases} \terminate{align} \]

We draw the same conclusion as from the phase portrait. We can ostend the Taylor assay using stability() to check the stability of each equilibrium signal:

Hence nosotros attain conclusion as above that $y^*=two$ is stable, and $y^*=-ii$ is unstable. Therefore, if $y(0)>2$ or $0<y(0)<ii$, then the solution will eventually approach $y=two$. However, if $y(0)<0$, $y\rightarrow-\infty$ equally $t\rightarrow\infty$.

Case 2. Two-Dimensional ODE

Consider the Lotka-Volterra model, a simple two species competition model, used in ecology, given past the following:

\[\frac{dx}{dt} = 10-xy, \ \frac{dy}{dt} = xy-y\]

a. Plot the velocity field with several trajectories:

$Figure 4: Plot of the velocity field with several trajectories for $\frac{dx}{dt} = x-xy, \frac{dy}{dt} = xy-y$.$

Figure iv: Plot of the velocity field with several trajectories for $\frac{dx}{dt} = x-xy, \frac{dy}{dt} = xy-y$.

Nullclines

Here, $x$-nullclines are defined where $f(10,y)=0$, while the $y$-nullclines are defined where $chiliad(10,y)=0$. Thus, the $x$- and $y$-nullclines define the locations where $10$ and $y$ practice not change with time $t$. When plotting a vector field, it is a good idea to plot the nullclines commencement, considering the line segments/vectors along the nullclines move parallel to the $ten$- and $y$-axes.

b. Calculate and Sketch the Nullclines:

\[ \begin{align} \mathbf{\dot x = 0} \\ & \mathrm{x-xy = 0} \ \Longrightarrow \ x(1-y) = 0 \\ & \Rightarrow \mathbf{x=0} \text{ or } \mathbf{y=1},\\ \mathbf{\dot y = 0}\\ & \mathrm{xy-y = 0} \Rightarrow y(ten-1) = 0 \\ & \Rightarrow \mathbf{x=1} \text{ or } \mathbf{y=0} \end{align} \]

Hence, the fixed points come up out to be $\mathrm{(0, 0)}$ and $\mathrm{(1, ane)}$.

$Figure 5: Sketch of the nullclines for the system of equations $\frac{dx}{dt} = x-xy, \frac{dy}{dt} = xy-y$.$

Figure five: Sketch of the nullclines for the system of equations $\frac{dx}{dt} = x-xy, \frac{dy}{dt} = xy-y$.

Equilibrium points and stability

Equilibrium points defined as the Fixed points are at $(x_*,y_*)$ where:

\[ f(x^*,y^*)=g(x^*,y^*)=0 \]

c. Taking fractional derivatives we compute the Jacobian at any equilibrium point $(x^∗,y^∗)$:

\[ \brainstorm{marshal} \mathbf{A} &= \brainstorm{bmatrix} \frac{\partial (x-xy)}{\fractional x} & \frac{\partial (x-xy)}{\partial y} \\ \frac{\partial (xy-y)}{\partial x} & \frac{\fractional (xy-y)}{\fractional y} \end{bmatrix} \\ \\ \mathbf{A}_{(ten,y)} &= \begin{pmatrix} \mathrm{1-y} & \mathrm{-x} \\ \mathrm{y} & \mathrm{x-1} \end{pmatrix} \end{align} \]

For the fixed betoken $\mathbf{(0,0)}$

\[ \mathbf{A}_{(0,0)} = \left.\begin{pmatrix} \mathrm{1} & \mathrm{0} \\ \mathrm{0} & \mathrm{-1} \finish{pmatrix}\right|_{(0,0)} \]

Then $\text{tr}(J)=T=0$ and $\text{det}(J)=\Delta=-1$; which from our table to a higher place makes $(0,0)$ a saddle point.

For the fixed signal $\mathbf{(ane,i)}$

\[ \mathbf{A}_{(1,1)} = \left. \begin{pmatrix} \mathrm{0} & \mathrm{-ane} \\ \mathrm{1} & \mathrm{0} \cease{pmatrix}\right|_{(one,1)} \]

Therefore, $\text{tr}(J)=T=0$ and $\text{det}(J)=\Delta=i$; which from our table above makes $(i,one)$ a centre.

If we look back at our earlier plot, nosotros can observe trajectories diverging away from $(0,0)$, simply traversing effectually $(ane,1)$.

d. Plot the Oscillating Nature

Here nosotros plot $x$ and $y$ trajectories against $t$. For the case of $(x_0,y_0)=(3,four)$ in this example organization this results in the following plot which we can view the aquiver nature of $x$ and $y$:

Figure 6: Plot $x$ and $y$ trajectories against $t$. In this example, we can view the oscillating nature of $x$ and $y$.

Figure 6: Plot $x$ and $y$ trajectories against $t$. In this example, we can view the aquiver nature of $ten$ and $y$.

Case three. Wilson-Cowan Model

The Wilson-Cowan System is a coupled, nonlinear, differential equation for the excitatory and inhibitory populations' firing rates of neurons:

\[\tau_e \frac{dE}{dt} = -E + (k_e - r_e East) \, S_e(c_1 E - c_2 I + P)\]

\[\tau_i \frac{dI}{dt} = -I + (k_i - r_i I) \, S_i(c_3 E - c_4 I + Q),\]

In our numerical search of the steady-state points, we study the cases where the derivatives of $Eastward$ and $I$ are zero. Since the system is highly nonlinear, we have to do it numerically. Kickoff, we describe the menstruum field so draw the nullclines of the organization. The $x$-nullclines are defined by $f(10,y)=0$, and the $y$-nullclines are defined by $yard(x, y)=0$. These are the locations where $x$ and $y$ exercise not change with time.

Limit Cycles. Non-linear systems tin can exhibit a type of behavior known as a limit bike. If there is only one steady-state solution, and if the steady-land solution is unstable, a limit cycle will occur. In the following, we define the parameters for satisfying conditions 18 and 20 as $c_1=xvi$, $c_2 = 12$, $c_3=fifteen$, $c_4=iii$, $a_e = 1.3$, $\theta_e=4$, $a_i=2$, $\theta_i = 3.7$, $r_e=1$ and $r_i=1$. Nosotros can make up one's mind a steady-land solution by the intersection of the nullclines equally follows.

$Figure 7: Phase Plane Analysis. Determine the steady-state solution by the nullclines' intersection. Parameters: $c_1=16$, $c_2 = 12$, $c_3=15$, $c_4=3$, $a_e = 1.3$, $\theta_e=4$, $a_i=2$, $\theta_i = 3.7$, $r_e=1$, $r_i=1$.$

Figure 7: Phase Plane Analysis. Make up one's mind the steady-state solution by the nullclines' intersection. Parameters: $c_1=16$, $c_2 = 12$, $c_3=xv$, $c_4=3$, $a_e = one.3$, $\theta_e=4$, $a_i=2$, $\theta_i = 3.7$, $r_e=one$, $r_i=1$.

In the phase airplane, the limit cycle is an isolated closed orbit, where "closed" means the periodicity of move, and "isolated" means the limit of motility, where nearby trajectories converge or deviate. We can change our initial values of $E_0$ and $I_0$ to obtain dissimilar paths in the phase space.

$Figure 8: Phase Plane Analysis showing limit cycle trajectory in response to constant simulation $P=1.25$. Dashed lines are nullclines. Parameters: $c_1=16$, $c_2 = 12$, $c_3=15$, $c_4=3$, $a_e = 1.3$, $\theta_e=4$, $a_i=2$, $\theta_i = 3.7$, $r_e=1$, $r_i=1$.$

Effigy 8: Phase Plane Analysis showing limit cycle trajectory in response to constant simulation $P=1.25$. Dashed lines are nullclines. Parameters: $c_1=16$, $c_2 = 12$, $c_3=15$, $c_4=3$, $a_e = 1.iii$, $\theta_e=4$, $a_i=2$, $\theta_i = 3.7$, $r_e=1$, $r_i=1$.

The phase aeroplane analysis illustrates a bounded steady-state solution that is classified as unstable; this is a typical feature of a limit cycle. The solution's oscillating behavior, shown in Figure 10, follows typical limit cycle beliefs:

Trajectories near the equilibrium point are pushed further abroad from the equilibrium.
Trajectories far from the equilibrium point motion closer toward the equilibrium.

We established the resting state $E=0, I=0$ to be stable in the absence of an outside force. Therefore, the neural population merely exhibits limit wheel activeness in response to a abiding external input (P or Q). All in all, the premise of Theorem three is that there is only one steady-state, and it must exist most the inflection indicate for a limit wheel to exist. Therefore, if nosotros study the limit behavior as a function of $P$, where $Q=0$, then it follows that:

There is a threshold value of P, and below this threshold, the limit cycle activity cannot occur.
There is a college value of P, and higher up this bound, the arrangement's limit cycle action will end.
Within the range defined to a higher place, both the limit cycle frequency and the boilerplate value of $Due east(t)$ increases monotonically with respect to $P$.

Figure 9: $I(t)$ and $E(t)$ for limit cycle shown above. The limit cycle depends on the value of $P$, i.e. $Q$ being set equal to zero.

Figure 9: $I(t)$ and $E(t)$ for limit cycle shown above. The limit cycle depends on the value of $P$, i.due east.$Q$ existence set equal to zip.

Summary

The above demonstrates the components necessary to perform qualitative analysis on a one-dimensional autonomous ODE:

plot the flow field
plot several trajectories ane the flow field
identify the equilibrium points
cull a method to make up one's mind stability of equilibrium points

Code Appendix

            library(knitr) library(phaseR) library(deSolve) library(graphics) library(captioner) library(latex2exp) knitr::opts_chunk$prepare(echo = Imitation, out.width = 400, fig.marshal = "middle", message = FALSE) fig_nums <- captioner(prefix = "Figure") body_cap <- fig_nums(proper name = "phase", caption = "Phase Portrait for $\\frac{dy}{dt}=iv-y^2$. The trajectories plotted shows that solutions converge towards $y=2$, just abroad from $y=-ii$.") example1_phasePortrait  <- phasePortrait(   example1, ylim = c(-5, five)) body_cap2 <- fig_nums(proper noun = "phase2", explanation = "The flow field and various trajectories, adding horizontal lines at equilibrium points.") example2_flowField  <- flowField(example2,                                  xlim = c(0, 2),                                  ylim = c(-ane, three),                                  organisation ="one.dim",                                  add = FALSE)  grid()  example2_nullclines <- nullclines(example2,                                   xlim = c(0, 2),                                   ylim = c(-1, iii),                                   organization = "ane.dim",                                   col=c("#ff5ccd"),                                   ltw=ii)  example2_trajectory <- trajectory(example2,                                   y0 = c(-0.5, 0.v, 1.five, 2.five),                                   tlim = c(0, 4),                                   system = "one.dim")  body_cap3 <- fig_nums(proper noun = "phase3", caption = "The period field and various trajectories, adding horizontal lines at equilibrium points.") example2_phasePortrait <- phasePortrait(   example2, ylim = c(-0.5, 2.5)) example2_stability_1 <- stability(   example2, ystar = 0, arrangement = "ane.dim")  example2_stability_2 <- stability(   example2, ystar = ane, system = "i.dim")  example2_stability_3 <- stability(   example2, ystar  = two, organisation = "i.dim") body_cap4 <- fig_nums(name = "phase4", caption = "Plot of the velocity field with several trajectories for $\\frac{dx}{dt} = ten-xy, \ \\frac{dy}{dt} = xy-y$.") lotkaVolterra_flowField    <- flowField(lotkaVolterra,                                         xlim       = c(0, 5),                                          ylim       = c(0, 5),                                         parameters = c(one, ane, one, 1),                                         add        = F) lotkaVolterra_trajectories <- trajectory(lotkaVolterra,                                          y0     = rbind(c(2, 2),                                                         c(0.5, 0.5),                                                         c(0.v, 1.5),                                                         c(ane.five, 0.five),                                                         c(three, 3)),                                          parameters = c(i, one, 1, 1),                                          col    = rep("black", 5),                                          tlim   = c(0, 100)) body_cap5 <- fig_nums(proper name = "phase5", explanation = "Sketch of the nullclines for the system of equations $\\frac{dx}{dt} = ten-xy, \ \\frac{dy}{dt} = xy-y$.") lotkaVolterra_flowField    <- flowField(lotkaVolterra,                                         xlim       = c(-1.v,1.5),                                          ylim       = c(-1.5,1.5),                                         parameters = c(1, 1, 1, 1),                                         add        = F) lotkaVolterra_trajectories <- nullclines(lotkaVolterra,                                          xlim      = c(-1.5,one.5),                                          ylim      = c(-i.v,ane.five),                                          col = c("aquamarine2", "blueviolet"),                                          parameters = c(1, i, 1, one),                                          points = 251) body_cap6 <- fig_nums(proper name = "phase6", explanation = "Plot $x$ and $y$ trajectories against $t$. In this example, nosotros can view the oscillating nature of $x$ and $y$.") lotkaVolterra_numericalSolution <- numericalSolution(lotkaVolterra,                                                      y0 = c(3, four),                                                      tlim = c(0, l),                                                      parameters = c(one, 1, 1, 1)) se <- office(ten){   ae=i.3   theta_e=four   1/(1+exp(-ae*(x-theta_e))) }  si <- part(10){   ai=ii   theta_i= 3.7   1/(ane+exp(-ai*(10-theta_i))) }  WilsonCowan2 <- function(t, y, parameters) {   # couplings   c1 = 16   c2 = 12   c3 = 15   c4 = 3      # Refractory periods   rE = 1   rI = i      # external inputs   P = one.25   Q = 0      ki=0.825   ke=0.88   I <- y[1]   Eastward <- y[ii]   dy <- c(      -I + (ki - rI * I) * si(c3 * E - c4 * I + Q),      -E + (ke - rE * E) * se(c1 * E - c2 * I + P))   listing(dy) } body_cap8 <- fig_nums(proper noun = "phase8", explanation = "Stage Plane Analysis. Determine the steady-land solution by the nullclines' intersection. Parameters: $c_1=16$, $c_2 = 12$, $c_3=15$, $c_4=three$, $a_e = 1.three$, $\\theta_e=four$, $a_i=two$, $\\theta_i = 3.7$, $r_e=1$, $r_i=1$.") example4_flowField  <- flowField(WilsonCowan2,                                  xlim = c(-0.01, .425),                                  ylim = c(0, .five),                                  add together  = Fake,                                  ylab = TeX("$Due east$"),                                  xlab= TeX("$I$"),                                  frac=i) grid() example4_nullclines <- nullclines(WilsonCowan2,                                   xlim = c(-0.01, .425),                                   ylim = c(0, .5),                                   lty = two, lwd = 2,                                   col=c("lightseagreen","aquamarine4"))  body_cap9 <- fig_nums(proper name = "phase9", caption = "Phase Plane Assay showing limit cycle trajectory in response to constant simulation $P=i.25$. Dashed lines are nullclines. Parameters: $c_1=16$, $c_2 = 12$, $c_3=15$, $c_4=3$, $a_e = 1.3$, $\\theta_e=4$, $a_i=ii$, $\\theta_i = 3.7$, $r_e=1$, $r_i=1$.") example4_flowField  <- flowField(WilsonCowan2,                                  xlim = c(-0.01, .425),                                  ylim = c(0, .v),                                  add  = Simulated,                                  ylab = TeX("$E$"),                                  xlab= TeX("$I$"),                                  frac=one) grid() example4_nullclines <- nullclines(WilsonCowan2,                                   xlim = c(-0.01, .425),                                   ylim = c(0, .5),                                   lty = ii, lwd = 2,                                   col=c("lightseagreen","aquamarine4"))  y0 <- matrix(c(0.1,0.3,                0,0,                0.1,0.2), 3, 2, byrow = TRUE) example4_trajectory <- trajectory(WilsonCowan2,                                   y0   = y0,                                   pch=16,                                   tlim = c(0, 100),                                   col="blackness",                                   add=T, ylab=TeX("$E, I$"), xlab=TeX("$t$")) grid() body_cap10 <- fig_nums(proper name = "phase10", caption = "$I(t)$ and $Due east(t)$ for limit cycle shown above. The limit cycle depends on the value of $P$, i.e. $Q$ being set equal to zero.")  example4_solution <- numericalSolution(WilsonCowan2,                         y0 = c(0.0, 0.1),                         type = "two",                         col=c("cornflowerblue", "aquamarine4"),                         add.legend = T,                         xlab = TeX("$t$"),                         ylab = c(TeX("$I$"), TeX("$E$")),                         add.grid = F,                         tlim = c(0,20),                         lwd=two,                         ylim=c(0,0.3),                         xlim=c(0,xviii))

References

[1] Grayling, M. J. (2014). phaseR: An R Parcel for Phase Plane Analysis of Democratic ODE Systems. The R Journal 6 43–51.

davisstagaind.blogspot.com

Source: https://hluebbering.github.io/phase-planes/