Stability Analysis of Equilibrium Point and Limit Cycle of Two-Dimensional Nonlinear Dynamical Systems—A Tutorial

: The equilibrium state of a dynamical system can be divided into the equilibrium point and limit cycle. In this paper, the stability analysis of the equilibrium point and limit cycle of dynamical systems are presented through different and all possible approaches, and those approaches are compared as well. In particular, the author presented the stability analysis of the equilibrium point through phase plane approach, Lyapunov–LaSalle energy-based approach, and linearization approach, respectively, for two-dimensional nonlinear system, while the stability analysis of the limit cycle is analyzed by using the LaSalle local invariant set theorem and Poincar é –Bendixson theorem, which is only valid in two-dimensional systems. Different case studies are used to demonstrate the stability analysis of equilibrium point and limit cycle.


Introduction
The Lyapunov theorem can be categorized into local theorem and global theorem.The local theorem can be used to prove whether an equilibrium point is stable or asymptotically stable [1][2][3][4][5][6], whereas the global theorem can be used to prove if an equilibrium point is globally asymptotically stable and if it has been used in several applications [7][8][9][10][11][12][13][14][15][16][17].The LaSalle theorem can also be categorized into a local version and a global version, and they are collectively called invariant set theorem [18,19].The purpose of using the LaSalle global invariant set theorem is to further prove that an equilibrium point is globally asymptotically stable, given the fact that the Lyapunov theorem is limited somehow to proving the global stability of an equilibrium point [20][21][22][23].For example, for the mass-spring-damper system, we usually choose the total energy (i.e., the kinetic and potential energy together) of the system as our Lyapunov candidate function, and it is known that if the total energy keeps dissipating, the system will approach an equilibrium point and stop there [24].So, by looking at the time derivative of the chosen Lyapunov function, we can judge if the energy increases or decreases as time goes on.In order to have a stable equilibrium point, we want to make sure that the time derivative of the chosen Lyapunov function at is least negative semi-definite so that the energy decreases as time goes on.If we use the Lyapunov theorem, it can only prove that the mass-spring-damper system has a stable equilibrium point because the time derivative of the Lyapunov function is negative semi-definite but not negative definite.However, the fact of the matter is that the mass-spring-damper system not only has a stable equilibrium point, but also that equilibrium point is asymptotically stable [25].By further using the LaSalle global invariant set theorem, we can then prove that the equilibrium point is asymptotically stable.The LaSalle local invariant set theorem is usually used to prove whether the limit cycle of a dynamical system is stable or not, assuming that there exists a limit cycle.
The Lyapunov and LaSalle theorems are not the only way to prove stability of an equilibrium point [26][27][28][29] and are not the best approaches to use in some cases [30][31][32][33] in order to prove stability of an equilibrium point, given the amount of calculations required and the method not being as intuitive compared to the graphical phase plane method [34][35][36][37].The graphical phase plane method is somehow straightforward as this is a picture-based method and it has been used in many applications [38][39][40][41][42][43][44][45][46][47][48].By looking at the vector field of the state, we can determine how the state evolves as time goes on, i.e., if the vector field goes towards the origin (suppose the equilibrium point is at the origin), then the equilibrium point is asymptotically stable, and otherwise it is not stable.However, the downside of the graphical phase plane method is that it can only be applied to the autonomous system, and it does not apply for non-autonomous systems, as we cannot draw the arrows of the vector field precisely if the system is also a function of time.In the following sections, the author will use two-dimensional system as a demonstration to show how to use the graphical phase plane method to determine the stability of an equilibrium point.
Yet another method to determine if an equilibrium point of a nonlinear dynamical system is stable or not is to linearize around the equilibrium point by using the Taylor formula and to subsequently look at the Jacobian matrix from the linearization and its corresponding eigenvalues, trace and determinant to see if an equilibrium point is stable or not and whether it has been used in different applications [49][50][51][52][53][54][55][56][57][58][59][60].This type of analytical method has limitations, as it can only truly show that the equilibrium point of the real nonlinear system is stable when the equilibrium point is a saddle, node, or spiral.For example, if the equilibrium point is a centre, degenerated node, star, or non-isolated, then the linearized result cannot truly tell you if the equilibrium point is stable or not [61].In the following sections, the author will use a two-dimensional system as a demonstration to show how to use the linearization method to determine the stability of an equilibrium point in detail.
A limit cycle is a closed trajectory in which the adjacent trajectories spiral toward or away from the limit cycle as time goes on [61].Regarding the stability of the limit-cycle of a nonlinear dynamical system, we can use either LaSalle's local invariant set theorem to see if the limit cycle is stable or not [21], or we can use the Poincaré-Bendixson theorem to determine if there exist a stable limit cycle [62].It is noted that the Poincaré-Bendixson theorem is only valid in two-dimensional systems.
The one-dimensional dynamical system is obviously very easy for all audiences to understand, so this study focuses on two-dimensional nonlinear dynamical systems, and the stabilities of the equilibrium point and limit cycle are analyzed.The stability of twodimensional nonlinear dynamical systems contains both the stability of the equilibrium point and the stability of the limit cycle, whereas the stability of one-dimensional nonlinear dynamical systems contains only the stability of the equilibrium point.Furthermore, the two-dimensional nonlinear dynamical systems can have very interesting behaviors, for example, they can have the periodical oscillations, damped oscillations or limit cycle behaviors, where the one-dimensional systems do not have such functions.Also, when we analyze the two-dimensional nonlinear dynamical systems by using the phase plane method, we have two variables to be analyzed as opposed to one variable.In terms of the three and more-than-three dimensional systems, they will lead to chaos and belong to the complex systems [63].The author will not present this topic, particularly in this paper.
This paper is organized in the following order.In Section 2, the stability analysis of the equilibrium points of two-dimensional nonlinear dynamical systems are presented by using the Lyapunov-LaSalle energy-based approach, graphical phase plane approach, and linearizing-around-an-equilibrium-point approach, respectively.Different examples are given.In Section 3, the stability analyses of the limit cycle of two-dimensional nonlinear dynamical systems are presented by using LaSalle local invariant set theorem and the Poincaré-Bendixson theorem, respectively.Comparisons among different approaches are presented in Section 4, and finally Section 5 gives the conclusion.

