phaseR: v2.1

Michael J Grayling (michael.grayling@newcastle.ac.uk)

2022-09-02

phaseR is an R package for the qualitative analysis of one- and two-dimensional autonomous ODE systems, using phase plane methods. Programs are available to identify and classify equilibrium points, plot the direction field, and plot trajectories for multiple initial conditions. In the one-dimensional case, a program is also available to plot the phase portrait. Whilst in the two-dimensional case, additionally programs are available to plot nullclines and the stable/unstable manifolds of saddle points. Many example systems are provided for the user.

Introduction

Contrary to what may seem to be the case when you first encounter ordinary differential equations (ODEs), the majority of ODE systems cannot be solved analytically. In such a scenario, there is usually no option but to resort to numerical solution. However, for certain classes of ODE systems it is possible to undertake a qualitative examination using phase plane methods, as introduced by Henry Poincare, amongst others, in the 19th century. These methods allow the analyser to circumvent the need for explicit solutions, via a highly graphical approach. And in fact, this qualitative analysis is often useful even when the system can be solved analytically. Specifically, it is usually possible to plot trajectories for various initial conditions, before obtaining information about stability and other motion patterns of the system.

This package, phaseR, allows the user to perform such analyses for one- and two-dimensional autonomous ODE systems. Programs are available to determine and classify equilibrium points, plot the flow or vector field, and plot trajectories for multiple initial conditions. In the one-dimensional case, a program is also available to plot the phase portrait. Whilst in the two-dimensional case, additionally programs are available to plot nullclines and the stable/unstable manifolds of saddle points. This accompanying guide has been written not only to provide further information on how to use phaseR, but also as a teaching utility for phase plane methods. In this way, phaseR can hopefully serve as both a package for both independent learning, and for group based teaching; assisting lecturers in explaining the supported techniques.

Thus, since it is an important skill to be able to perform phase plane analysis by hand, and also as a background to the package, this guide will proceed by introducing mathematically the systems that it can examine as well as of the techniques for analysis that it supports. The complexity is most similar to a first year undergraduate mathematics course, however it should hopefully also be useful to those from natural science and engineering backgrounds as well. Following the mathematical descriptions, an explanation of the usage of the programs in phaseR will be given. For this, good knowledge of R is useful, however the programs are not difficult to use. Worked examples will then be provided for both one- and two-dimensional systems. Further example systems available in the package will next be described, before finally, exercises are provided for the user to undertake should they wish. Throughout, to keep things simple, we will stick to using the letters \(x\), \(y\) and \(t\) only as variables, as these are the variable names used by the programs. In practice however, it is obviously not difficult to deal with cases where alternative notation is used.

