Phase plane analysis

As mentioned in the section on modeling approaches, in many cases there is no analytical solution to systems with nonlinear interacting dynamics. An alternative to numerical simulation is phase plane analysis. It is an important technique for studying the behavior of nonlinear systems, consisting of the following analytical steps:

Steady states (aka fixed points)

Steady states (or fixed points) are states in the development of a system at which the dynamics do not show any further changes. A way to find steady states is to set the differential equation expressing the change of the system per time step to 0 (i.e. to no change):

Consider the system \(\{𝑓(π‘₯,𝑦),𝑔(π‘₯,𝑦)\}\) with \(x\) and \(y\) being two interacting stocks changing in time according to the following differential equations:

\[\frac{dx}{dt}=x(2-x-y)      =x'=f(x,y)\] \[\frac{dy}{dt}=y(-1+x)      =y'=g(x,y)\]

Setting these terms to zero:

\[x(2-x-y)=0\] \[y(-1+x)=0\]

leaves us with

\[y=2-x\] \[x=\frac{y}{y}=1\]

This system hence has two trivial steady states \((0, 0)\) and \((2, 0)\) and a non-trivial steady state \((1, 1)\), indicating that the dynamics will not change if either \(x = 0\) and \(y = 0\), or \(x = 2\) and \(y = 0\) or \(x\) and \(y\) both are \(1\).

Alternatively, the non-trivial steady state can be seen in the following plots, indicated as the state (or point) to which the curves seem to run up to, i.e. their attractor

System in time
Phase plane

The left plot is a temporal representation of the system's development, with time \(t\) being represented on the horizontal axis. The right plot is a phase plane (or phase space or state space) portrait of the system. It shows the size of the stock \(x\) on the horizontal axis parametrically with the size of the stock \(y\) on the vertical axis. Each point on this curve refers to one particular state of the system.

A phase plane portrait can provide a global way of looking at a system. You might know the example of the pair of curves \((x(t); y(t)) = (sin t; cos t)\), which , when plotted against time, do not necessarily reveal that they define a circle. If they are plotted parametrically in the plane however, the geometry is clear.

The term phase derives from the representation of phase spaces such as the one of temperature and pressure of water, showing the three phases solid, fluid and gaseous. It would be more correct however, to speak of state spaces since the points in these spaces indicate states of the system.

The stability of steady states

Steady states (or fixed points) can have different stabilities. A population exhibiting limited growth of \(\frac{dN}{dt}=\gamma N\left ( 1-\frac{N}{K} \right )\) for instance has a stable steady state at its carrying capacity \(K\) and an unstable steady state at \(N=0\).

This can be seen if the system is started from a value close to a steady state. Starting it at 0.001 will induce growth. It will drive the development away from zero. Zero works as a repellor in this case. It is an unstable equilibrium.

Starting the same dynamics slightly below or slightly above the carrying capacity \(K\) will lead it back to \(K\). \(K\) works as an attractor. It is a stable equilibrium.

Tangents and derivatives