Stability Analysis by Lyapunov-LaSalle Energy-Based Approach
Before getting into demonstrating how to use the Lyapunov theorem to prove the stability of an equilibrium point of a two-dimensional nonlinear dynamical system, it is necessary to reiterate the Lyapunov theorem here [64,65].
For the Lyapunov local theorem, it states as follows: "for a chosen Lyapunov candidate function V(x ) is positive definite and .

−
) is negative semi definite, the equilibrium point being analyzed is stable.Furthermore, if the .

−
) is negative definite, the equilibrium point is asymptotically stable".It is noted that the underscore represents a vector.
For the Lyapunov global theorem, it states as follows: "for a chosen Lyapunov candidate function ) is positive definite, and .

−
) is negative definite, and V(x − ) tends to infinity as the state tends to infinity, the equilibrium point being analyzed is globally asymptotically stable".
Consider the pendulum in Figure 1, and the dynamical equation of the pendulum is written below:

Stability Analysis by Lyapunov-LaSalle Energy-Based Approach
Before getting into demonstrating how to use the Lyapunov theorem to prove the stability of an equilibrium point of a two-dimensional nonlinear dynamical system, it is necessary to reiterate the Lyapunov theorem here [64,65].
For the Lyapunov local theorem, it states as follows: "for a chosen Lyapunov candi- ( ) V x − tends to infinity as the state tends to infinity, the equilibrium point being analyzed is globally asymptotically stable".
Consider the pendulum in Figure 1, and the dynamical equation of the pendulum is written below: In order to show the stability of the equilibrium point 0 0 function is chosen as below: is not negative definite.
The fact of the matter is that the equilibrium point 0 0  is actually globally asymptotically stable, which can be determined just by observing how the pendulum behaves as time goes on [66,67].As we can imagine, the pendulum will eventually come ally asymptotically stable, which can be determined just by observing how the pendulum behaves as time goes on [66,67].As we can imagine, the pendulum will eventually come back to the equilibrium point, that is to say, the equilibrium point is at least asymptotically stable.However, by using the Lyapunov local theorem, we cannot prove this, which is why the LaSalle theorem comes in.θ = −g • θ, which means the angular acceleration is not zero, except when we reach the equilibrium point at θ = 0. Therefore, the pendulum always has an acceleration, except when it reaches to the equilibrium point and stops there [69].
Hence, the equilibrium point θ As one can see from the above analysis because we made the small-angle assumption, we proved that equilibrium point is globally asymptotically stable.If, however, we do not make this assumption, i.e., sin θ = θ, then we cannot prove that the equilibrium point is globally asymptotically stable by using the LaSalle global invariant set theorem.The reason is that when .. θ = −g • sin θ = 0, this indicates θ = kπ instead of zero.This is the flaw in the LaSalle global invariant set theorem.
In the subsequent section, we will show that by, using the graphical phase plane method, we do not have to make this small-angle assumption in order to prove that the equilibrium point is globally asymptotically stable.In the meantime, the graphical phase plane method is also able to demonstrate how the pendulum behaves when there is a very large angle [61], which cannot be demonstrated by using the analytical Lyapunov-LaSalle approach.
The LaSalle local invariant set theorem states that: "for a chosen Lyapunov candidate function Typically, the LaSalle local invariant set theorem is used to prove whether a limit cycle is stable or not.Consider the following two different examples.
Case 1: consider the following dynamical system, and the dynamical equation is given as follows: . .
where the state is x 2 , and we want to find the equilibrium state (equilibrium point and/or limit cycle if there exists one) of this dynamical system.First of all, we can choose the Lyapunov candidate function of V(x 2 , and by looking at the time derivative of the V(x − ), we can determine how V(x

−
) changes as time goes on.
Appl.Sci.2023, 13, 1136 ) increases as time goes on.At this point, we could not use the LaSalle local invariant set theorem as the conditions of the LaSalle local invariant set theorem are not satisfied, and it might also mean that the limit cycle or the equilibrium point is unstable.In order to determine exactly whether the equilibrium point or limit cycle is stable or not, we will need to choose a different Lyapunov function as follows, V(x If we look at the time derivative of the V(x From here, we can immediately see that the given dynamical system has an equilibrium point x 1 = 0 x 2 = 0 and a limit cycle x 1 2 + x 2 2 = 1.Also, based on the LaSalle global invariant set theorem, we can conclude that the trajectories either tend to the equilibrium point at However, it is not possible to tend to the limit cycle x 1 2 + x 1 2 = 1, as we had the previous assumption that x 1 2 + x 2 2 < 1 and therefore, all the trajectories tend to the equilibrium point at is asymptotically stable.In terms of the limit cycle, based on the LaSalle local invariant set theorem and the Lyapunov candidate function we chose V(x 2 , even though we could not use the LaSalle local invariant set theorem to conclude the result as the conditions of the LaSalle local invariant set theorem are not satisfied.However, we can still use the one of the conditions .

V(x
) increases, which indicates that the limit cycle is unstable.We can also use the computer simulation to verify the result.Now, if we plot the 2 as a function of the state as shown in Figure 2, we can see that all the trajectories start nearby the limit cycle will move away from the limit cycle because .

−
) > 0, as shown in Figure 2c.Please note that the author used the blue arrows in Figure 2c below to indicate how V(x ) changes as time goes on.The valley in Figure 2a is the limit cycle, and it can be seen from the top view, as indicated in the Figure 2b.Again, according to Figure 2, the equilibrium point is stable, as we can see from Figure 2c by the fact that the trajectories (the blue arrows) are approaching the equilibrium point, not moving away from it.The limit cycle is, however, unstable.As we can see from Figure 2c, the trajectories (the blue arrows) are moving away from the limit cycle, not approaching it.
that the author used the blue arrows in Figure 2c below to indicate how ( ) V x − changes as time goes on.The valley in Figure 2a is the limit cycle, and it can be seen from the top view, as indicated in the Figure 2b.Again, according to Figure 2, the equilibrium point is stable, as we can see from Figure 2c by the fact that the trajectories (the blue arrows) are approaching the equilibrium point, not moving away from it.The limit cycle is, however, unstable.As we can see from Figure 2c, the trajectories (the blue arrows) are moving away from the limit cycle, not approaching it.
(a) (b) (c) Case 2: consider the following dynamical system as an example, and the dynamical equation is given as follows: Case 2: consider the following dynamical system as an example, and the dynamical equation is given as follows: . .
where the state is x 2 , and we want to find the equilibrium state (equilibrium point and/or limit cycle if there exists one) of this dynamical system.First of all, we can choose the Lyapunov candidate function as V(x ) changes as time goes on.After calculation, we have below: .
Based on the LaSalle global invariant set theorem, we can conclude that .

−
) = 0 as time goes on, which further indicates that x 1 = 0 x 2 = 0 or x 1 2 + 2x 2 2 = 9 as time goes on.This means, starting from anywhere, the trajectories will tend to Also, from here, we can see that the given dynamical system has an equilibrium point x 1 = 0 x 2 = 0 and a limit cycle x 1 2 + 2x 2 2 = 9.However, the question is that which one of the above (i.e. the equilibrium point or the limit cycle) the trajectories will go to.Just based on the LaSalle global theorem, we cannot further give conclusions.By further using the LaSalle local theorem, we can then conclude that all trajectories will tend to the limit cycle unless the trajectory starts exactly from the equilibrium point at x 1 = 0 x 2 = 0 , which indicates that the equilibrium point at origin is unstable whereas the limit cycle is stable.The proof is as follows: as .

V(x
) decreases as time goes on.So, starting from the nearby limit cycle, all trajectories will come back to the limit cycle.Additionally, starting from nearby origin, all trajectories will move away from the origin.We can also use the computer simulation to verify the result.The plot of the Lyapunov candidate function V(x 2 is illustrated as below in Figure 3.We can see that all the trajectories starting nearby the limit cycle will come back to the limit cycle because .

−
) < 0, as shown in Figure 3c.The valley in Figure 3a is the limit cycle, and it can be seen from the top view, as indicated in the Figure 3b.We can also extend the above method to other dynamical systems.First of all we will need to choose a Lyapunov energy like function ( ) V x − and look at the ( ) and then subsequently make a conclusion about the stability of an equilibrium point.

Stability Analysis by Graphical Phase Plane Approach
Another way to demonstrate the stability of an equilibrium point and a limit cycle is to use the graphical phase plane approach.Here, the graphical phase plane approach will be used to address the same pendulum case as we discussed in the previous section, and the advantages of using such a graphical approach will be shown.
As shown before, we know that the dynamical equation for the pendulum in the previous section is given as Now the state or the input is , and the output is . We will first draw the vector field for the case when  We can also extend the above method to other dynamical systems.First of all we will need to choose a Lyapunov energy like function V(x ) and look at the .

−
) and then subsequently make a conclusion about the stability of an equilibrium point.