Acknowledgment goes to Professors Daniel Kaplan and Daniel Flath at Macalester College who established the initial code for phase plane analysis of two dimensional systems (http://www.macalester.edu/~kaplan/math135/pplane.pdf) upon which this package was based. Additionally, I would like to thank Stephen Ellner and John Guckenheimer who graciously provided code which allowed the creation of the findEquilibrium(), drawManifolds(), and phasePlaneAnalyser() functions, as well as Gerhard Burger who helped to substantially advance the package’s inner workings. I welcome any corrections or comments on both the programs and these notes.

First order one-dimensional dynamical systems

Autonomous one-dimensional ordinary differential equations

A first order dynamical system of one variable, \(y(t)\) say, can be written in the following form: \[ \frac{dy}{dt}=f(y,t). \] In many cases (usually the ones found in introductory calculus texts) this ODE can be solved analytically; with several techniques, such as integrating factors and separation of variables, at hand to help. However, more often than not, when a differential equation is written down to describe a real life system, it cannot be solved analytically. This is particularly true of non-linear ODEs, for which numerical solution would frequently have to be utilised. As a result, many computer packages are today available to support such calculations.

However, an alternative approach to numerical integration is sometimes possible. This approach is usually termed the phase plane method, or phase plane analysis. This methodology is concerned with determining qualitative features of the solutions to ODEs, without the need for explicit solution. Whilst such analysis may be more germane to systems we cannot solve analytically, the methods are just as valid to systems we can.

Although this qualitative analysis is indeed possible for ODEs of the type above, in this package we restrict ourselves slightly to the case of ‘autonomous ODEs’, for reasons that should hopefully become clear later. This class of first order ODEs can be written in the following form: \[ \frac{dy}{dt}=f(y), \] i.e., it is the case where this is no dependence upon the independent variable \(t\) in the functional form of \(f\). Moreover, formally, we also assume that \(f\) is a continuous, differentiable function. Whilst this may seem a strong restriction, many real life models can be written in this form.

Now, within this framework of qualitative analysis there are several important concepts that we will proceed to discuss. Namely, the flow field, equilibrium points, and the phase portrait.

The flow field

We begin with a discussion of the flow, or direction, field. Consider again the autonomous ODE above, and imagine making a sketch in the \(t\)-\(y\) plane by drawing at every point \((t,y)\) a small line of slope \(f(y)\). This resulting picture of line segments is the flow field. The use of such a picture lies in that solution curves to the ODE must be tangent to the directions of the line segments. Thus we can construct approximate graphical solutions to this ODE system by beginning at any point \((0,y_0)\) (i.e., \(y(0)=y_0\)) and sketching a curve that proceeds through the plane in the direction of the flow field. In this way, if we start at many different initial points, we can generate a family of solution curves that qualitatively describe the behaviour specified by the ODE.

It is important to note, however, that whilst the flow field method is incredibly useful for plotting trajectories by hand, it is an approximate method: since we can only plot a finite number of line segments some approximation will always be introduced. Usually however, solutions accurate enough to gain a reasonable understanding of the ODE can be achieved, and in general, the more line segments we plot the more accurate our sketches will be. By hand this can be time consuming, utilising a computer however, it is not so difficult.

Note also that some texts on phase plane methods would here discuss the concept of isoclines, defined as lines across which \(dy/dt\) is constant, i.e., where \(f(y)=\alpha\) for different values of \(\alpha\). These lines are used in the same manner as the small line segments of the flow field, since we know the angle at which solution curves should cut them. They however are more useful in the setting of non-autonomous ODEs, and thus we will make little further mention of them.

Additionally, some texts advise to plot the line segments at lengths reflecting the rate of change of \(y\). However, by hand this will almost always be a very laborious task, whilst even with a computer if \(f(y)\) takes a large range of values the resulting plot can become somewhat uninformative with obscuring arrows of great length and other arrows of length too short to be useful. Thus, it is usually best to plot all line segments at some small arbitrary length.

As is often the case in mathematics, concepts can be more easily understood through an example. As such, consider now the ODE: \[ \frac{dy}{dt}=4-y^2, \] provided in the package as example1(). More information will be provided later on how to utilise the programs in phaseR, as well as how to specify your own systems. For now though, simply note the flow field produced below, and the multiple trajectories that follow it:

Equilibrium points and stability

We now turn our attention to the so-called equilibrium points of our autonomous ODE. These points are defined by the locations where: \[ f(y)=0. \] It is easy to understand why they are termed equilibrium points. Beginning at a point \(y_*\) where \(f(y_*)=0\), the system if unperturbed will remain at \(y_*\) throughout its evolution. The great importance of these points lies in determining the long-term behaviour of the ODE.

Considering our example ODE \(dy/dt=4-y^2\) again, it is a simple matter to find its equilibrium points: \[ f(y_*)=0 \Rightarrow 4-y_*^2=0 \Rightarrow (2-y_*)(2+y_*)=0 \Rightarrow y_*=-2,2.\] For the equilibrium points however, just as much as we are interested in their location, we are interested in whether they are stable or unstable. Here, informally, being stable means that if the system is placed a small distance away from the equilibrium point, it will remain close to this equilibrium point. Whilst being unstable means a small perturbation away from the equilibrium point causes the solution to diverge large distances away. More precisely, the definition of stability can be stated as: if for every \(\epsilon>0\), there exists \(\delta>0\) such that whenever \(|y(0)-y_*|<\delta\) then \(|y(t)-y_*|<\epsilon\) for all \(t\).

Classically, to determine the stability of any located equilibrium points, we have two options. The first method is the phase portrait. Indeed, our earlier decision to restrict our attention to autonomous systems was motivated by the condition required for phase portrait analysis: when we remove time dependence from our systems evolution, it allows us to collapse our qualitative analysis from the \(t\)-\(y\) plane to simply considering how \(f(y)\) varies with \(y\).

So, in phase portrait analysis, we first plot \(f(y)\) against \(y\). From the definition of our autonomous ODE it should be easy to see that whenever \(f(y)>0\), \(y\) will increase. Whilst whenever \(f(y)<0\), \(y\) will decrease. Moreover, the locations where \(f(y)\) cross the \(y\)-axis are exactly the equilibrium points (thus this plot can be useful for locating equilibrium points). Therefore, we can represent the evolution of \(y\) in this plot by simply placing arrows along the \(y\)-axis indicating whether \(y\) would be increasing or decreasing. Then, the cases where arrows either side of an equilibrium point towards each other denote stability, whilst when they point away they denote instability.

Again, as an example we consider the system \(dy/dt=4-y^2\). Plotting \(f(y)=4-y^2\) against \(y\) and adding arrows as described we acquire the following graph:

Thus, we can see that the equilibrium point \(y_*=2\) is stable, whilst \(y_*=-2\) is unstable. Indeed, looking back at the trajectories we plotted earlier, we can observe that solutions do converge towards \(y=2\), but away from \(y=-2\).

Moreover, we now note an important consequence of requiring \(f\) to be continuous and differentiable: that the solution curves cannot touch each other (except to converge at equilibrium points). This is because these conditions on \(f\) guarantee solutions to autonmous ODEs are unique. We can observe in our earlier plot of several trajectories of the system \(dy/dt=4-y^2\) that this is indeed the case.

Our second option to perform such stability analyses comes from utilising the Taylor Series expansion of \(f\). We begin by supposing we are a small distance \(\delta(0)\) away from our fixed point \(y_*\), i.e., \(y(0)=y_*+\delta(0)\), and in general that \(y(t)=y_*+\delta(t)\). Then we can write the Taylor Series of \(f\) as: \[ f(y_*+\delta)=f(y_*)+\delta\frac{\partial f}{\partial y}(y_*)+o(\delta), \] assuming higher order terms can be neglected. Recalling \(f(y_*)=0\), our autonomous ODE becomes: \[ \frac{d}{dt}(y_*+\delta)=\delta\frac{\partial f}{\partial y}(y_*) \Rightarrow \frac{d\delta}{dt}=\delta\frac{\partial f}{\partial y}(y_*)=k\delta. \] This ODE for \(\delta\) can be solved easily to give \(\delta(t)=\delta(0)e^{kt}\). Then stability can be found based upon whether \(\delta(t)\) grows or decays as \(t\) increases, i.e., we have: \[ k = \frac{\partial f}{\partial y}(y_*)\begin{cases} > 0 &: \ \text{Stable},\\ < 0 &: \ \text{Unstable}. \end{cases} \] Here, \(k\) is sometimes referred to as the discriminant, whilst this approach is often referred to as perturbation analysis.

Returning to our example ODE again, we can perform such analysis easily: \[ \frac{\partial f}{\partial y}(y_*) = -2y_*= \begin{cases} -4 &: \ y_*=2, \\ +4 &: \ y_* = -2. \end{cases}\] Thus we draw the same conclusion as before; \(y_*=2\) is stable, and \(y_*=-2\) is unstable. We will see later how one of the programs in phaseR can perform this stability analysis for us.

It should now be clear that we can clearly state if \(y(0)>2\) or \(0<y(0)<2\) then the solution will eventually approach \(y=2\). However, if \(y(0)<0\), \(y\rightarrow-\infty\) as \(t\rightarrow\infty\). Such general statements can often be made as a result of the phase plane analyses.

It is worthwhile noting here that if we find: \[ \frac{\partial f}{\partial y}(y_*)=0, \] then to this order of the Taylor Series no conclusion can be drawn about stability.

So, we have now observed all of the key components required to perform a qualitative analysis upon a one-dimensional autonomous ODE. We begin by plotting the flow field, and from this several trajectories. We then identify the equilibrium points and choose a method to determine their stability. All such techniques are available in this package, and we will later discuss how to implement them. First, however, we will discuss how these methods can be generalised to coupled ODEs.

First order two-dimensional dynamical systems

Autonomous two-dimensional ordinary differential equations

As may well be expected, things get substantially more complex in the world of coupled ODEs; very rarely can such systems be solved analytically. Unfortunately, the analysis of many real life systems does involve interacting variables, and so these systems are not uncommon. Here, the first restriction we make is to the case of two-dimensional (or two variable) systems; a necessity for the following techniques to be implementable (this is often considered a disadvantage of phase plane methods; that they cannot be generalised to more than two dimensions. Fortunately however, many systems can be approximated to two dimensions). These systems can be written in their most general form as: \[ \frac{dx}{dt}=f(x,y,t),\ \frac{dy}{dt}=g(x,y,t), \] for \(x=x(t)\) and \(y=y(t)\). In this most general case numerical solution would almost certainly be the only way forward. However, if we again make the restriction to autonomous systems, the phase plane methods from one dimension can be generalised to avoid the need for numerical integration. Following the same route as in the one-dimensional case, an autonomous system can be written for two coupled ODEs as: \[ \frac{dx}{dt}=f(x,y),\ \frac{dy}{dt}=g(x,y). \] As before, the definition of the flow field (more commonly, and from here on out, referred to as the velocity field) and equilibrium points, as well as their stability, will be important. Here however, we also meet the concept of nullclines and stable/unstable manifolds. Again, we formally require that \(f\) and \(g\) are continuous, (and now) partially differentiable functions.

Before we proceed to discuss the generalisation of our earlier techniques to such two-dimensional systems, it is useful to note that certain second order ODEs can indeed be re-cast by variable substitution into a system of the type above. Indeed, consider the second order ODE given by: \[ a(y)\frac{d^2y}{dt} + b(y)\frac{dy}{dt} + c(y) = 0. \] We make the substitution \(x=dy/dt\) and re-write our system as: \[ \frac{dy}{dt} = x,\ \frac{dx}{dt} = \frac{1}{a(y)}\{-b(y)x-c(y)\}. \] In this way, it is actually possible to analyse the behaviour of certain second order ODEs using the methods for coupled first order ODEs.

The velocity field

We observed earlier that the restriction to autonomous ODEs in the one-dimensional case allowed us to restrict attention to the phase portrait; the plot of \(f(y)\) against \(y\). In the two-dimensional case, this restriction allows us to limit our attention to the plane produced by the two dependent variables. Using our notation from the two-dimensional autonomous system prescribed above, this is the \(x\)-\(y\) plane, and is often referred to in this context as the phase plane. Representation in this manner proves to be the most convenient way to visualise the system.

In this plane, we can produce a plot analogous to the flow field discussed earlier for one-dimensional autonomous systems: at many points \((x,y)\) we plot a small line segment (a vector) in the direction given by the rates of change of \(x\) and \(y\), which are provided by \(f(x,y)\) and \(g(x,y)\). This plot is usually referred to as the velocity field, or sometimes the direction field, and perhaps confusingly, the phase portrait. We can then again for any point trace out the trajectory of a solution by using the fact that it must pass through our line segments in a tangential manner. Repeating this procedure for several starting points, we can again build up a family of solutions and a good picture of the behaviour of solutions to our two-dimensional autonomous ODE system. As before, however, it is important to understand that using this method is only an approximation to performing numerical integration, and things can here become very ambiguous around certain points (the equilibria).

To illustrate the concept of the velocity field, we again turn to an example. This time consider the system given by (this is an example of a Lotka-Volterra model, which will be discussed more later): \[ \frac{dx}{dt} = x-xy, \ \frac{dy}{dt} = xy-y. \] Using phaseR we can produce the following plot of the velocity field along with several trajectories:

Analogous to the one-dimensional analysis performed earlier, we observe how our restriction to continuous partially differentiable \(f\) and \(g\) ensures that trajectories cannot cross (though they can again converge at equilibria).

What is more, as before, some texts advise to plot the vectors at lengths reflecting the magnitudes of the rates of change of \(x\) and \(y\). However, some small arbitrary length usually still remains the best option.

Finally, as was the case in the one-dimensional analysis, some texts here again refer to the method of isoclines for tracing out trajectories. Isoclines are now defined as curves of constant gradient in the \(x\)-\(y\) plane, i.e.: \[ \frac{dy}{dx} = \frac{g(x,y)}{f(x,y)} = \alpha, \] for different values of \(\alpha\). Once more, trajectories would be produced by using the fact that we know the angle they should cut each isocline. We will make no further reference to isoclines in these notes; hopefully for reasons discussed below it should become clear why certain tricks make the need for to plot isoclines very unusual.

Nullclines

A new important concept in the case of two-dimensional systems is that of nullclines. Here, \(x\)-nullclines are defined by the locations where \(f(x,y)=0\), whilst the \(y\)-nullclines are defined by the locations where \(g(x,y)=0\). Thus, the \(x\)- and \(y\)-nullclines define the locations where \(x\) and \(y\), respectively, do not change with \(t\). As a consequence, when plotting a vector field by hand it is usually wise to plot the nullclines first, as the line segments (or vectors) along them move parallel to the \(x\)- and \(y\)-axes.

Returning to our example given above, we can find its nullclines as follows: \[ \begin{align} x &∶ x-xy = 0 \Rightarrow x(1-y) = 0 \Rightarrow x=0 \text{ or } y=1,\\ y &∶ xy-y = 0 \Rightarrow y(x-1) = 0 \Rightarrow x=1 \text{ or } y=0. \end{align} \] We can then plot these nullclines along with the velocity field:

Equilibrium points and stability

Equilibrium points maintain their importance in two dimensions. Here, generalisation defines them to be the locations \((x_*,y_*)\) where: \[ f(x_*,y_*)=g(x_*,y_*)=0. \] Thus, another utility of nullclines immediately becomes apparent; the locations where \(x\)- and \(y\)-nullclines cross are the equilibria. However, it is important to note that locations where \(x\)-nullclines or \(y\)-nullclines cross each other, are not equilibria. For this reason it is usually useful to plot \(x\)- and \(y\)-nullclines in different colours.

Revisiting our example system from earlier, it is easy to find either analytically, or from the nullcline plot, that two equilibria are present; the points \((0,0)\) and \((1,1)\).

We now note a useful fact about equilibria and nullclines, which can be observed in the plot above. On opposite sides of an equilibria, along a nullcline, the orientation of the velocity arrows is reversed. This is a property shared by the majority of systems (with the exception being certain singular cases where the Jacobian that we meet later is zero). Because trajectories must be continuous, the direction vectors must vary continuously from one point to another everywhere else. So in most cases when seeking to plot the velocity field, it suffices to determine direction vectors at a few select locations and deduce the rest by preserving continuity and switching orientation when an equilibrium point is crossed. It is this trick that makes plotting numerous isoclines often unnecessary.

Again, we must now turn our attention to the stability of the equilibrium points. In two dimensions the definition of stability remains the same, but as well as determining whether a point is stable or unstable, we can additionally classify the nature in which trajectories move away or towards it. Ultimately, as in the one-dimensional case, we aim to identify the long term behaviour of solutions in different regions of the plane.

In this case, we must make use of a mathematical argument; to gain a full understanding use of a graph is not enough (though it can be useful in certain singular cases). Here, many texts distinguish between the case of linear and non-linear systems, and so we will also first make such a distinction, though ultimately we will treat these two types of system equally.

The linear version of a two-dimensional autonomous ODE system is given by: \[ \frac{dx}{dt} = ax+by,\ \frac{dy}{dt} = cx+dy, \] or in matrix form: \[ \frac{d}{dt}\begin{pmatrix} x \\ y\end{pmatrix} = \frac{d}{dt}\boldsymbol{x} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} x \\ y\end{pmatrix} = A \boldsymbol{x}. \] Now, provided \(\text{det}(A)\neq0\), a unique solution exists to this system (we will return to discuss the case \(\text{det}(A)=0\) later), and it can be written in the form: \[ \boldsymbol{x}=C_1e^{\lambda_1t}\boldsymbol{e}_1 + C_2e^{\lambda_2t}\boldsymbol{e}_2, \] where \(\lambda_1\) and \(\lambda_2\) are the eigenvalues of \(A\), \(\boldsymbol{e}_1\) and \(\boldsymbol{e}_2\) are their corresponding eigenvectors, and \(C_1\) and \(C_2\) are arbitrary constants. From here we can determine stability based on the values of the eigenvalues. However, the procedure from here on out is the same as that for non-linear systems and so we will move to the analysis required for the more complex case of non-linearity.

