Strong Lyapunov Functions for Stability

• Strong Lyapunov functions are refined constructs that guarantee not only stability but also precise decay rates and strict descent outside recurrent sets.
• They are applied in continuous systems, differential inclusions, PDEs, and optimization, using rate-coded inequalities and computational methods like SOS programming.
• Practical implementations include certifying exponential or algebraic decay rates and ensuring safety in control systems via a single Lyapunov–barrier formulation.

Strong Lyapunov functions are Lyapunov constructions that do more than certify qualitative stability. Depending on the setting, they encode exact convergence rates, strict decay outside recurrent sets, dissipation inequalities that hold for every solution of a differential inclusion, or a single certificate for both safety and asymptotic stability. This suggests that the term does not denote one universal definition, but a family of strengthened Lyapunov notions in which derivative inequalities, level-set geometry, or global dynamical information are sharpened beyond the classical requirement $\dot V\le 0$ (Jagt et al., 4 Jan 2026, Giesl et al., 2021, Adly et al., 2018, Bernardi et al., 2017, Quartz et al., 15 Sep 2025).

1. Terminological scope and unifying themes

A recurring theme is the replacement of a merely qualitative statement—such as asymptotic stability or monotonic decay—by a quantitative or structural one. In continuous optimization, the “strong Lyapunov condition” is the inequality

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$

which yields exponential decay when $q=1$ and algebraic decay when $q>1$ (Chen et al., 2021). In differential inclusions, a strong $a$-Lyapunov pair $(V,W)$ requires

$$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$

for every solution, not merely along selected trajectories (Adly et al., 2018). In strong chain recurrence theory, a strong Lyapunov function is a Lipschitz or continuous Lyapunov function that is strict outside the strong chain recurrent set $\mathcal{SCR}(\phi)$ or $\mathcal{SCR}_d(\phi)$ (Bernardi et al., 2017, Bernardi et al., 2020).

A plausible common denominator is that “strong” refers to one or more of four upgrades. The first is strictness, meaning $\dot V<0$ or $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$0 away from the target set. The second is quantitative rate coding, where the derivative of $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$1 reproduces a prescribed decay law. The third is global dynamical encoding, where the function separates recurrent and gradient-like regions of phase space. The fourth is universal quantification, where the Lyapunov inequality must hold for every admissible trajectory or successor state.

This family resemblance is visible across control, dynamical systems, nonsmooth analysis, PDEs, and optimization. The resulting notions are not interchangeable, but they are structurally aligned: each sharpens Lyapunov theory from mere existence of decrease to decrease with rate, geometry, robustness, or exact set representation.

2. Rate-coded strong Lyapunov functions for nonlinear ordinary differential equations

The most explicit rate-based formulation is developed for nonlinear autonomous ODEs

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$2

through the notion of $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$3-stability. With a continuous state measure $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$4 and a normalized, time-invariant $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$5, the system is $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$6-stable on $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$7 with rate $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$8 and gain $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$9 if

$q=1$0

Here the rate performance is the largest $q=1$1 for which the bound holds for some $q=1$2 (Jagt et al., 4 Jan 2026).

The central result is a converse Lyapunov theorem showing that $q=1$3-stability with a prescribed rate $q=1$4 is equivalent to the existence of Lyapunov functions satisfying inequalities parametrized by the scalar comparison law $q=1$5, where $q=1$6. For every $q=1$7, there exists $q=1$8 such that

$q=1$9

and

$q>1$0

almost everywhere, or the corresponding Dini derivative inequality in the finite-time setting. This is not merely sufficient: it is necessary and sufficient for the exact rate notion encoded by $q>1$1 (Jagt et al., 4 Jan 2026).

Three standard decay laws appear as special cases. For exponential decay,

$q>1$2

and the strong Lyapunov inequality becomes

$q>1$3

For rational decay,

$q>1$4

giving

$q>1$5

For finite-time stability,

$q>1$6

and the inequality becomes

$q>1$7

These are “strong” because the derivative of $q>1$8 is rate-coded and gain-normalized, and because the corresponding inequalities characterize the exact achievable rate, not merely asymptotic convergence (Jagt et al., 4 Jan 2026).

The same framework yields computational certificates. For polynomial vector fields, the inequalities can be enforced by SOS programming, and the numerical examples show that the SOS tests can achieve tight bounds on rate performance with accurate inner bounds on associated regions of performance. In the exponential case, the paper reports an SOS-certified rate $q>1$9 for a modified Lorenz system, while direct simulation gives $a$0 (Jagt et al., 4 Jan 2026). This suggests that, in this setting, a strong Lyapunov function is one whose derivative reproduces the exact scalar comparison dynamics $a$1.

3. Global and Conley-type strong Lyapunov functions

A different usage of “strong” arises in global dynamical systems, where Lyapunov functions are required to organize the full recurrent structure of the flow. For an autonomous ODE

$a$2

Giesl, Hafstein, and Suhr study complete Lyapunov functions $a$3. Such a function is $a$4, satisfies

$a$5