Stability Analysis by Graphical Phase Plane Approach
Another way to demonstrate the stability of an equilibrium point and a limit cycle is to use the graphical phase plane approach.Here, the graphical phase plane approach will be used to address the same pendulum case as we discussed in the previous section, and the advantages of using such a graphical approach will be shown.
As shown before, we know that the dynamical equation for the pendulum in the previous section is given as ml 2 ..

θ + b
. θ + mgl sin θ = 0. To generalize the dynamical equation and for the ease of analysis, let us replace θ with x and assume that the length of the pendulum is 1, so the equation of motion is reduced to m Now the state or the input is x y , and the output is .
x .y .We will first draw the vector field for the case when .
x = 0, and then plot the case when x = 0, and finally we will combine the above two cases to get the entire vector field.So, when .
x = 0, .y = −g sin x, the vector field is plotted as follows, as shown in Figure 4. We then draw the vector field for the case when Let us first draw the case when the damping coefficient b is zero, as shown in Figure 5, then we can easily plot the case when b is not zero.x.Let us first draw the case when the damping coefficient b is zero, as shown in Figure 5, then we can easily plot the case when b is not zero.We then draw the vector field for the case when .
Let us first draw the case when the damping coefficient b is zero, as shown in Figure 5, then we can easily plot the case when b is not zero.If we combine the above two plots together, we have the following vector field, as shown in Figure 6, and immediately we can sense the direction of the whole vector field.The above vector field is for the case when the damping coefficient b is zero, as we can see that the vectors form an infinite number of closed orbits, and each of those closed orbits represents how the pendulum behaves as time goes on.Additionally, it can be seen that the trajectories will not move back to the origin but also not move away from the origin, which indicates that the equilibrium point 0 0  is stable but not asymptotically stable.Now when the damping coefficient b is not zero, in this case, all the arrows on the y-axis that are above the origin will point downwards and all the arrows on the yaxis that are below the origin will point upwards, as shown in Figure 7, and the slope of If we combine the above two plots together, we have the following vector field, as shown in Figure 6, and immediately we can sense the direction of the whole vector field.We then draw the vector field for the case when .
Let us first draw the case when the damping coefficient b is zero, as shown in Figure 5, then we can easily plot the case when b is not zero.If we combine the above two plots together, we have the following vector field, as shown in Figure 6, and immediately we can sense the direction of the whole vector field.The above vector field is for the case when the damping coefficient b is zero, as we can see that the vectors form an infinite number of closed orbits, and each of those closed orbits represents how the pendulum behaves as time goes on.Additionally, it can be seen that the trajectories will not move back to the origin but also not move away from the origin, which indicates that the equilibrium point 0 0  is stable but not asymptotically stable.Now when the damping coefficient b is not zero, in this case, all the arrows on the y-axis that are above the origin will point downwards and all the arrows on the yaxis that are below the origin will point upwards, as shown in Figure 7, and the slope of The above vector field is for the case when the damping coefficient b is zero, as we can see that the vectors form an infinite number of closed orbits, and each of those closed orbits represents how the pendulum behaves as time goes on.Additionally, it can be seen that the trajectories will not move back to the origin but also not move away from the origin, which indicates that the equilibrium point x .
x = 0 0 is stable but not asymptotically stable.Now when the damping coefficient b is not zero, in this case, all the arrows on the y-axis that are above the origin will point downwards and all the arrows on the y-axis that are below the origin will point upwards, as shown in Figure 7, and the slope of the downwards and upwards arrows depends on the damping coefficient b.The larger the value of b, the larger the slope will be, and vice versa.In this case, as one can see that when the damping coefficient b is not zero, the trajectories will tend towards the origin as time goes on, which indicates that the equilibrium point x .
On the contrary, if the damping coefficient b is not just nonzero, but also negative, in this case the trajectories will move away from the origin, which means the equilibrium point x .

