Artificial Intelligence 🤖
Differential equations
Ordinary differential equations
Systems of linear and nonlinear ODEs
Nonlinear systems in the phase plane

Nonlinear systems in the phase plane

In Sections 14.2 and 14.3, we considered the behaviour of 2D linear systems in the phase plane, x˙=Ax\dot{\boldsymbol{x}}=A \boldsymbol{x} where AA is a real, constant-coefficient matrix. Based on the eigenvalues and eigenvectors of AA, we were able to draw the qualitative behaviour of trajectories in the phase plane and classified the types of fixed points that can occur in 2D linear systems according to the trace and determinant of AA.

27{ }^{27} The subject of bifurcations is not discussed in these notes. For information on bifurcations, refer to the textbook in footnote 23.

Figure 14.16: The trace-determinant plane showing the classification of fixed points for 2D linear systems. Note that the borderline cases are centres, stars/degenerate nodes and lines of fixed points.

We now consider the following nonlinear system:

x˙=f(x),y˙=g(x);\begin{aligned} & \dot{x}=f(\boldsymbol{x}), \\ & \dot{y}=g(\boldsymbol{x}) ; \end{aligned}

where either ff or gg, or both are nonlinear. Our objective here is very similar to what we have been doing so far: given x˙=f(x)\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x}), we want to draw the phase portrait (made up of all the qualitatively different trajectories) and extract information that would allow us to give some physical interpretation for a given system. To do this we need to determine fixed points and assess their stability. This is much harder in 2D2 \mathrm{D} compared to 1D1 \mathrm{D}; the main idea is to find the fixed points, linearise around them and use our knowledge of linear 2D systems to help us draw the phase portrait.

Fixed points in nonlinear systems

Consider the following system,

x˙=f(x,y) and y˙=g(x,y).\dot{x}=f(x, y) \text { and } \dot{y}=g(x, y) .

We look for fixed points at x˙=0\dot{x}=0 and y˙=0\dot{y}=0 or, equivalently,

f(x,y)=0 and g(x,y)=0.f\left(x^{*}, y^{*}\right)=0 \text { and } g\left(x^{*}, y^{*}\right)=0 .

The types of fixed points that occur in nonlinear systems are very similar to those seen in linear systems (see Subsec. 14.3), i.e. we can have saddles, nodes, spirals etc. Suppose that by solving Eqs. (14.98) we find that x=(x,y)\boldsymbol{x}^{*}=\left(x^{*}, y^{*}\right) is a fixed point of the system.

Stability: linearisation around fixed points

Consider small deviations or perturbations u(t)u(t) and v(t)v(t) such that:

x(t)=x+u(t),y(t)=y+v(t),x(t)=x^{*}+u(t), \quad y(t)=y^{*}+v(t),

where uu and vv are small and time dependent. Then, plugging Eq. (14.99) in Eq. (14.97), we have

x˙=f(x+u(t),y+v(t)) and y˙=g(x+u(t),y+v(t)).\dot{x}=f\left(x^{*}+u(t), y^{*}+v(t)\right) \text { and } \dot{y}=g\left(x^{*}+u(t), y^{*}+v(t)\right) .

Since uu and vv are small, we can expand the RHS of x˙\dot{x} and y˙\dot{y} in Eq. (14.100) using a Taylor series expansion (in 2D). For the x˙\dot{x} equation, we have:

x˙f(x,y)+u(t)(fx(x,y))+v(t)(fy(x,y))+ H.O.T \dot{x} \approx f\left(x^{*}, y^{*}\right)+u(t)\left(\frac{\partial f}{\partial x}\left(x^{*}, y^{*}\right)\right)+v(t)\left(\frac{\partial f}{\partial y}\left(x^{*}, y^{*}\right)\right)+\text { H.O.T }

where, H.O.T denotes higher order terms in uu and vv, i.e. uv,u2,v2u v, u^{2}, v^{2}, etc. Further, the first term on the RHS vanishes since ff evaluated at the fixed point is zero while the terms in blue yield a number since these are the partial derivatives of the function ff evaluated at the fixed point. Similarly, for the y˙\dot{y} equation we have,

y˙g(x,y)+u(t)(gx(x,y))+v(t)(gy(x,y))+ H.O.T ,\dot{y} \approx g\left(x^{*}, y^{*}\right)+u(t)\left(\frac{\partial g}{\partial x}\left(x^{*}, y^{*}\right)\right)+v(t)\left(\frac{\partial g}{\partial y}\left(x^{*}, y^{*}\right)\right)+\text { H.O.T },