So, from the above it should be obvious that provided \(\text{det}(A)\neq0\), linear systems have only one equilibrium point; \((0,0)\). Non-linear systems, however, are much more complicated; they can have multiple equilibria and even display limit cycle behaviour (as defined later). However, close to an equilibrium point, behaviour can be usually understood by linearising the model about the equilibria.

To do this we proceed in a similar fashion to the Taylor Series method for one-dimensional autonomous ODEs. We suppose we have an equilibrium point given by \((x_*,y_*)\) and that our system lies initially slightly away from this point at \((x_*+\delta(0),y_*+\epsilon(0))\), and in general at \((x_*+\delta(t),y_*+\epsilon(t))\). Then, using the Taylor expansion for \(f\), our differential equation for \(\boldsymbol{x}\) becomes: \[ \begin{align} \frac{d\delta}{dt} &= f(x_*+\delta,y_*+\epsilon),\\ &= f(x_*,y_*) + \delta \frac{\partial f}{\partial x}(x_*,y_*) + \epsilon\frac{\partial f}{\partial y}(x_*,y_*) + o(\delta) + o(\epsilon),\\ &= \delta \frac{\partial f}{\partial x}(x_*,y_*) + \epsilon\frac{\partial f}{\partial y}(x_*,y_*) + o(\delta) + o(\epsilon). \end{align} \] Similarly, our differential equation for \(y\) becomes: \[ \frac{d\epsilon}{dt} = \delta \frac{\partial g}{\partial x}(x_*,y_*) + \epsilon\frac{\partial g}{\partial y}(x_*,y_*) + o(\delta) + o(\epsilon). \] Here we have again assumed terms of second order and higher are negligible.