Stability Analysis by Linearizing around an Equilibrium Point
The third approach to determine the stability of an equilibrium point is to linearize around the equilibrium point to see how its neighbours behave as time goes on [70,71].Here, we will use the previous pendulum case as an example, and we will show the case with damping and the case without damping.This method needs to be used with care as only the following cases can the linearization give qualitatively correct dynamics near the equilibrium point: a saddle, node, and spiral or sometimes referring as focus [61], whereas if the equilibrium point is a centre, star, or non-isolated point, then the linearized result cannot tell you if the equilibrium point is truly stable or not, which is the potential problem of using the linearizing-around-an equilibrium-point method to determine the stability of an equilibrium point.
as time goes on, which indicates that the equilibrium point 0 0 stable.On the contrary, if the damping coefficient b is not just nonzero, but also nega- tive, in this case the trajectories will move away from the origin, which means the equilibrium point 0 0

Stability Analysis by Linearizing around an Equilibrium Point
The third approach to determine the stability of an equilibrium point is to linearize around the equilibrium point to see how its neighbours behave as time goes on [70,71].Here, we will use the previous pendulum case as an example, and we will show the case with damping and the case without damping.This method needs to be used with care as only the following cases can the linearization give qualitatively correct dynamics near the equilibrium point: a saddle, node, and spiral or sometimes referring as focus [61], whereas if the equilibrium point is a centre, star, or non-isolated point, then the linearized result cannot tell you if the equilibrium point is truly stable or not, which is the potential problem of using the linearizing-around-an equilibrium-point method to determine the stability of an equilibrium point.
With damping present, we know that the dynamical equation is .In order to linearize around the equilibrium point 0 0 using the Taylor formula, we will first need to convert the above equation into the following equivalent and assume the length 1 l = for the ease of calculation: Considering a small deviation from the equilibrium point 0 0 With damping present, we know that the dynamical equation is ml 2 ..
In order to linearize around the equilibrium point x − = x .
x = 0 0 using the Taylor formula, we will first need to convert the above equation into the following equivalent and assume the length l = 1 for the ease of calculation: .
Considering a small deviation from the equilibrium point x − = x .
x = 0 0 as follows: In order to look at how the small deviation changes as time goes on, we need to look at the time derivative of u(t) and v(t) as follows: . .
By using the Taylor expansion formula, we can linearize the above equations, which leads to the following: .
We noticed that the higher-order terms are neglected here as they are extremely small compared to the first few terms.By further converting the above equations into a matrix format, we have the following: where J is the Jacobian matrix.As we can see, whether the equilibrium point is stable or not depends on the Jacobian matrix, and particularly depends on the eigenvalues of the Jacobian matrix.If all eigenvalues are negative, that indicates the small deviation will decrease as time goes on, which further indicates that the equilibrium point is stable.To further generalize the whole scenarios including, when eigenvalues are not negative, we can use the trace and determinant of the Jacobian matrix to determine if the equilibrium point is stable or not.By looking at the Jacobian matrix that evaluated at the equilibrium point, i.e., at . As we can see that when b = 0, i.e., no damping, the Jacobian matrix is J = 0 1 −g 0 , and the trace of the Jacobian matrix is 0 and the determinant of the Jacobian matrix is g, which is larger than 0. This indicates that we have a linear centre, and this linear centre indicates that the original nonlinear system is indeed also a centre because when b = 0 the mass-spring system is a conservative system, i.e., there is no energy being dissipated.The following vector field is plot as shown in Figure 8 to demonstrate the centre.The origin represents the angular position and angular velocity of the pendulum both equal to 0. The next circle that surrounds the origin represents the small-angle swing, and the one after the next circle represents the swing with a large angle, and so on and so forth, and this makes sense because there is no damping.
Appl.Sci.2023, 13, 1136 12 of 21 the next circle represents the swing with a large angle, and so on and so forth, and this makes sense because there is no damping.With damping present, i.e., 0 b ≠ , the Jacobian matrix has become the following: , and the trace of the Jacobian matrix is b m − and the determinant of the Jacobian matrix is g , which is greater than 0. Based on the trace and the determinant information, we can conclude that the equilibrium point is stable.
From above analysis, we can see that after linearizing around the equilibrium point, and based on the linearized result, we can conclude whether the equilibrium point of the original nonlinear system is stable or not.The linearization method works for the above case because the equilibrium point of the pendulum with damping is a spiral and without the damping case it is centre but the energy of the system is conserved.However, the linearization method does not work for the cases where the equilibrium point is centre, degenerated node, star, or non-isolated.