where, again, g(x,y)=0g\left(x^{*}, y^{*}\right)=0 and the terms in blue are numbers. From Eq. (14.99), differentiating with respect to time gives,

x˙=u˙,y˙=v˙.\dot{x}=\dot{u}, \quad \dot{y}=\dot{v} .

Using Eqs. (14.101) - (14.103), setting f(x,y)=g(x,y)=0f\left(x^{*}, y^{*}\right)=g\left(x^{*}, y^{*}\right)=0 and, ignoring H.O.T, we have the following 2 differential equations for the time evolution of the perturbations uu and vv

u˙=u(fx(x,y))+v(fy(x,y)),v˙=u(gx(x,y))+v(gy(x,y)).\begin{aligned} & \dot{u}=u\left(\frac{\partial f}{\partial x}\left(x^{*}, y^{*}\right)\right)+v\left(\frac{\partial f}{\partial y}\left(x^{*}, y^{*}\right)\right), \\ & \dot{v}=u\left(\frac{\partial g}{\partial x}\left(x^{*}, y^{*}\right)\right)+v\left(\frac{\partial g}{\partial y}\left(x^{*}, y^{*}\right)\right) . \end{aligned}

In matrix form, Eqs. (14.104), are given by

(u˙v˙)=(fx(x,y)fy(x,y)gx(x,y)gy(x,y))(uv),\left(\begin{array}{c} \dot{u} \\ \dot{v} \end{array}\right)=\left(\begin{array}{cc} f_{x}\left(x^{*}, y^{*}\right) & f_{y}\left(x^{*}, y^{*}\right) \\ g_{x}\left(x^{*}, y^{*}\right) & g_{y}\left(x^{*}, y^{*}\right) \end{array}\right)\left(\begin{array}{c} u \\ v \end{array}\right),

where the matrix in green is known as the Jacobian of the system, evaluated at the fixed point (x,y)\left(x^{*}, y^{*}\right). Note that the subscripts ' xx ' and ' yy ' denote partial differentiation with respect to xx and yy, respectively. So, now, what we have managed to do is to express the equations for the perturbations uu and vv near the fixed point (x,y)\left(x^{*}, y^{*}\right) as a linear, constant-coefficient system. It follows that the perturbations satisfy this linear system. From the system given in (14.105), we can assess the stability of the fixed points found using Eq. (14.98) simply by considering the eigenvalues and eigenvectors of the Jacobian of the system and using our knowledge on linear 2D systems. Let us look at an example of this next.

Example 14.1 Consider the following system

x˙=x+x3,y˙=2y\dot{x}=-x+x^{3}, \dot{y}=-2 y

Find and classify the fixed points.

Solution \quad Fixed points occur at x˙=0\dot{x}=0 and y˙=0\dot{y}=0, i.e.

x(x21)=0,y=0.x\left(x^{2}-1\right)=0, y=0 .

Equation (14.107) gives (0,0),(1,0)(0,0),(1,0) and (1,0)(-1,0) as the fixed points. To assess the stability, we need to use the Jacobian of the system which is computed as:

J=(1+3x2002)J=\left(\begin{array}{cc} -1+3 x^{2} & 0 \\ 0 & -2 \end{array}\right)

Next, we evaluate the Jacobian at each fixed point and consult our trace-determinant plane (see Fig. 14.16) to classify it. At (0,0), we have

J=(1002)J=\left(\begin{array}{cc} -1 & 0 \\ 0 & -2 \end{array}\right)

the eigenvalues are λ1=1\lambda_{1}=-1 and λ2=2\lambda_{2}=-2 and so (0,0)(0,0) is a stable node. Similarly, by evaluating the Jacobian at the other 2 fixed points we find that both (1,0)(1,0) and (1,0)(-1,0) are saddles.

Example 14.1 shows how linearising around the location of the fixed points helps us classify the type of fixed points occurring in the nonlinear 2D system.

Structural stability

A natural question is the following: by ignoring the H.O.T in Eq. (14.105), does the resulting linear system predict the stability of the fixed points correctly? It turns out that whether or not we trust the result specified by the eigenvalues of the linearisation, depends on the prediction itself. If the linearisation predicts saddles, spirals or nodes, then we can trust the result. But, if the linearisation predicts one of the so-called borderline cases, i.e. stars, degenerate nodes, centres or non-isolated fixed points, then it is possible that the nonlinear terms that we ignored in our linearisation (i.e. the H.O.T) may become important and affect the results. It helps to think about where each class of fixed points lies in the trace-determinant plane. In the case of saddles, spirals and nodes, the fixed points are determined by strict inequalities represented by entire regions in the trace-determinant plane (i.e. Δ<0\Delta<0 for saddles). A small perturbation of these strict inequalities (such as that represented by u2u^{2} or v2v^{2} or uvu v which are smaller than uu and vv ) does not change the inequalities. On the contrary, the borderline cases, as their name suggests, fall exactly on some curve (e.g. τ24Δ=0\tau^{2}-4 \Delta=0 for degenerate nodes or stars). It follows, that it is possible that small perturbations may affect the results.

