How to find a Lyapunov function in this case?

Question

We have the system of differential equations $$ \begin{aligned} \frac{dx}{dt} &= y + \sin{x}\\ \frac{dy}{dt} &= -5x-2y. \end{aligned} $$

It's necessary to prove that the system is stable using a Lyapunov function or else show that it's not, following Lyapunov's/Chetaev's theorem.

The first thing I want to ask about is whether it's appropriate to solve the problem for using the fact that $\sin{x} \approx x$ around $x=0$: $$ \begin{aligned} \frac{dx}{dt} &= y + x\\ \frac{dy}{dt} &= -5x-2y. \end{aligned} $$

If that is possible, I would usually check a few functions such as $V(x, y) = ax^2 + by^2$ or $V(x, y) = ax^4 + by^2$, or $V(x, y) = ax^4 + by^4$ or even $V(x,y) = ax^2 + by^2 + cxy$.

The problem is that, unlike simpler problems I haven't yet managed to find such a function that the total derivative of $V$ is strictly positive/negative and the function itself is strictly negative/positive respectively for all pairs of $(x,y)$ except for $(0, 0)$.

Probably, I am trying to prove something that is not true and this is obvious from the beginning. I have tried a few simulation in Python to iterate over different suitable functions and values of $a, b$ to match the criterion, but there was no match.

Answer 1

Taking a left eigenvector of the linearized system and then the real and imaginary parts as coordinate transform suggests to take $$ V=(3x+2y)^2+11x^2=4[5x^2+3xy+y^2] $$ as Lyapunov function. Then $$ \dot V =4[(10x+3y)(x+y)-(3x+2y)(5x+2y)] \\=4[-5x^2-3xy-y^2]=-V $$

Answer 2

For the linearized system $d\mathbf{x}/dt=A \mathbf{x}$, there is an algorithmic way to find Lyapunov functions, namely take any positive definite matrix $Q$ that you like, and solve for $P$ in the Lyapunov equation $A^T P + PA = -Q$. Then if (and only if) the system is asymptotically stable (all eigenvalues of $A$ have negative real part), $P$ will be positive definite, and $V(\mathbf{x}) = \mathbf{x}^T P \, \mathbf{x}$ will be a Lyapunov function with $\dot V(\mathbf{x}) = -\mathbf{x}^T Q \, \mathbf{x}$.

For example, with your $A$ and with $Q=\operatorname{diag}(2,4)$, you get (with the help of Wolfram Alpha) $$ P = \begin{pmatrix} 19 & 4 \\ 4 & 3 \end{pmatrix} , $$ so $V(x,y) = 19 x^2 + 8 xy + 3 y^2$ works; you can check for yourself that it's positive definite and that $\dot V(x,y) = -2 x^2 - 4 y^2$.

The same Lyapunov function works for the original nonlinear system, since $\dot V = -2 x^2 - 4 y^2 + O\bigl( (x^2+y^2)^2 \bigr)$ if you use $\sin x = x + O(x^3)$, so that $\dot V$ is negative definite in a sufficiently small neighbourhood of the origin. If you actually want to use the Lyapunov function to find a domain of stability (which is the whole point, since stability alone is proved in a much easier way from just linearization), then you would of course have to be more careful with keeping track of the size of the remainder in the expansion, so that you can tell how close to the origin you need to be in order to guarantee $\dot V < 0$, but I haven't worked out any details.

Answer 3

As others have already pointed out, one can check the local stability by looking at the linearization around the equilibrium point. Namely, if the system is not locally stable then it also would not be globally stable. Once you have established that your system is locally stable one could try to find a candidate quadratic Lyapunov function by solving a continues Lyapunov equation $A^\top P + P\,A + Q = 0$, with $V(x) = x^\top P\,x$. This should guarantee that near the equilibrium point you have a valid Lyapunov function. However, depending on your choice for $Q$ it might not yield a Lyapunov function with $\dot{V}(x) \leq 0\ \forall\,x\neq0$.

---

A more systematic approach can be obtained using the circle criterion. More detailed steps on how to apply this can be found in another answer of mine, here. It can be noted that the circle criterion is a sufficient but not necessary condition. So if the criterion is not satisfied does not imply the system is unstable.

For example when using the circle criterion on your system allows one to obtain

$$ P = \begin{bmatrix} 2.7160 & 0.8451 \\ 0.8451 & 0.6192 \end{bmatrix}, \tag{1} $$

for $\epsilon = 0.7$, $W = 1.4142$ and $L = \begin{bmatrix}-1.0576 & -0.5976\end{bmatrix}$.

---

Another sufficient way of finding a quadratic Lyapunov function is by formulating the dynamics as $\dot{x} = A(x)\,x$, with $x \in \mathbb{R}^n$ and $A(x) \in \mathbb{R}^{n \times n}$. If the matrix $A(x)$ can also be written as

$$ A(x) = \sum_{i = 1}^N \alpha_i\,A_i, \tag{2} $$

with for each $x$ weights $\alpha_i \in \mathbb{R}$ such that $0 \leq \alpha_i \leq 1$ and $\sum \alpha_i = 1$, then one can find a Lyapunov function by finding a common quadratic Lyapunov function for all $A_i$. In your case the nonlinear term $\sin(x)$ can be bounded by $1\,x$ and $-0.22\,x$, yielding

$$ A_1 = \begin{bmatrix} 1 & 1 \\ -5 & -2 \end{bmatrix}, \quad A_2 = \begin{bmatrix} -0.22 & 1 \\ -5 & -2 \end{bmatrix}. $$

The common quadratic Lyapunov function can be found by solving the following linear matrix inequalities

$$ A_i^\top P + P\,A_i \prec 0,\ \forall\,i \in 1, \cdots, N, \tag{3} $$

for $P \succ 0$. Namely, for the Lyapunov function $V(x) = x^\top P\,x$ one gets

$$ \dot{V}(x) = x^\top \left(A(x)^\top P + P\,A(x)\right) x, \tag{4} $$

which, after substituting in $(2)$, yields

$$ \dot{V}(x) = x^\top \left[\sum_{i = 1}^N \alpha_i \left(A_i^\top P + P\,A_i\right)\right] x. \tag{5} $$

Now by using $(3)$ in $(5)$ yields that the summation in $(5)$ is a sum of negative definite matrices, which always return another negative definite matrix and thus $\dot{V}(x) \leq 0\ \forall\,x\neq0$.

It can be noted that the $P$ from $(1)$ also satisfies $(3)$ using the earlier defined $A_1$ and $A_2$, though a solution could also be obtained using any LMI solver.