If we write this system in matrix form we acquire: \[ \frac{d}{dt}\boldsymbol{\delta} = \left. \begin{pmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{pmatrix}\right|_{(x_*,y_*)}\boldsymbol{\delta} = \left. \begin{pmatrix} f_x & f_y \\ g_x & g_y \end{pmatrix}\right|_{(x_*,y_*)}\boldsymbol{\delta} = J\boldsymbol{\delta},\] where \(J\) is called the Jacobian of the system, and: \[ \boldsymbol{\delta} = \begin{pmatrix} \delta \\ \epsilon \end{pmatrix}. \] If we let the eigenvalues of \(J\) be denoted by \(\lambda_1\) and \(\lambda_2\), with corresponding eigenvectors \(\boldsymbol{e}_1\) and \(\boldsymbol{e}_2\), then the general solution to the above is: \[ \boldsymbol{\delta} = C_1e^{\lambda_1t}\boldsymbol{e}_1+C_2e^{\lambda_2t}\boldsymbol{e}_2, \] where \(C_1\) and \(C_2\) are arbitrary constants.

Considering the linear system from earlier we find that \(J=A\). Thus, stability of \((0,0)\) in the linear case can be determined by the same classification rules as below for the non-linear case. Specifically, we have:

Fortunately, it is not actually necessary to find the exact values of the eigenvalues (though computationally this is not a difficult task, by hand it can be time consuming). We only require the signs of the eigenvalues, or of their real parts, to perform the classification. To this end, consider the characteristic equation of \(J\): \[ (f_x-λ)(g_y-λ)-f_y g_x=0. \] Observing that \(\text{tr}(J)=T=f_x+g_y\) and \(\text{det}(J)=\Delta=f_x g_y-f_y g_x\), we can write the characteristic equation of \(J\) as: \[ \begin{align} \lambda^2 - T\lambda+\Delta &= 0,\\ \Longrightarrow \lambda &= \frac{T \pm \sqrt{T^2-4\Delta}}{2}. \end{align}\] From this we can draw up the following table that allows us to classify the equilibria using the signs of \(T\), \(\Delta\), and \(T^2-4\Delta\):

\(\Delta\) \(T^2-4\Delta\) Eigenvalues of \(J\) \(T\) Classification
\(<0\) \(>0\) Real, opposite signs N/A Saddle
\(>0\) \(>0\) Real, same signs \(<0\) Stable node
\(>0\) \(>0\) Real, same signs \(>0\) Unstable node
\(>0\) \(<0\) Complex conjugate pair \(<0\) Stable focus
\(>0\) \(<0\) Complex conjugate pair \(=0\) Centre
\(>0\) \(<0\) Complex conjugate pair \(>0\) Unstable focus
\(=0\) N/A N/A Indeterminate
N/A \(=0\) Real, equal \(<0\) Stable node
N/A \(=0\) Real, equal \(>0\) Unstable node

Note: Focus’ are often referred to as spirals.

From here, we will refer to \(T^2-4\Delta\) as the discriminant.

To be more precise, for the case of \(\Delta=0\); we would have to consider second-order terms in the Taylor Series approximation made earlier in order to determine stability. Alternatively, in this case, use of the velocity field and traced trajectories can allow us to identify if the point is stable or not.

Returning to our example system from earlier, taking partial derivatives we can compute the Jacobian at any equilibrium point \((x_*,y_*)\): \[ J = \begin{pmatrix} 1 - y_* & -x_* \\ y_* & x_* - 1\end{pmatrix}. \] Thus, at \((0,0)\), we have: \[ J = \left.\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\right|_{(0,0)}. \] So \(\text{tr}(J)=T=0\) and \(\text{det}(J)=\Delta=-1\); which from our table above makes \((0,0)\) a saddle point. For \((1,1)\), however, we have: \[ J = \left.\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}\right|_{(1,1)}. \] Therefore, \(\text{tr}(J)=T=0\) and \(\text{det}(J)=\Delta=1\); which from our table above makes \((1,1)\) a centre. Indeed, if we look back at our earlier plot, we can observe trajectories diverging away from \((0,0)\), but traversing around \((1,1)\). Again, we will see later how this analysis can be performed for us in phaseR.