These results can be summarised in terms of the eigenvalues of the Jacobian evaluated at the fixed points: the fixed point x\boldsymbol{x}^{*} is said to be hyperbolic if the real part of all the eigenvalues of the linearisation (these are the eigenvalues of the Jacobian matrix evaluated at x\boldsymbol{x}^{*} ) have nonzero real parts, i.e. Re(λi)0\operatorname{Re}\left(\lambda_{i}\right) \neq 0 for both eigenvalues; here, i=1,2i=1,2 (but this can also apply to higher order systems as well, i.e. i=1,,n)i=1, \ldots, n). The fixed point is nonhyperbolic if Re(λi)=0\operatorname{Re}\left(\lambda_{i}\right)=0 for some ii. Hyperbolic fixed points are said to be robust: they are structurally stable and unaffected by small perturbations or nonlinearities and their stability is correctly predicted by the linearisation. Nonhyperbolic fixed points are generally not structurally stable and to determine their stability correctly, one needs to consider the nonlinear terms. Going back to the trace-determinant plane and considering the aforementioned structurally unstable cases, we can see that a predicted linear centre lies on the Δ>0,τ=0\Delta>0, \tau=0 line. Small perturbations may push the point up or down resulting in an unstable or stable spiral, respectively. In degenerate cases, the points lie on the τ24Δ=0\tau^{2}-4 \Delta=0 curve; if τ>0\tau>0, small perturbations may result in the fixed point switching from an unstable degenerate node (or unstable star) to an unstable node or unstable spiral. Similarly, if τ<0\tau<0, the fixed point may switch from a stable degenerate node (or stable star) to an stable node or stable spiral. The Hartman-Grobman theorem essentially states that the local phase portrait near a hyperbolic fixed point is successfully predicted by the linearisation.

Out of all the borderline cases on the trace-determinant plane, the most important one by far is the case of centres, i.e. where τ=0\tau=0 and Δ>0\Delta>0. The reason why is twofold. First, in the other borderline cases, if the predicted linear result is affected by small perturbations, the stability does not change; for instance consider the stars or degenerate nodes that lie on the τ24Δ=0\tau^{2}-4 \Delta=0 line where small perturbations can make the point switch to either an unstable spiral or an unstable node (provided that τ>0\tau>0 ) or, either a stable spiral or a spiral node (provided that τ<0\tau<0 ). In those cases however, the stability remains the same: i.e. an unstable degenerate node cannot switch to a stable fixed point given small perturbations. Second, trajectories given by centres represent oscillatory motion; it is important to know whether certain parameters can cause the system to exhibit oscillations; for instance, large, stable oscillations might be undesirable such as those observed in aircraft systems. In 2D systems, there exist cases where there is more to a phase portrait than just fixed points; it is possible to have structures known as limit cycles. The topic is beyond the scope of the course but we give a brief introduction here: limit cycles are isolated closed trajectories that do not begin or end at a fixed point and their existence is a nonlinear phenomenon only. These are different to the centres since neighbouring trajectories are not closed: they either spiral away (limit cycle is unstable) or toward a limit cycle (stable). While limit cycles can be stable or unstable, the stable one is much more commonly encountered in applications. Considering the stable limit cycle, the solution represents a periodic orbit that the system is trying to achieve. Stable limit cycles are more robust than centres. Take a region in the phase plane where the solution is described by concentric periodic orbits surrounding a centre (see Fig. 14.10). Each one of the orbits is neutrally stable: the solution oscillates with the same amplitude forever. Suppose that given initial conditions, the solution to our system is represented by one of those orbits. A disturbance can push us up (or down) to another orbit where, again, we move on a periodic trajectory but the system is now neutrally stable at this new orbit: the dynamics do not act to bring the system back to its old state \Rightarrow perturbations persist. Contrast this with a stable limit cycle representing a periodic orbit. The system forgets the effect of a perturbation causing us to come off the cycle, since the trajectories near the limit cycle (inside and out) always want to approach it. The various limit cycles encountered in 2D nonlinear systems are shown in Fig. 14.17. Panel (a) shows a stable limit cycle where all neighbouring trajectories approach it. Panel (b) shows an unstable limit cycles where neighbouring trajectories move away from it. Panel (c) shows a semistable limit cycle; these are present in bifurcations where a change in system parameters can trigger qualitative changes in the behaviour of the system.