Stability of Limit Cycle
The stability of the limit cycle of dynamical systems can be analyzed by using the invariant set theorem and Poincaré-Bendixson theorem.The former is mainly used to prove the stability of a limit cycle and the latter is mainly used to see if there exists a limit cycle.

Stability of Limit Cycle by Invariant Set Theorem
Sometimes the equilibrium state of a dynamical system only has an equilibrium point, and sometimes the equilibrium state of a dynamical system contains both the equilibrium point and the limit cycle [72].In the previous Section 2.1, we touched upon the limit cycle's stability when we analyzed the stability of an equilibrium point.We looked at two different cases, i.e., stable limit cycle and unstable limit cycle.Here we will use a With damping present, i.e., b = 0, the Jacobian matrix has become the following: , and the trace of the Jacobian matrix is − b m and the determinant of the Jacobian matrix is g, which is greater than 0. Based on the trace and the determinant information, we can conclude that the equilibrium point is stable.
From above analysis, we can see that after linearizing around the equilibrium point, and based on the linearized result, we can conclude whether the equilibrium point of the original nonlinear system is stable or not.The linearization method works for the above case because the equilibrium point of the pendulum with damping is a spiral and without the damping case it is centre but the energy of the system is conserved.However, the linearization method does not work for the cases where the equilibrium point is centre, degenerated node, star, or non-isolated.

Stability of Limit Cycle
The stability of the limit cycle of dynamical systems can be analyzed by using the invariant set theorem and Poincaré-Bendixson theorem.The former is mainly used to prove the stability of a limit cycle and the latter is mainly used to see if there exists a limit cycle.