Analytically, whether a steady state is stable or unstable can be checked by drawing the tangent \(T\) of the system's dynamics close to a steady state. The tangent is generated by finding the derivative \(f'(x)\) of \(f(x)\). In general, the tangent of a function\(f(x)\) at the point \(x=a\) has the form \(f(a)+f'(a)(x-a)\)


Exponential growth of \(N_t= N_0(1 + \gamma)^t\) with its derivative \(\frac{dN}{dt}=(1 + \gamma)^t N ln{(1 + \gamma)}\) and a tangent \(y=N_{50}+N'_{50}(t-50)\) in \(t = 50\), for \(0≀t≀100\)

A table for finding derivatives can be accessed here

To make this a bit more practical, consider the speed of a car on a highway cruising from A to B. To find the car's speed at the point \(x\) of this cruise we need to find the slope of the car’s dynamics at \(x\), which geometrically is the tangent in the point \(x\) of the curve that plots these dynamics.

The second derivative \(f''\) of the car's dynamics in \(x\) would indicate the car's acceleration in point \(x\).

If the tangent has an upward slope (\(f'(t)>0\), like the one in the image above) the steady state or fixed point is unstable - the dynamics are fleeing the fixed point. The fixed point is a repellor.

If the tangent would have a downward slope \((f'(t)<0)\) the steady state (or fixed point) would be stable - the dynamics are attracted by the fixed point. The fixed point is an attractor.


Another possibility to assess the dynamics of a system without simulating it, is to sketch the nullclines \(f(x, y) = 0\) and \(g(x, y) = 0\). Nullclines are curves that allow breaking the phase plane into regions of different qualitative behavior.

The \(x\)-nullcline is the set of points in the plane where \(f(x,y)=0\). This naturally divides the plane into regions where \(f(x,y)>0\) and \(f(x,y)<0\). Similarly, the \(y\)-nullcline is the set of points where \(g(x,y)=0\). It indicates where \(y\) is increasing, decreasing, and remaining constant. Intersections of nullclines are places where \(f(x,y)=0=g(x,y)\), hence steady states.

Stream plot

Phase plane portraits of the system \(f(x,y)=x(2-x-y)\) and \(g(x,y)=y(-1+x)\). Left: its steady states (red dots) and nullclines, right: its stream plot.
A stream plot (or quiver plot in MATLAB) is generated by starting the system from a large number of initial conditions and including all dynamics into one plot.

Geometrically, the \(x\)-nullcline is the set of points where the vectors are either straight up or straight down. Algebraically, the \(x\)-nullcline is found by solving \(f(x,y)=0\).

The \(y\)-nullcline is the set of points where the vectors are horizontal, going either to the left or to the right. Algebraically, the \(y\)-nullcline is found by solving \(g(x,y)=0\).

The directions of vectors vary continuously from one point to the next except at equilibria. A change in orientation (up or down, right or left) can take place only at equilibria.

When a nullcline passes through an equilibrium, the orientation of an  arrow is reversed (if the Jacobian matrix (see below) satisfies \(det(J)\neq 0\)). Local changes in behavior can take place close to equilibria. Thus, it is important to know the local behavior of the system near each equilibrium.

The Jacobian matrix and eigenvalues

The Jacobian matrix \(A\) of the partial derivates of the system {\(f(x,y),g(x,y,)\)} is:

\[A = \left( \begin{array}{ll} \frac{df}{dx} & \frac{df}{dy} \\ \frac{dg}{dx}& \frac{dg}{dy} \\ \end{array} \right) = \left( \begin{array}{ll} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{array} \right) \]

Remember that the determinant of \(A\), \(detA = a_{11}a_{22} - a_{12}a_{21}\) and the trace of \(A\), \(trA=a_{11} + a_{22}\)

The Eigenvalues of \(A\) are \(\lambda = \frac{trA \pm \sqrt{(trA)^2-4detA}}{2}\), of which the real part \((Re)\) informs about stability.

We get a:

  • Sink (stable node): if \(Re (\lambda_1) < 0\) and \(Re (\lambda_2) < 0\)
  • Source (unstable node): if \(Re (\lambda_1) > 0\) and \(Re (\lambda_2) >0\)
  • Saddle (unstable): if \(Re (\lambda_1) < 0\) and \(Re (\lambda_2) > 0\)
  • Spiral or vortex: \(\lambda_1\) and \(\lambda_2\) are complex conjugates: if \(Re (\lambda_1) < 0\) then stable, if \(Re (\lambda_1) > 0\) then unstable.

The above example system

We construct the Jacobian matrix of the system \(\frac{dx}{dt}=x(2-x-y)\) and \(\frac{dy}{dt}=y(-1+x)\)

\[A=\left( \begin{array}{ll} 2-2x-y & -x \\ y & -1+x \\ \end{array} \right) \]

and evaluate them at the fixed points \((0,0), (2,0)\) and \((1,1)\).

For \((0,0)\) we get \(A=\left( \begin{array}{ll} 2 & 0 \\ 0 & -1 \\ \end{array} \right) \) implying that \(detA=-2, TrA=1, \lambda_1=2, \lambda_2=-1\). The fixed point is unstable.

For \((2,0)\) we get \(A=\left( \begin{array}{ll} -2 & -2 \\ 0 & 1 \\ \end{array} \right) \) implying that \(detA=-2, TrA=-1, \lambda_1=-2, \lambda_2=1\). The fixed point is unstable as well.

For \((1,1)\) we get \(A=\left( \begin{array}{ll} -1 & -1 \\ 1 & 0 \\ \end{array} \right) \) implying that \(detA=1, TrA=-1, \lambda_1=-(-1)^{\frac{1}{3}}=-0.5-0.866i, \lambda_2=(-1)^{\frac{2}{3}}=-0.5+0.866i\). The fixed point is stable.

Check out these examples too.