As a last point of interest, note that it is sometimes interesting to plot \(x\) and \(y\) trajectories against \(t\). For the case of \((x_0,y_0)=(3,4)\) in our example system this results in the following plot where we can witness the oscillating nature of \(x\) and \(y\):

The utility of such plots becomes more apparent in cases where trajectories can be seen to converge upon an equilibrium point; indicating its stability and often whether it is a node or focus.

So, we have now discussed all of the techniques required to a perform phase plane analysis of a two-dimensional autonomous ODE system. We begin by locating and plotting nullclines, using these to create the velocity field. From this we can plot numerous trajectories. We then identify any equilibria and classify them according to the table above. Before we move to discuss how to use phaseR however, we will give brief attention to two additional considerations of interest in the analysis of two-dimensional autonomous ODE systems.

Stable and unstable manifolds

The stable and unstable manifolds of an equilibrium point give a formal mathematical definition to the general notions embodied in the idea of a point being an attractor or repellor. Precisely, the stable and unstable manifolds are defined as the set of all points in the plane which tend to the equilibrium point as time goes to positive, and respectively negative, infinity. For example, a system with a single stable node will have the whole plane as its stable manifold and just the node itself as its unstable manifold.

Things are slightly more interesting for saddle points, however. Recall that saddle points have a positive eigenvalue and a negative eigenvalue. The stable (unstable) manifold of a saddle point is the eigenvector corresponding to the positive (negative) eigenvalue. Moreover, the stable and unstable manifolds of saddle points are very important to understanding the global behaviour of the orbits of the whole planar system. They form separatrices; partitioning the plane into invariant regions of differing dynamics. We will see evidence of such separatrices in an example later when we observe how phaseR can be used to plot the stable and unstable manifolds of any saddle point.

Limit cycles

Non-linear systems can also exhibit a type of behaviour known as a limit cycle. In the phase plane, a limit cycle is defined as an isolated closed orbit. Closed here denotes the periodic nature of the motion and isolated denotes the limiting nature of the cycle; with nearby trajectories converging to, or diverging away from, it. Limit cycles have a complex mathematical theory behind them, which we will not go into here. We will however observe an example of limit cycle behaviour later on.

phaseR usage