is strictly decreasing outside the chain-recurrent set $a$6,

$a$7

has $a$8 on $a$9, and its critical values separate chain transitive components (Giesl et al., 2021).

The same work proves that the orbital derivative can be prescribed on any compact set $(V,W)$0. If $(V,W)$1 is $(V,W)$2, then there exists a complete $(V,W)$3 Lyapunov function $(V,W)$4 such that

$(V,W)$5

This sharpens the complete Lyapunov notion from qualitative monotonicity to quantitative shaping of the descent rate on selected compact subsets (Giesl et al., 2021).

A related metric refinement replaces Conley’s chain recurrence by strong chain recurrence. For a continuous flow on a compact metric space, there exists a Lipschitz continuous Lyapunov function $(V,W)$6 such that

$(V,W)$7

equivalently,

$(V,W)$8

Moreover,

$(V,W)$9

where $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$0 is the set of all Lipschitz continuous Lyapunov functions for the flow (Bernardi et al., 2017). Bernardi and Florio establish the analogous Conley-type result for $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$1, constructively proving the existence of a continuous Lyapunov function strictly decreasing outside the strong chain recurrent set and expressing $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$2 as an intersection of unions $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$3 over strongly stable sets $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$4 (Bernardi et al., 2020).

In smooth flow theory, quadratic infinitesimal Lyapunov functions play a comparable role. For star flows, a $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$5 vector field is characterized by the existence of a field of quadratic forms $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$6 such that the flow is strictly $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$7-separated on the preperiodic set, the linear Poincaré flow is strictly $$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$8-monotone, and at singularities

$$e^{a t} V\bigl(x(t;x_0)\bigr) + \int_0^t W\bigl(x(\tau;x_0)\bigr)\,d\tau \le V(x_0)$$9

These cone and infinitesimal monotonicity conditions encode strong homogeneity, star property, and singular hyperbolicity (Salgado, 2017).

Taken together, these results show a global sense in which a Lyapunov function is “strong”: it is not merely local, and it does not merely imply stability of one equilibrium. It can separate recurrent components, prescribe descent outside them, and, in the Conley framework, provide a gradient-like description of the entire dynamics.

4. Strong Lyapunov functions for nonsmooth, set-valued, and safety-critical dynamics

For differential inclusions of the form

$\mathcal{SCR}(\phi)$0

with $\mathcal{SCR}(\phi)$1 a Lipschitz Cusco multifunction and $\mathcal{SCR}(\phi)$2 a maximal monotone operator, the relevant object is a strong $\mathcal{SCR}(\phi)$3-Lyapunov pair $\mathcal{SCR}(\phi)$4. Here $\mathcal{SCR}(\phi)$5 are proper lower semicontinuous functions, $\mathcal{SCR}(\phi)$6, and

$\mathcal{SCR}(\phi)$7

must hold for every solution starting from every $\mathcal{SCR}(\phi)$8. The paper gives necessary and sufficient characterizations in terms of proximal or Fréchet subdifferentials, for example

$\mathcal{SCR}(\phi)$9

with $\mathcal{SCR}_d(\phi)$0 (Adly et al., 2018). Strength here lies in the universal quantifier: the inequality must hold for every admissible trajectory.

For convex processes $\mathcal{SCR}_d(\phi)$1, a strong Lyapunov function $\mathcal{SCR}_d(\phi)$2 is defined relative to the feasible set $\mathcal{SCR}_d(\phi)$3. It must be positive definite with respect to $\mathcal{SCR}_d(\phi)$4, and there must exist $\mathcal{SCR}_d(\phi)$5 such that

$\mathcal{SCR}_d(\phi)$6

This differs from earlier definitions that imposed positivity and decay on all of $\mathcal{SCR}_d(\phi)$7; the revised formulation reflects stability properties of nonstrict convex processes better because it is tied to feasible trajectories rather than arbitrary points of the domain (Eising et al., 2020). The same paper proves a weak-to-strong duality result: if $\mathcal{SCR}_d(\phi)$8 is a weak Lyapunov function for a closed convex process $\mathcal{SCR}_d(\phi)$9 and

$\dot V<0$0

then

$\dot V<0$1

is a strong Lyapunov function for the dual process $\dot V<0$2 (Eising et al., 2020).

In safety-critical control, the strengthened notion is a single control Lyapunov–barrier function $\dot V<0$3 for the control-affine system

$\dot V<0$4

The function must be continuously differentiable and positive definite, satisfy

$\dot V<0$5

and represent the safe set exactly: $\dot V<0$6 The paper proves that the existence of a strictly compatible pair of a control Lyapunov function and a control barrier function is equivalent to the existence of such a single smooth $\dot V<0$7, and gives a PDE characterization with prescribed boundary conditions on the safe set (Quartz et al., 15 Sep 2025). This is a particularly strong formulation because one scalar function simultaneously certifies asymptotic stability and safety, and because failure of such a unified certificate implies an unavoidable CLF–CBF conflict in the strict compatibility sense.

5. PDEs, optimization, and learning-based formulations

