Lyapunov Function To Analyze System Stability And Absence Of Closed Orbits
Introduction
In the study of dynamical systems, understanding the stability of solutions is paramount. One powerful tool for analyzing stability is the Lyapunov function. This article delves into the construction of a Lyapunov function for a given system of ordinary differential equations, demonstrating how it can be used to prove the absence of closed orbits and, consequently, the absence of limit cycles. We will explore the theoretical underpinnings of Lyapunov stability theory and then apply it to a specific example, providing a step-by-step guide to constructing a suitable Lyapunov function.
Lyapunov Stability Theory
Lyapunov stability theory provides a framework for assessing the stability of equilibrium points in dynamical systems without explicitly solving the differential equations. At the heart of this theory lies the Lyapunov function, a scalar function that serves as a measure of the "energy" of the system. A Lyapunov function, denoted by V(x), must satisfy certain conditions: it must be positive definite (V(x) > 0 for all x ≠0 and V(0) = 0), and its time derivative along the system's trajectories must be negative semi-definite (dV/dt ≤ 0) or negative definite (dV/dt < 0 for all x ≠0). If dV/dt ≤ 0, the equilibrium point is stable, and if dV/dt < 0, the equilibrium point is asymptotically stable. The intuition behind this is that as the system evolves, the Lyapunov function decreases, indicating that the system's state is moving towards the equilibrium point. Conversely, if we can construct a Lyapunov function whose derivative is positive in some region, it suggests instability in that region. The formal definition of a Lyapunov function hinges on its properties within a neighborhood of an equilibrium point. For a system described by the differential equation dx/dt = f(x), where x is a vector representing the system's state, a Lyapunov function V(x) is a scalar function that satisfies V(0) = 0 and V(x) > 0 for x ≠0 in some neighborhood of the origin. The crucial condition is that the time derivative of V(x) along the trajectories of the system must be non-positive, i.e., dV/dt ≤ 0. This condition implies that the system's state remains bounded and does not move away from the equilibrium point. If, in addition, dV/dt is negative definite (dV/dt < 0 for x ≠0), then the equilibrium point is asymptotically stable, meaning that solutions not only stay close to the equilibrium point but also converge to it as time goes to infinity. Constructing Lyapunov functions is often challenging and may require ingenuity and a good understanding of the system's dynamics. There is no universal method for finding Lyapunov functions, but some common approaches include using the system's energy or a quadratic form. The choice of Lyapunov function depends on the specific system being analyzed, and often involves a process of trial and error. The effectiveness of a particular Lyapunov function is determined by its ability to demonstrate stability or instability over a sufficiently large region of the state space. In some cases, a Lyapunov function may only be valid within a certain domain, and the stability conclusions are only guaranteed within that domain. The absence of a Lyapunov function does not necessarily imply instability; it simply means that this particular method has failed to provide conclusive results. Other methods may still be used to analyze the stability of the system.
System Description
Consider the following system of ordinary differential equations:
dx/dt = x(y^2 + 1) + y
dy/dt = x^2y + x
Our goal is to construct a Lyapunov function to analyze the stability of this system and, more specifically, to show that it has no closed orbits (limit cycles).
Constructing a Lyapunov Function
To construct a Lyapunov function, we need to find a scalar function V(x, y) that satisfies the conditions mentioned earlier. A common starting point is to try a quadratic function of the form:
V(x, y) = ax^2 + by^2
where a and b are positive constants. The choice of a quadratic form is motivated by its simplicity and the fact that it is positive definite (assuming a and b are positive). However, this is just a starting point, and we may need to modify the function based on the system's dynamics. The next step is to compute the time derivative of V(x, y) along the trajectories of the system. This is done using the chain rule:
dV/dt = (∂V/∂x)(dx/dt) + (∂V/∂y)(dy/dt)
Substituting the expressions for dx/dt and dy/dt from the given system, we get:
dV/dt = (2ax)[x(y^2 + 1) + y] + (2by)[x^2y + x]
= 2ax^2(y^2 + 1) + 2axy + 2bx^2y^2 + 2bxy
Now, we want to choose a and b such that dV/dt is negative semi-definite or negative definite. Ideally, we want to eliminate terms that make dV/dt positive. Notice that the terms 2axy and 2bxy have opposite signs depending on the signs of x and y. To make dV/dt negative, we can choose a = -b. However, this would make V(x, y) indefinite, violating the condition that V(x, y) must be positive definite. Therefore, a simple quadratic function might not work in this case. We need to try a different approach. Let's try a function that includes a cross-term:
V(x, y) = ax^2 + by^2 + cxy
where a, b, and c are constants. This form allows for more flexibility in shaping the function to match the system's dynamics. Again, we compute the time derivative:
dV/dt = (∂V/∂x)(dx/dt) + (∂V/∂y)(dy/dt)
= (2ax + cy)[x(y^2 + 1) + y] + (2by + cx)[x^2y + x]
= 2ax^2(y^2 + 1) + 2axy + cxy(y^2 + 1) + cy^2 + 2bx^2y^2 + 2bxy + cx^3y + cx^2
Rearranging the terms, we have:
dV/dt = 2ax^2(y^2 + 1) + 2bx^2y^2 + cx^3y + cx^2 + cxy^3 + (2a + 2b)xy + cy^2
Now, we want to choose a, b, and c such that dV/dt is negative semi-definite. We can try to eliminate some terms by carefully selecting the constants. For example, if we choose c = 0, we are back to the quadratic form we tried earlier, which didn't work. Let's try to eliminate the xy terms. We can set 2a + 2b = 0, which implies a = -b. But again, this would make V(x, y) indefinite. Another approach is to try to make dV/dt a sum of squares, which is always non-positive. This often involves some trial and error and requires a good understanding of the system's dynamics. After some experimentation, let's consider the function:
V(x, y) = x^2 + y^2
This is a simple quadratic function, but let's see if it works. We compute the time derivative:
dV/dt = (∂V/∂x)(dx/dt) + (∂V/∂y)(dy/dt)
= (2x)[x(y^2 + 1) + y] + (2y)[x^2y + x]
= 2x^2(y^2 + 1) + 2xy + 2x^2y^2 + 2xy
= 2x^2y^2 + 2x^2 + 4xy + 2x^2y^2
= 4x^2y^2 + 2x^2 + 4xy
This expression for dV/dt is not obviously negative semi-definite, as it contains both positive and mixed terms. However, it gives us some insight into how we might modify V(x, y) further. The term 4xy is problematic because it can be positive or negative depending on the signs of x and y. To eliminate this term, we might consider adding a term to V(x, y) that cancels out the effect of 4xy in dV/dt. Let's try a different function:
V(x, y) = x^2 - 2xy + 2y^2
This function is positive definite because it can be rewritten as V(x, y) = (x - y)^2 + y^2, which is always non-negative and zero only when x = y = 0. Now, we compute the time derivative:
dV/dt = (∂V/∂x)(dx/dt) + (∂V/∂y)(dy/dt)
= (2x - 2y)[x(y^2 + 1) + y] + (-2x + 4y)[x^2y + x]
= 2x^2(y^2 + 1) + 2xy - 2xy^2 - 2y^2 - 2x^3y - 2x^2 + 4x^2y^2 + 4xy
= 2x^2y^2 + 2x^2 + 2xy - 2xy^2 - 2y^2 - 2x^3y - 2x^2 + 4x^2y^2 + 4xy - 4xy^2 + 4y^2
= 6x^2y^2 - 2x^3y - 6xy^2 + 6xy + 2y^2
This expression for dV/dt is still not obviously negative semi-definite, but it's a step closer. Let's simplify the derivative:
dV/dt = 2x^2(y^2 + 1) + 2xy - 2xy^2 - 2y^2 + (-2x + 4y)(x^2y + x)
= 2x^2y^2 + 2x^2 + 2xy - 2xy^2 - 2y^2 - 2x^3y - 2x^2 + 4x^2y^2 + 4xy
= 6x^2y^2 - 2x^3y - 2xy^2 + 6xy + 2x^2
After careful consideration, let's revisit the function:
V(x, y) = x^2 + y^2
And its derivative:
dV/dt = 2x[x(y^2 + 1) + y] + 2y[x^2y + x]
= 2x^2y^2 + 2x^2 + 2xy + 2x^2y^2 + 2xy
= 4x^2y^2 + 2x^2 + 4xy
We can rewrite dV/dt as:
dV/dt = 2x^2(2y^2 + 1) + 4xy
This form is still not negative semi-definite. Let's try another approach. We need a function that, when differentiated, cancels out positive terms. A key observation is the presence of the term x(y^2 + 1) in dx/dt and x in dy/dt. These terms suggest that we might want to include a term in V(x, y) that involves x^2. Let's consider:
V(x, y) = x^2 + f(y)
where f(y) is some function of y that we need to determine. Now, we compute the time derivative:
dV/dt = (∂V/∂x)(dx/dt) + (∂V/∂y)(dy/dt)
= 2x[x(y^2 + 1) + y] + f'(y)[x^2y + x]
= 2x^2(y^2 + 1) + 2xy + f'(y)x^2y + f'(y)x
We want to choose f'(y) such that dV/dt is negative semi-definite. Let's try to eliminate the 2xy term by setting f'(y) = -2y. This gives us:
f(y) = ∫ -2y dy = -y^2
So, our candidate Lyapunov function becomes:
V(x, y) = x^2 - y^2
However, this function is not positive definite, as it can be negative when |y| > |x|. This means it is not a valid Lyapunov function. Let's try a different approach. Let's add a term that depends on the integral of the terms in dy/dt with respect to x. We have:
dy/dt = x^2y + x
Integrating with respect to x gives us:
∫ (x^2y + x) dx = (1/3)x^3y + (1/2)x^2 + C(y)
Let's consider a Lyapunov function of the form:
V(x, y) = ax^2 + by^2
dV/dt = 2ax[x(y^2+1)+y] + 2by[x^2y+x]
= 2ax^2y^2 + 2ax^2 + 2axy + 2bx^2y^2 + 2bxy
= 2x^2(ay^2+a+by^2) + 2xy(a+b)
If we let a = -b, then we have
dV/dt = 2ax^2y^2+2ax^2-2ax^2y^2 = 2ax^2
This is not negative definite. It seems a simple quadratic Lyapunov function will not work. Let's analyze the system more closely. The nullclines are:
x(y^2 + 1) + y = 0 => x = -y/(y^2 + 1)
x^2y + x = 0 => x(xy + 1) = 0 => x = 0 or y = -1/x
Let's try a function of the form
V(x, y) = x^2 + rac{y^2}{2}
Then
dV/dt = 2x(x(y^2 + 1) + y) + y(x^2y + x)
= 2x^2y^2 + 2x^2 + 2xy + x^2y^2 + xy
= 3x^2y^2 + 2x^2 + 3xy
This is still not negative definite. A key strategy is to consider functions that incorporate the terms present in dx/dt and dy/dt in a way that leads to cancellation or negative terms in dV/dt. After several attempts, there is no traditional Lyapunov function that can prove stability.
Showing No Closed Orbits
To show that the system has no closed orbits (limit cycles), we can use Bendixson's criterion. Bendixson's criterion states that if, in a simply connected region D of the phase plane, the expression:
∂P/∂x + ∂Q/∂y
has a constant sign (either always positive or always negative), where P(x, y) = dx/dt and Q(x, y) = dy/dt, then the system has no closed orbits lying entirely in D.
In our case:
P(x, y) = x(y^2 + 1) + y
Q(x, y) = x^2y + x
We compute the partial derivatives:
∂P/∂x = y^2 + 1
∂Q/∂y = x^2
Thus:
∂P/∂x + ∂Q/∂y = y^2 + 1 + x^2
Since y^2 + 1 + x^2 is always positive for all (x, y), Bendixson's criterion tells us that the system has no closed orbits in the entire phase plane.
Conclusion
In this article, we attempted to construct a Lyapunov function for the given system of ordinary differential equations. While we were unsuccessful in finding a Lyapunov function to prove stability, we were able to use Bendixson's criterion to demonstrate that the system has no closed orbits. This implies that the system does not have any limit cycles. The construction of Lyapunov functions can be challenging, and sometimes other methods, like Bendixson's criterion, provide a more straightforward approach to analyzing the system's behavior. The absence of closed orbits is a significant result, as it provides insights into the long-term dynamics of the system, indicating that solutions do not oscillate indefinitely but either converge to an equilibrium point or diverge to infinity. Understanding the stability properties of dynamical systems is crucial in various fields, including engineering, physics, and biology, where mathematical models are used to describe real-world phenomena. The tools and techniques discussed in this article provide a foundation for further exploration of stability analysis in dynamical systems.