To perform all of the techniques discussed above, the package contains nine key functions. Below is a description of each ones utility, as well as the user specifiable input variables, and the output you can expect. The description of the inputs and outputs is repetitive on purpose to reflect how many are common across programs. For the most part, they are also equal to the descriptions seen in the function help files in R. In addition, we will continue to use \(x\), \(y\), and \(t\) as our variables.

drawManifolds()

drawManifolds(deriv, y0 = NULL, parameters = NULL, tstep = 0.1, tend = 1000,
              col = c("green", "red"), add.legend = TRUE,
              state.names = c("x", "y"), ...)

This function allows the user to plot the stable and unstable manifolds of a saddle point. The following inputs can be set:

Returned by drawManifolds() is a list containing all of the input variables as well as following components:

findEquilibrium()

findEquilibrium(deriv, y0 = NULL, parameters = NULL, system = "two.dim",
                tol = 1e-16, max.iter = 50, h = 1e-6, plot.it = FALSE,
                summary = TRUE,
                state.names = if (system == "two.dim") c("x", "y") else "y")

This function allows the user to search for an equilibrium point. The following inputs can be set:

Returned by findEquilibrium() is a list object containing all of the input variables as well as following components (the exact make up is dependent on the value of system):

flowField()

flowField(deriv, xlim, ylim, parameters = NULL, system = "two.dim", points = 21,
          col = "gray", arrow.type = "equal", arrow.head = 0.05, frac = 1,
          add = TRUE,
          state.names = if (system == "two.dim") c("x", "y") else "y",
          xlab = if (system == "two.dim") state.names[1] else "t",
          ylab = if (system == "two.dim") state.names[2] else state.names[1],
          ...)

This function allows the user to plot the flow or velocity field for a one- or two-dimensional autonomous ODE system. The following inputs can be set:

Returned by flowField() is a list object containing all of the input variables as well as following components (the exact the exact make up is dependent on the value of system):

nullclines()

nullclines(deriv, xlim, ylim, parameters = NULL, system = "two.dim",
           points = 101, col = c("blue", "cyan"), add = TRUE, add.legend = TRUE, 
           state.names = if(system == "two.dim") c("x", "y") else "y", ...)

This function allows the user to plot nullclines for two-dimensional autonomous ODE systems. Or it can be used to plot horizontal lines at equilibrium points for one dimensional autonomous ODE systems. The following inputs can be set:

Returned by nullclines() is a list object containing all of the input variables as well as following components (the exact the exact make up is dependent upon the value of system):

numericalSolution()

numericalSolution(deriv, y0 = NULL, tlim, tstep = 0.01, parameters = NULL,
                  type = "one", col = c("red", "blue"), add.grid = TRUE, 
                  add.legend = TRUE, state.names = c("x", "y"), xlab = "t",
                  ylab = state.names)

Used for two-dimensional systems, this function numerically solves the autonomous ODE system for a given initial condition. It then plots the dependent variables against the independent variable. The following inputs can be set:

Here, the numerical integration is performed by the function deSolve::ode().

Returned by numericalSolution() is a list object containing all of the input variables as well as following:

phasePlaneAnalysis()

phasePlaneAnalysis(deriv, xlim, ylim, tend = 100, parameters = NULL,
                   system = "two.dim", add = FALSE,
                   state.names = if (system == "two.dim") c("x", "y") else "y")

This function allows the user to perform a basic phase plane analysis and produce a simple plot without the need to use the other functions directly. Specifically, a list is provided and the user inputs a value to the console to decide what is added to the plot. The following inputs can be set:

phasePortrait()

phasePortrait(deriv, ylim, ystep = 0.01, parameters = NULL, points = 10,
              frac = 0.75, arrow.head = 0.075, col = "black", add.grid = TRUE, 
              state.names = "y", xlab = state.names,
              ylab = paste0("d", state.names), ...)

For a one-dimensional autonomous ODE, it plots the phase portrait, i.e., the derivative against the dependent variable. In addition, along the dependent variable axis it plots arrows pointing in the direction of dependent variable change with increasing value of the independent variable. From this stability of equilibrium points (i.e., locations where the horizontal axis is crossed) can be determined:

Returned by phasePortrait() is a list object containing all of the input variables as well as the following:

stability()

stability(deriv, ystar = NULL, parameters = NULL, system = "two.dim", h = 1e-07, 
          summary = TRUE,
          state.names = if (system == "two.dim") c("x", "y") else "y")

Uses stability analysis to classify equilibrium points. Uses the Taylor Series approach (also known as Perturbation Analysis) to classify equilibrium points of a one dimensional autonomous ODE system, or the Jacobian approach to classify equilibrium points of a two dimensional autonomous ODE system. The following inputs can be set:

Returned by stability() is a list object containing all of the input variables as well as following components (the exact the exact make up is dependent upon the value of system):

trajectory()

trajectory(deriv, y0 = NULL, n = NULL, tlim, tstep = 0.01, parameters = NULL,
           system = "two.dim", col = "black", add = TRUE,
           state.names = if (system == "two.dim") c("x", "y") else "y", ...)

This function allows the user to plot trajectories by performing numerical integration of the chosen ODE system, for a user specifiable range of initial conditions. The following inputs can be set:

Here, the numerical integration is performed by the function deSolve::ode().

Returned by trajectory() is a list object containing all of the input variables as well as following components (the exact the exact make up is dependent on the value of system):

Derivative specification

In addition to the above functions, phaseR contains multiple example one- and two-dimensional autonomous ODE systems; these are the focus of later sections. Here, we discuss how the user can create their own system.

In order to be compatible with phaseR, system functions need to be compatible with deSolve::ode(), which is used internally to perform numerical integration. Thus, detailed guidance on how to write such functions is available in the deSolve documentation. A brief overview is provided here.