Stability of Limit Cycle by Invariant Set Theorem
Sometimes the equilibrium state of a dynamical system only has an equilibrium point, and sometimes the equilibrium state of a dynamical system contains both the equilibrium point and the limit cycle [72].In the previous Section 2.1, we touched upon the limit cycle's stability when we analyzed the stability of an equilibrium point.We looked at two different cases, i.e., stable limit cycle and unstable limit cycle.Here we will use a different case, i.e., a semi-stable limit cycle, to demonstrate how to use the invariant set theorem to determine the stability of the limit cycle.
Consider the following dynamical system: .
where the state is In order to see the entire equilibrium state of this system, we can use the global invariant set theorem.First of all, we can choose the Lyapunov function as follows: V(x Based on the global invariant set theorem, we can conclude that .

−
) = 0 as time goes on, which indicates that x 1 2 + x 2 2 = 1 or x 1 = 0, x 2 = 0.This means that, starting from outside of the cycle x 1 2 + x 2 2 = 1, all trajectories tend towards the cycle or the origin.If, however, x 1 2 + x 2 2 < 1, geometrically, it means that, if we are inside of the cycle ) increases as time goes on.If we plot the Lyapunov function as a function of the state, as shown in Figure 9, we can see that inside of the cycle, the trajectories tend to the origin if there is a small perturbation.Conversely, outside of the cycle, all trajectories tend to the limit cycle, as shown in Figure 9c.The valley in Figure 9a is the limit cycle, and it can be seen from the top view as indicated in the Figure 9b and therefore, the limit cycle is semi-stable or half stable.Again, from Figure 9c, we can see that the blue arrows are approaching the equilibrium point from two sides, and so the origin is stable; whereas one set of blue arrows are approaching to the limit cycle and the other set of the blue arrows are getting away from the limit cycle, so the limit cycle is semi-stable.
Appl.Sci.2023, 13, 1136 13 of 21 where the state is . In order to see the entire equilibrium state of this system, we can use the global invariant set theorem.First of all, we can choose the Lyapunov function as follows: and If > , geometrically, it means that, if we are outside of the cycle Based on the global invariant set theorem, we can conclude that ( ) 0 as time goes on, which indicates that This means that, starting from outside of the cycle = , all trajectories tend to- wards the cycle or the origin.
If, however, < , geometrically, it means that, if we are inside of the cycle − increases as time goes on.If we plot the Lyapunov function as a function of the state, as shown in Figure 9, we can see that inside of the cycle, the trajectories tend to the origin if there is a small perturbation.Conversely, outside of the cycle, all trajectories tend to the limit cycle, as shown in Figure 9c.The valley in Figure 9a is the limit cycle, and it can be seen from the top view as indicated in the Figure 9b and therefore, the limit cycle is semi-stable or half stable.Again, from Figure 9c, we can see that the blue arrows are approaching the equilibrium point from two sides, and so the origin is stable; whereas one set of blue arrows are approaching to the limit cycle and the other set of the blue arrows are getting away from the limit cycle, so the limit cycle is semi-stable.Now, if we choose the below Lyapunov candidate function, as shown in Equation ( 18), it will not be able to capture the right dynamics of the system and it will not be able to demonstrate the existence of the limit cycle.Now, if we choose the below Lyapunov candidate function, as shown in Equation ( 18), it will not be able to capture the right dynamics of the system and it will not be able to demonstrate the existence of the limit cycle. V(x The .

−
) is always negative.Based on the LaSalle global invariant set theorem, we can conclude that .
V(x − ) = 0 as time goes on, which indicates that x 1 2 + x 2 2 = 1, x 1 = 0, x 2 = 0 as time goes on.However, when we plot the Lyapunov candidate function V(x ), there is no existence of the limit cycle, as shown in the Figure 10.The dashed line is the limit cycle.
The plot shows the 3D view, the side view and the top view.The trajectories are supposed to converge to the limit cycle from outside of the cycle x 1 2 + x 2 2 = 1.However, because of the way we chose the Lyapunov function, which does not capture the limit cycle and therefore, the trajectories do not converge to the limit cycle.So, choosing the right Lyapunov candidate function is critical when we use the Lyapunov-LaSalle energy-based approach.
Appl.Sci.2023, 13, 1136 14 of 21 And The ( ) , there is no existence of the limit cycle, as shown in the Figure 10.The dashed line is the limit cycle.The plot shows the 3D view, the side view and the top view.The trajectories are supposed to converge to the limit cycle from outside of the cycle However, because of the way we chose the Lyapunov function, which does not capture the limit cycle and therefore, the trajectories do not converge to the limit cycle.So, choosing the right Lyapunov candidate function is critical when we use the Lyapunov-LaSalle energy-based approach.

