# How to write that Lyapunov function?

For didactic purpose I investigate the stability of the ODEs system {x'[t] == x[t] + y[t], y'[t] == -5 x[t] - 2 y[t]}.

Its stream plot

StreamPlot[{x + y, -5 x - 2 y}, {x, -10, 10}, {y, -10, 10}] clearly shows a spiral sink at the origin.

How to prove it, making use of Lyapunov function? First of all, how to explicitly write down that Lyapunov function for the system under consideration with help of LyapunovSolve? I don't find it here. The first example in the Applications section of the documentation to LyapunovSolve is useful, but it produces the conclusion, not a Lyapunov function.

Maybe this?:

ClearAll[v, x, y];
a = CoefficientArrays[{x + y, -5 x - 2 y}, {x, y}] // Last;
p = LyapunovSolve[Transpose[a], -IdentityMatrix];
v[x_, y_] =  p . {x, y} . {x, y} // Simplify (* Lyapunov function *)

(* Check "energy dissipation" dv/dt *)
Dt[v[x, y]] /. Thread[Dt[{x, y}] -> {x + y, -5 x - 2 y}] // Simplify

(*  1/6 (32 x^2 + 14 x y + 5 y^2)  -- Lyapunov function *)

(*  -x^2 - y^2  -- dv/dt is neg. def. *)

• As usually, very good. Jun 2, 2021 at 19:05
• @Michael E2 could you do this to find a Lyapunov function for "any" system?
– Math
Mar 17, 2022 at 14:47
• @Math I doubt it. It fails on this: {x'[t] == x[t]^3 + y[t], y'[t] == -5 x[t] - 2 y[t]^3}. But it succeeds here: ode = {x'[t] == Sin[x[t]] + y[t], y'[t] == -5 x[t] - 2 Sin[y[t]]}; vf = First@SolveValues[ode, {x'[t], y'[t]}]; a = D[vf, {{x[t], y[t]}}] /. {x[t] -> 0, y[t] -> 0}; p = LyapunovSolve[Transpose[a], -IdentityMatrix]; v[x_, y_] = p . {x, y} . {x, y} // Simplify. I don't think there's any guarantee that there exists a quadratic Lyapunov function. Mar 17, 2022 at 16:51
• @MichaelE2 Thank you! Can you please have a look at this question: mathematica.stackexchange.com/questions/265288/…
– Math
Mar 18, 2022 at 11:35

Stability of {x',y'}]== f[x,y] can be investigated by searching for critical points of f (see e.g: https://en.wikipedia.org/wiki/Stability_theory). In your case this is obvious {0,0}. Then in the general case f is linearized at the critical points, again in your case that is not needed, as f is already linear. Then the Eigensystem is

mat= m = {{1, 1}, {-5, -2}};
{x'[t], y'[t]} = mat . {x[t], y[t]}
{eigenval, eigenvec}=Eigensystem[{{1, 1}, {-5, -2}}] // N

{-0.5 + 1.65831 I, -0.5 - 1.65831 I}, {{-0.3 - 0.331662 I,
1.}, {-0.3 + 0.331662 I, 1.}}}


The stability depends on the real parts of the eigenvalues. In your case mat is real, what makes the 2 eigenvalues conjugate complex and there is only one real part of -0.5.

If the real part is negative, as in your case, we have an attractive critical point, called a sink. As the imaginary part is different from zero, it is a spiral sink. The trajectory will approach zero the faster, the larger the absolute value of the eigenvalue.

• Sorry, this is not it. The question is about the construction of Lyapunov function with help of Mathematica. Hope I am and was clear. Jun 2, 2021 at 18:40
• You have a very simple problem that does not need the Liapunov coefficient. This coefficient is needed if you want to investigate the rate of separation of infinitesimally close trajectories. .The Liapunov theory is more evolved than the simple stability analysis. Jun 2, 2021 at 18:47