(a)

(b)

(c)

Figure 14.17: (a) Stable limit cycle: trajectories from the inside and outside of the cycle approach the limit cycle. (b) Unstable limit cycle: all trajectories move away from the limit cycle. (c) Semi-stable limit cycle: trajectories approach the limit cycle from one direction and move away from it from the other.

Generally, to prove the existence or absence of limit cycles in a 2D2 \mathrm{D} nonlinear systems is not a trivial task. Simply studying the local behaviour of fixed points (i.e. through linearisation as shown above) does not offer enough information to conclude that a limit cycle does/does not exist. This raises an important question in dynamical systems in 2D: given a nonlinear system, do limit cycles exist? Additionally, suppose the linearisation results predict a linear centre for a fixed point, are there situations for which the fixed point really is a centre? We can easily think of physical systems for which the last statement isi s true. For example, we expect a system of equations describing the nonlinear pendulum in the absence of damping to admit (undamped) oscillatory solutions; the latter are represented by centres and therefore centres can be observed in physical systems. It can be shown that the predicted linear centre is robust, i.e. it exists in nonlinear systems, if the system is conservative e.g. in systems where we have conservation of energy.

Application: competition models

We now consider the Lotka-Volterra model for the dynamics of 2 population species competing for shared resources. We let x(t)x(t) denote the size of the rabbit population and y(t)y(t) be the size of the sheep population. We start off with a logistic growth for both species, i.e.

x˙=x(3x),y˙=y(2y).\dot{x}=x(3-x), \quad \dot{y}=y(2-y) .

We then add a death rate to each equation; note that this represents death due to depletion of the resources because the food supplies are being consumed by the competing species and not because the rabbits eat the sheep or vice-versa. This results in the following nonlinear competition model,

x˙=x(3x)2xy,y˙=y(2y)xy.\dot{x}=x(3-x)-2 x y, \quad \dot{y}=y(2-y)-x y .

To analyse the system, we proceed in the usual way: we find and classify fixed points. The fixed points occur at x˙=0\dot{\boldsymbol{x}}=0 and are calculated as

  • (0,0)(0,0) - extinction of both species;
  • (0,2)(0,2) - no rabbits;
  • (3,0)(3,0) - no sheep;
  • (1,1)(1,1) - coexistence.

To classify their stability, we compute the Jacobian of the system as follows,

J=(32x2y2xy2x2y).J=\left(\begin{array}{cc} 3-2 x-2 y & -2 x \\ -y & 2-x-2 y \end{array}\right) .

Evaluating the Jacobian at each fixed point gives the behaviour locally, around the fixed point. This reduces the nonlinear problem to a linear one and we can draw the phase portrait close to the fixed point as if it behaves like a linear system. As indicated in Subsec. 14.4.1, we can trust the linearisation for hyperbolic fixed points (saddles, nodes and spirals) and we need to be careful with the rest. For this example, at (0,0)(0,0), we have an unstable node (since the eigenvalues of the Jacobian are calculated as λ1=3\lambda_{1}=3 and λ2=2\lambda_{2}=2 ). This is a hyperbolic fixed point since λ10\lambda_{1} \neq 0 and λ20\lambda_{2} \neq 0. Let us examine the behaviour of the trajectories close to the fixed point.