Functions that detail a system need to be coded to return a list, and must take three inputs: t, y, and parameters. Thus the basic skeleton for a one- or two-dimensional system (with the function named derivative) is as follows:

derivative <- function(t, y, parameters) {
  # Enter derivative computation here
  list(dy)
}

All that needs to be done is to set the value of dy, using t, y, and parameters (note here that t should not generally be used explicitly as we work only with autonomous systems). The approach should change slightly depending on whether you are setting up a one- or two-dimensional system, as dy should be a vector that has length equal to the number of dimensions.

Thus, for a system such as: \[ \frac{dx}{dt}=3y,\qquad \frac{dy}{dt}=2x, \] we could use the following code:

derivative <- function(t, y, parameters) {
  x  <- y[1]
  y  <- y[2]
  dy <- c(3*y, 2*x)
  list(dy)
}

As a more complex example, consider instead changing the system above to: \[ \frac{dx}{dt} = \alpha y,\qquad \frac{dy}{dt}=\beta x, \] with \(\alpha\) and \(\beta\) parameters. The code would then proceed as follows:

derivative <- function(t, y, parameters) {
  alpha <- parameters[1]
  beta  <- parameters[2]
  x     <- y[1]
  y     <- y[2]
  dy    <- c(alpha*y, beta*x)
  list(dy)
}

Things are slightly simpler for one-dimensional systems. We could for example create a derivative function for the system: \[ \frac{dy}{dt}=a(b-3-y)^2, \] using the following code:

derivative <- function(t, y, parameters) {
  a  <- parameters[1]
  b  <- parameters[2]
  dy <- a*(b – 3 - y)^2
  list(dy)
}

Examples

Within phaseR, numerous example systems are available. Here we will analyse some of them in order to indicate how to perform phase plane analyses by hand or with the help of phaseR. The language is again at times deliberatively repetitive; here to indicate how a general procedure can be used when performing analyses. It is not a useful exercise for me to provide inferior hand drawn plots, only ones produced by phaseR are shown. It is in addition useful to note that the small circles on the trajectory plots indicate initial conditions specified by the user.

example2()

We begin with the one-dimensional autonomous ODE: \[ \frac{dy}{dt}=y(1-y)(2-y), \] provided in the package as example2(). We begin by plotting the flow field and several trajectories using the following code, adding horizontal lines at any equilibrium points to indicate their presence as well:

example2_flowField  <- flowField(example2,
                                 xlim    = c(0, 4),
                                 ylim    = c(-1, 3),
                                 system  = "one.dim",
                                 add     = FALSE)
grid()
example2_nullclines <- nullclines(example2,
                                  xlim   = c(0, 4),
                                  ylim   = c(-1, 3),
                                  system = "one.dim")
example2_trajectory <- trajectory(example2,
                                  y0     = c(-0.5, 0.5, 1.5, 2.5),
                                  tlim   = c(0, 4),
                                  system = "one.dim")
#> Note: col has been reset as required

Thus, three equilibrium points have been identified; appearing to be \(y_*=0,1,2\) Indeed if we set the RHS of our ODE to zero we can identify these three points as the equilibrium points analytically: \[ y_* (1-y_* )(2-y_* )=0 \Longrightarrow y_*=0,1,2. \] Plotting the phase portrait, we find that \(y_*=0\) and \(y_*=2\) are unstable, whilst \(y_*=1\) is stable (as is also apparent from the flow field and trajectories above):

example2_phasePortrait <- phasePortrait(example2,
                                        ylim = c(-0.5, 2.5))

Alternatively, using the Taylor Series approach we have: \[ \frac{d}{dy}\left.\left(\frac{dy}{dt}\right)\right|_{y=y_*} = 3y_*^2 - 6y_* + 2 = \begin{cases} 2 &: \ y_*=0,\\ -1 &: \ y_*=1,\\ 2 &: \ y_*=2. \end{cases} \] Thus we draw the same conclusion as from the phase portrait.

Finally, we can confirm our Taylor analysis using stability() and the following code:

example2_stability_1 <- stability(example2,
                                  ystar  = 0,
                                  system = "one.dim")
#> discriminant = 2, classification = Unstable
example2_stability_2 <- stability(example2,
                                  ystar  = 1,
                                  system = "one.dim")
#> discriminant = -1, classification = Stable
example2_stability_3 <- stability(example2,
                                  ystar  = 2,
                                  system = "one.dim")
#> discriminant = 2, classification = Unstable

logistic()

The logistic growth model is frequently used in biology to model the growth of a population under density dependence. It is given by: \[ \frac{dy}{dt} = \beta y\left(1-\frac{y}{K}\right). \] With the following code, we can plot the flow field and several trajectories (for the case \(\beta=1\) and \(K=2\)), as well as adding horizontal lines at any equilibrium points to indicate their presence:

logistic_flowField  <- flowField(logistic,
                                 xlim       = c(0, 5),
                                 ylim       = c(-1, 3),
                                 parameters = c(1, 2),
                                 system     = "one.dim",
                                 add        = FALSE)
grid()
logistic_nullclines <- nullclines(logistic,
                                  xlim       = c(0, 5),
                                  ylim       = c(-1, 3),
                                  parameters = c(1, 2),
                                  system     = "one.dim")
logistic_trajectory <- trajectory(logistic,
                                  y0         = c(-0.5, 0.5, 1.5, 2.5),
                                  tlim       = c(0, 5),
                                  parameters = c(1, 2),
                                  system     = "one.dim")
#> Note: col has been reset as required

Again, two equilibrium points have been identified. We can confirm their location in the general case analytically by setting the RHS of the ODE to zero: \[ \beta y_* \left(1-\frac{y_*}{K}\right)=0 \Longrightarrow y_*=0,K. \] Plotting the phase portrait we can observe that \(y_*=0\) is unstable and \(y_*=K\) stable, for the case \(\beta=1\) and \(K=2\):