For hyperbolic systems of balance laws, strong Lyapunov functions appear as energy functionals with strict differential inequalities. In a $\dot V<0$8 linear hyperbolic system with source terms, classical exponential weights certify stability only under

$\dot V<0$9

Replacing exponential weights by the hyperbolic family

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$00

leads to Lyapunov functionals satisfying

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$01

under less conservative conditions, with optimized sufficient threshold

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$02

The spectral obstruction remains $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$03, but the new weights move Lyapunov-based certification substantially closer to that limit (Gugat, 2024). In the quasilinear setting, analogous hyperbolic-weight functionals yield

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$04

under smallness assumptions on source terms, wave-speed variation, and boundary gains (Gugat, 2024).

In first-order convex optimization, the strengthening is explicit in the definition. For the ODE

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$05

a strong Lyapunov condition is

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$06

This yields exponential decay when $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$07 and algebraic decay when $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$08, and the same template is applied to gradient flow, scaled gradient flow, heavy-ball flow, Nesterov-type flows, and Hessian-driven Nesterov accelerated gradient flow (Chen et al., 2021). The discrete counterpart uses inequalities such as

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$09

or

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$10

producing geometric or $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$11 convergence rates for proximal gradient, heavy-ball, and accelerated schemes (Chen et al., 2021).

A separate discrete-time construction based on Schur decomposition yields a particularly transparent strong Lyapunov function for first-order methods on strongly convex quadratics. Under the conjugate-pair spectral condition, the function

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$12

decreases monotonically for several methods, including Heavy Ball, Nesterov accelerated gradient, and NAG-GS. The construction is derived from the last coordinate in a Schur-transformed basis, where that coordinate is multiplied by a factor of modulus less than one at each step (Merkulov et al., 2023). This yields a strict discrete-time Lyapunov inequality tied directly to the spectral radius of the iteration matrix.

Neural approximations extend strong Lyapunov design to high-dimensional nonlinear systems. Lyapunov-Net defines

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$13

which guarantees positive definiteness by construction, and trains $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$14 by minimizing

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$15

Under the assumption that there exists $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$16 with

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$17

the paper proves approximation results showing that Lyapunov-Net can approximate such functions in $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$18, and that empirical risk minimization yields a $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$19-accurate Lyapunov function on $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$20 (Gaby et al., 2021). In this context, strength lies in simultaneous positive definiteness and a derivative bound of the form $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$21.

6. Computation, certification, and necessary conditions

Algorithmic construction of strong Lyapunov functions proceeds through several distinct verification paradigms. For polynomial ODEs, SOS programming enforces rate-coded inequalities and can maximize the certified rate $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$22 while producing invariant sublevel sets that serve as inner bounds on regions of performance (Jagt et al., 4 Jan 2026). For flows and maps near hyperbolic equilibria or fixed points, computer-assisted quadratic constructions use

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$23

together with strict negative definiteness of

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$24

for flows, or

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$25

for maps, over an explicit Lyapunov domain. These matrix inequalities are rigorously checked by interval arithmetic and are equivalent to cone-based procedures describing enclosures of stable and unstable manifolds (Matsue et al., 2016).

Formal verification replaces numerical sufficiency by symbolic certification. In SMT-based synthesis, a candidate $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$26 is checked against the negation of the Lyapunov conditions. For continuous-time systems, the verifier searches for

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$27

if the formula is unsatisfiable, then

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$28

hold throughout the modeled domain. The same framework extends to discrete-time systems and state-dependent switching systems, yielding formally correct strict Lyapunov or common Lyapunov certificates (Munser et al., 2021).

At the opposite end, one can rule out candidates before any dynamical verification. For a candidate $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$29 corresponding to a globally asymptotically stable equilibrium $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$30, define on each unit direction $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$31 the radial profile

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$32

and let $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$33 be the smallest $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$34 such that $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$35, or $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$36 if no such point exists. If $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$37 has a finite local minimum at $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$38, then

$$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$39

so $$-\nabla \mathcal L(x)\cdot \mathcal G(x)\ge c(x)\,\mathcal L(x)^q + p(x)^2,$$40 cannot be a Lyapunov function for a globally asymptotically stable equilibrium. The generalized steepest descent method proposed in that setting is therefore a first-level rejection test independent of the system dynamics (Athalye et al., 2014).

The main limitations are domain-specific. SOS formulations require polynomial vector fields or suitable transformations (Jagt et al., 4 Jan 2026). Complete Lyapunov prescription theorems prescribe the orbital derivative only on compact subsets of the non-chain-recurrent region (Giesl et al., 2021). Hyperbolic-weight PDE constructions improve stabilizability bounds but remain sensitive to source size, delay, and boundary gains (Gugat, 2024). Schur-based optimization Lyapunov functions rely on quadratic objectives and specific spectral structure (Merkulov et al., 2023). This suggests that strong Lyapunov theory is best understood as a collection of sharpened Lyapunov paradigms, each matched to a particular analytical or computational regime, rather than as a single universal doctrine.