The corresponding eigenvectors are (10)\left(\begin{array}{ll}1 & 0\end{array}\right)^{\top} and (01)\left(\begin{array}{ll}0 & 1\end{array}\right)^{\top} for λ1\lambda_{1} and λ2\lambda_{2}, respectively. These give us the half (or straight) solutions and their corresponding direction. Note that these straight solutions are only valid near the fixed point. Around the origin therefore, all trajectories move away from the fixed point and, as mentioned before, they do so tangent to the slow direction. This information is shown in Fig. 14.18 where the boxes around each fixed point represent regions wherein the equilibrium point linearisation is valid. Outside of these regions, nonlinear effects become important. At (3,0)(3,0) and (0,2)(0,2) we have stable nodes and all nearby trajectories move toward these stable points. Double and single arrowheads on the half solutions represent fast and slow directions, respectively. Again, nearby trajectories approach the fixed point tangent to the slow directions. Finally, at (1,1)(1,1), we have a saddle point and nearby trajectories (within the 'linear' box) look like hyperbolas moving away from the point. Trajectories around the saddle point approach one of the stable points depending on initial conditions. For example, trajectories below y=1y=1 and above some threshold x=xcx=x_{c}, want to approach the stable point (3,0)(3,0). The latter represents a state where we have only rabbits surviving. This makes sense: if we start with a low initial sheep population (i.e. y<1y<1 ) and a not-so-low rabbit population x>xcx>x_{c}, then the rabbit population dominates thus approaching the equilibrium state of (3,0)(3,0). We have a similar argument for the sheep. If x<1x<1 and y>ycy>y_{c} (=some threshold in yy ) then the sheep population dominates. The exact values of xcx_{c} and ycy_{c} may be determined by solving for the curves that pass through the saddle point; these curves are referred to as the separatrices of the saddle point and are shown in Fig. 14.19). Now that we understand the local behaviour of the fixed points by carrying out the linearisation, we can complete the phase portrait. The phase portrait generated by pplane (the pplane java software can be downloaded here (opens in a new tab)) is shown in allows the solution of system of ODEs in the phase plane and plots the phase portraits indicating location and stability of fixed points.

Figure 14.18: Linear behaviour around fixed points. The boxes around each fixed point represent regions wherein the linearisation captures the behaviour of the system adequately. The fixed points are shown at (0,0),(3,0),(0,2)(0,0),(3,0),(0,2) and (1,1)(1,1). Lines with double and single arrowheads represent the fast and slow half solutions, respectively. Dotted lines represent half solutions that have no physical significance (i.e. negative numbers for rabbits and sheep!).

Figure 14.19: The separatrices of the saddle point. These curves separate the different qualitative behaviour observed for the nonlinear system given by Eq. (14.111) (see also Fig. 14.20).

Fig. 14.20. We observe that coexistence, represented by (1,1)(1,1) is a highly unlikely scenario. To approach this equilibrium state where we have both rabbits and sheep existing in the long-term we would have to be exactly on the stable manifold of the saddle point. As the saddle point is unstable, small fluctuations (i.e. one rabbit dying or one sheep dying) immediately lead to the extinction of one species and survival of the other; these are equilibria represented by the 2 stable states: (3,0)(3,0) and (0,2)(0,2).

Figure 14.20: The phase portrait and direction field for the nonlinear system (14.111).

Equations of trajectories

In some cases, it is possible to obtain explicit equations for trajectories of coupled systems. Consider for example a simpler competition model that the one in Eq. (14.111),

x˙=x(2y),y˙=y(x1).\dot{x}=x(2-y), \quad \dot{y}=y(x-1) .

Since the negative competition term (i.e. xy-x y ) appears in the xx equation, we conclude that xx represents the prey and yy is its predator. Solving for fixed points we find that there are two: (0,0)(0,0) and (1,2)(1,2). To assess their stability, we compute the Jacobian of the system,

J=(2yxyx1).J=\left(\begin{array}{cc} 2-y & -x \\ y & x-1 \end{array}\right) .

Evaluated at the origin (14.114) yields Δ<0\Delta<0 implying that (0,0)(0,0) is a saddle. Evaluated at (1,2)(1,2), we have τ=0\tau=0 and Δ>0\Delta>0 which represents a centre. Now as alluded to above, in nonlinear systems, a predicted linear centre is not a trustworthy result. However, as mentioned above in the section on structural stability, there are certain nonlinear systems which do exhibit periodic behaviour. While this is generally nontrivial to prove, in this case we can proceed as follows to show that the system actually admits periodic solutions. Recall that we are looking to plot the phase portrait in the xyx-y plane. Dividing the yy equation by the xx equation, we get

dydx=y(x1)x(2y)\frac{d y}{d x}=\frac{y(x-1)}{x(2-y)}

which we recognise as a first order, separable ODE. Upon integration, for x>0,y>0x>0, y>0, we have

2lnyy+lnxx=c,2 \ln y-y+\ln x-x=c,

where cc is some constant. For different values of cc, Eq. (14.116) represents a family of curves in the xyx-y plane; see Fig. 14.21. Specifically, Fig. 14.21 shows contours of the function z(x,y)=2lnyy+lnxxz(x, y)=2 \ln y-y+\ln x-x from which we can see that the solution [Eq. 14.116] to the nonlinear system (14.113) admits periodic solutions characterised by the closed curves.

Figure 14.21: Some contours of the function z(x,y)=2lnyy+lnxxz(x, y)=2 \ln y-y+\ln x-x are shown. This shows that the solution (14.116) is represented by closed curves and hence the predicted centre in the nonlinear system (14.113) persists.