logistic_phasePortrait <- phasePortrait(logistic,
                                        ylim       = c(-0.5, 2.5),
                                        parameters = c(1, 2))

Finally, if we use our Taylor Series method we can draw this same conclusion: \[ \frac{d}{dy}\left.\left(\frac{dy}{dt}\right)\right|_{y=y_*} = \beta - \frac{2\beta y_*}{K} = \begin{cases} \beta &: \ y_*=0,\\ -\beta &: \ y_*=K. \end{cases} \] So for \(\beta=1\) and \(K=2\), we have a stable point at \(y=2\). Moreover, from this we can see that the point \(y=K\) will in general be stable provided \(\beta>0\).

The following code verifies our findings for the specific case studied above:

logistic_stability_1 <- stability(logistic,
                                  ystar      = 0,
                                  parameters = c(1, 2),
                                  system     = "one.dim")
#> discriminant = 1, classification = Unstable
logistic_stability_2 <- stability(logistic,
                                  ystar      = 2,
                                  parameters = c(1, 2),
                                  system     = "one.dim")
#> discriminant = -1, classification = Stable

example4()

We now turn our attention to linear two-dimensional autonomous ODE systems. Here we consider the coupled system given by: \[ \frac{dx}{dt}=-x,\qquad \frac{dy}{dt}=4x. \] This is provided in the package as example4(). We can find the \(x\)- and \(y\)-nullclines by setting the derivatives to zero as follows: \[ \begin{align} x &∶ -x=0 \Longrightarrow x=0,\\ y &∶ 4x=0 \Longrightarrow x=0. \end{align}\] Thus the nullclines are the same. This means we have a line of equilibrium points given by \(x=0\). To see why this is the case, and there is no unique solution, we examine the Jacobian of our system: \[ J = \begin{pmatrix} -1 & 0 \\ 4 & 0 \end{pmatrix}. \] Thus the determinant of \(J\) is zero, and we have a singular case of the general linear two-dimensional system; the Taylor approach cannot be used to classify \((0,0)\).

Thus, here, in order to determine whether the points along the line \(x=0\) are stable or not, we plot the nullclines, the velocity field, and add several trajectories:

example4_flowField  <- flowField(example4,
                                 xlim = c(-3, 3),
                                 ylim = c(-5, 5),
                                 add  = FALSE)
grid()
example4_nullclines <- nullclines(example4,
                                  xlim = c(-3, 3),
                                  ylim = c(-5, 5))
y0                  <- matrix(c(1, 0, -1, 0, 2, 2, -2, 2, 3, -3, -3, -4),
                              6, 2, byrow = TRUE)
example4_trajectory <- trajectory(example4,
                                  y0   = y0,
                                  tlim = c(0, 10))
#> Note: col has been reset as required

Thus we observe the trajectories moving towards the line \(x=0\); indicative of stability. This example illustrates that plotting trajectories can be useful when the Taylor approach fails.

example5()

We will now examine a further example of a linear two-dimensional system, given by: \[ \frac{dx}{dt}=2x+y,\qquad \frac{dy}{dt}=2x-y. \] It is provided in the package as example5(). Again we begin by setting the derivatives to zero to identify the nullclines: \[ \begin{align} x &∶ 2x+y=0 \Longrightarrow y=-2x,\\ y &∶ 2x-y=0\Longrightarrow y=2x. \end{align} \] From these two equations it is easy to see that the only equilibrium point is at \((0,0)\). We begin by plotting the nullclines, velocity field, and several trajectories (along with the stable and unstable manifolds of the saddle, we discuss below):

example5_flowField  <- flowField(example5,
                                 xlim = c(-3, 3), 
                                 ylim = c(-3, 3),
                                 add  = FALSE)
grid()
example5_nullclines <- nullclines(example5,
                                  xlim = c(-3, 3), 
                                  ylim = c(-3, 3))
y0                  <- matrix(c(1, 0, -1, 0, 2, 2, -2, 2, 0, 3, 0, -3),
                              6, 2, byrow = TRUE)
example5_trajectory <- trajectory(example5,
                                  y0   = y0,
                                  tlim = c(0, 10))
#> Note: col has been reset as required
example5_manifolds  <- drawManifolds(example5,
                                     y0 = c(0, 0))

From the trajectories it appears that \((0,0)\) is a saddle point. To verify this we compute the Jacobian of the system: \[ J = \begin{pmatrix} 2 & 1 \\ 2 & -1 \end{pmatrix}. \] Thus we have \(T=1\), \(\Delta=-4\) and \(T^2-4\Delta=17\). From our classification rules this confirms that \((0,0)\) is indeed a saddle. Finally, we verify this analysis using stability and the following code:

example5_stability <- stability(example5,
                                ystar = c(0, 0))
#> tr = 1, Delta = -4, discriminant = 17, classification = Saddle

example8()

As a final example of a linear two-dimensional system, we will examine: \[ \frac{dx}{dt}=y,\qquad \frac{dy}{dt}=-x-y, \] available in the package as example8(). Setting the derivatives to zero, we first locate the nullclines: \[ \begin{align} x &∶ y = 0,\\ y &∶ -x-y=0\Longrightarrow y=-x. \end{align} \] From this, again we can identify the one equilibrium is at \((0,0)\). We now plot the nullclines and velocity field, along with several trajectories:

example8_flowField  <- flowField(example8,
                                 xlim = c(-3, 3),
                                 ylim = c(-3, 3),
                                 add  = FALSE)
grid()
example8_nullclines <- nullclines(example8,
                                  xlim = c(-3, 3), 
                                  ylim = c(-3, 3))
y0                  <- matrix(c(1, 0, 0, 0.5, 2, -2, -2, -2),  
                              4, 2, byrow = TRUE)
example8_trajectory <- trajectory(example8,
                                  y0   = y0,
                                  tlim = c(0, 10))
#> Note: col has been reset as required