Stability of Limit Cycle by Poincaré-Bendixson (PB) Theorem
It is noted that the PB theorem is mainly used to prove whether there exists a limit cycle instead of proving the stability of the limit cycle, and that the PB theorem is only valid in two-dimensional systems.By using the PB theorem, we are able to determine whether there exists a limit cycle [73], and once there exists a limit cycle, we then need to go back to use the invariant set theorem to show if the limit cycle is stable or not.The Poincaré-Bendixson theorem is described as follows [61]: "If we have the following conditions satisfied: (1) if we can find a closed and bounded region, (2) and there is no equilibrium point in the region, (3) and there exists a trapped trajectory that lies in the region from the beginning and stays in the region for all future times, then the trajectory is either a closed trajectory itself or it spirals toward a closed trajectory as time goes on".
The standard way to find this region is to find a closed region so that the vector field points into the region on its boundary [61, [74][75][76].Consider the following system as an example: .
The question is that can we find a region such that there exists a limit cycle in this region.We can plot the system when .
x = 0 and .y = 0, as shown in Figure 11a.The red and black lines divide the entire region into four different small regions, as shown in Figure 11b.In region 1, we know that .
x > 0 and .y < 0, which indicates that x increases and y decreases, so the vector field in the region 1 points toward the southeast.In region 2, .
x < 0 and .y < 0, which indicates that x decreases and y decreases as well, so the vector field in region 2 points southwest.In region 3, .
x < 0 and .y > 0, which indicates that x decreases and y increases, so the vector field in region 3 point towards northwest.In region 4, .
x > 0 and .y > 0, which indicates that x increases and y increases as well, so the vector field in region 4 points towards northeast.So, overall the vector field in the entire field is in a clockwise direction.By observing the vector field in the entire region, we know that there might be a limit cycle as the vector field forms a cycle.To further verify if indeed there exists a limit cycle, we can construct a closed region, i.e., the region between the right trapezoid and the small circle, as shown in Figure 11b.In order to see if there exists a limit cycle in that region, we need to look at the vector field, and as we already knew that the vector field indeed points into the right trapezoid region on its boundary, all we know to see next is whether the cycle is a repelling cycle or not, i.e., if the vector field in the cycle points into the right trapezoid region or not.We can calculate the Jacobian matrix at the equilibrium point, i.e., the interaction point between and the red line and black line, as follows: At the equilibrium point x y = 0.45 2.31 , the trace of the Jacobian matrix is smaller than 0, which indicates that the equilibrium point is a saddle point.This means that the vector field surrounding the circle does not point into the right trapezoid region and therefore, based on the PB theorem, we cannot conclude that there exists a limit cycle.
At the equilibrium point 0 , the trace of the Jacobian matrix is smaller than 0, which indicates that the equilibrium point is a saddle point.This means that the vector field surrounding the circle does not point into the right trapezoid region and therefore, based on the PB theorem, we cannot conclude that there exists a limit cycle.

Comparison among Different Approaches
Comparisons among the approaches to analyze the stability of the equilibrium point of two-dimensional nonlinear dynamical systems, i.e., Lyapunov-LaSalle energy-based approach, graphical phase plane approach, and linearizing-around-an-equilibrium-point approach, are presented in the Table 1.The advantages and disadvantages of each method are illustrated.Similarly, the comparisons among the approaches to analyze the stability of the limit cycle of two-dimensional nonlinear dynamical system, i.e., Poincaré-Bendixson theorem and LaSalle local invariant set theorem, are also presented in Table 1.

Comparison among Different Approaches
Comparisons among the approaches to analyze the stability of the equilibrium point of two-dimensional nonlinear dynamical systems, i.e., Lyapunov-LaSalle energy-based approach, graphical phase plane approach, and linearizing-around-an-equilibrium-point approach, are presented in the Table 1.The advantages and disadvantages of each method are illustrated.Similarly, the comparisons among the approaches to analyze the stability of the limit cycle of two-dimensional nonlinear dynamical system, i.e., Poincaré-Bendixson theorem and LaSalle local invariant set theorem, are also presented in Table 1.

Advantages Disadvantages
Approaches for stability analysis of a limit cycle --Poincaré-Bendixson theorem 1.The theorem is intuitive, and it is one of the few theorems that can prove whether there exists a limit cycle 1.In order to use this theorem we need to construct a closed region, and it is difficult to construct a closed region that there exists a limit cycle in the region; From the above comparison, we can see that the Lyapunov-LaSalle energy-based analytical method can give quantitative result about the stability of an equilibrium point, and the Lyapunov-LaSalle energy-based approach is applicable to all dynamical systems [77][78][79][80][81][82].However, in order to use this method, we need to choose the right Lyapunov function, and usually it is difficult to choose the right Lyapunov candidate function.If one did not manage finding one, it does not mean the equilibrium point being analyzed is unstable, it could be stable or it could be unstable.Meanwhile, this method involves with large calculations when we calculate the time derivative of the Lyapunov candidate function [83][84][85][86][87][88][89].The graphical phase plane approach is a better choice to analyze the stability of an equilibrium point if quantitative result is not needed [90][91][92], and this approach is a visualization-based method, which can give the conclusion without being getting into calculation, and it is faster than the Lyapunov method when determining the stability of an equilibrium point.However, this method most of the time is limited to the two-dimensional systems and can only be applied to autonomous systems, not non-autonomous systems.Regarding the linearizing-around-an equilibrium-point approach, this involves a very simple calculation because the nonlinear system has been linearized, and this method can give us the qualitatively correct dynamics information quickly near the equilibrium point.Use with care is needed when using this method, as it is not applicable to all dynamical systems.It can only give us the correct dynamics information quickly near the equilibrium point when the equilibrium point is saddle, node, or focus/spiral.
Regarding the comparison among different approaches for stability analysis of a limit cycle, it shows that LaSalle local invariant set theorem is an obvious and easy method to determine if there exists a limit cycle once we determined the time derivative of the Lyapunov candidate function.However, it is difficult to choose the right Lyapunov candidate function.Some literature is available for finding the Lyapunov function [93][94][95].Whereas the Poincaré-Bendixson theorem is intuitive, and it is one of the few theorems [96] that can prove whether there exists a limit cycle, in order to use this theorem we need to construct a closed region, and it is difficult to construct a closed region where there exists a limit cycle in the region, and the method is only valid for two-dimensional systems.

Conclusions
In this paper, the stability analysis of equilibrium point and limit cycle of twodimensional nonlinear dynamical systems are presented.The author summarized that there are three major approaches to analyzing the stability of an equilibrium point, i.e., Lyaponov-LaSalle energy-based approach, graphical phase plane approach, and linearizing-aroundan-equilibrium-point approach.The author also compared among them and presented the advantages and disadvantages of each approach, using several examples to demonstrate each method.In the meantime, the author looked at the stability of a limit cycle by using two different approaches, i.e., LaSalle local invariant set theorem and Poincaré-Bendixson theorem.The advantages and disadvantages of each method are also presented, and different examples are used to demonstrate each method.The author believes that this paper can provide a systematic guideline and overall framework to study the stability of equilibrium points and limit cycles of nonlinear dynamical systems.
As we can see, each method has its drawbacks and limitations.Therefore, for future work, it is interesting to look at the stability analysis of equilibrium points and limit cycles of two-dimensional nonlinear dynamical systems by proposing and using a common approach or framework, for example, how to use phase plane approach to handle the non-autonomous systems.

.
It can be shown that the Lyapunov function is positive definite as
and there exist a value h where V(x − ) is in a bounded region defined by V(x − ) < h, then all trajectories tend to the largest invariant set x stay within this set for all future times".
geometrically it indicates location inside of the circle with a radius 1, then .V(x − ) < 0. Based on the LaSalle global invariant set theorem, we can conclude that .time goes on, which further indicates that

Figure 2 .
Figure 2. Plot of the Lyapunov candidate function of case 1.

Figure 2 .
Figure 2. Plot of the Lyapunov candidate function of case 1.
and by looking at the time derivative of V(x − ), we can determine how V(x −

Figure 3 .
Figure 3. Plot of the Lyapunov candidate function of case 2.

..
To generalize the dynamical equation and for the ease of analysis, let us replace θ with x and assume that the length of the pendulum is 1, so the equation of motion is reduced to The equivalent model for the dynamical equation can be written as below: plot the case when 0 x = , and finally we will combine the above two cases to get the entire vector field.So, when vector field is plotted as follows, as shown in Figure4.

Figure 3 .
Figure 3. Plot of the Lyapunov candidate function of case 2.
m sin x.The equivalent model for the dynamical equation can be written as below: .x = y .y = − b m y − g sin x (6)

Figure 4 .
Figure 4. Vector field when .x = 0. We then draw the vector field for the case when x = 0.When x = 0, .y = − b m

Figure 6 .
Figure 6.Entire vector field when damping coefficient b is zero.

Figure 6 .
Figure 6.Entire vector field when damping coefficient b is zero.

Figure 6 .
Figure 6.Entire vector field when damping coefficient b is zero.

Figure 7 .
Figure 7. Entire vector field when damping coefficient b is not zero.

Figure 7 .
Figure 7. Entire vector field when damping coefficient b is not zero.

Figure 9 .
Figure 9. Plot of the Lyapunov candidate function that shows the limit cycle.

Figure 9 .
Figure 9. Plot of the Lyapunov candidate function that shows the limit cycle.

Figure 10 .
Figure 10.Plot of the Lyapunov candidate function that does not show the limit cycle.Figure 10.Plot of the Lyapunov candidate function that does not show the limit cycle.

Figure 10 .
Figure 10.Plot of the Lyapunov candidate function that does not show the limit cycle.Figure 10.Plot of the Lyapunov candidate function that does not show the limit cycle.

Figure 11 .
Figure 11.Plot of the system when 0 x = 

Table 1 .
Comparison among different approaches for stability analysis of an equilibrium point and limit cycle.

Table 1 .
Comparison among different approaches for stability analysis of an equilibrium point and limit cycle.This method has limitations, it is not applicable to all dynamical systems.It can only give us the correct dynamics information fast near the equilibrium point when the equilibrium point is saddle, node, or focus/spiral.It could give us the false information when the equilibrium point is other than saddle, node, or focus.
2. It is only valid for two-dimensional systems