Discrete Lyapunov Equation

Concept Overview

For the discrete-time linear dynamical system xk+1 = A xk, a quadratic candidate Lyapunov function V(x) = xT P x changes over one time step by V(Ax) - V(x) = xT(ATPA - P)x. Requiring this one-step energy change to equal -||x||2 = -xTIx for all states x leads directly to the discrete Stein/Lyapunov equation ATPA - P = -I.

By Stein's theorem, the system is asymptotically stable (spectral radius ρ(A) < 1) if and only if there exists a unique symmetric positive-definite matrix P satisfying ATPA - P = -I.

Mathematical Formulation

For an upper-triangular 2×2 matrix A = [[a, b], [0, d]], the linear operator P ↦ ATPA - P has eigenvalues a2 - 1, ad - 1, and d2 - 1 on the space of symmetric 2×2 matrices. Consequently, the equation has a unique symmetric solution if and only if (1 - a2)(1 - ad)(1 - d2) ≠ 0.

When a unique solution P exists, it satisfies V(Ax) - V(x) = -||x||2 for every probe vector x, even if A is unstable! What distinguishes stability from mere algebraic solvability is whether P is positive definite (p11 > 0 and p11p22 - p122 > 0).

Worked Calculation

Let us consider a 2x2 upper triangular matrix A where a = 0.5, b = 1, and d = 0.5. The discrete Lyapunov equation to solve is AT P A - P = -I.

Since A is upper triangular, we can parameterize the symmetric matrix P as having entries p11, p12, and p22.

1. First, from the (1,1) entry of the matrix equation, we have a2 p11 - p11 = -1. Solving for p11 yields p11 = 1 / (1 - a2). For a = 0.5, this gives p11 = 1 / 0.75 = 4/3 ≈ 1.3333333333333333.
2. Next, using the (1,2) entry, we find a b p11 + a d p12 - p12 = 0. Solving for p12 gives p12 = (a b p11) / (1 - a d). Substituting our values gives p12 = (0.5 × 1 × 4/3) / (1 - 0.25) = (2/3) / (0.75) = 8/9 ≈ 0.8888888888888888.
3. Finally, the (2,2) entry gives b2 p11 + 2 b d p12 + d2 p22 - p22 = -1. Solving for p22 yields p22 = (1 + b2 p11 + 2 b d p12) / (1 - d2). Substituting the values gives p22 = (1 + 1 × 4/3 + 2 × 0.5 × 8/9) / (1 - 0.25) = (1 + 4/3 + 8/9) / 0.75 = (29/9) / 0.75 = 116/27 ≈ 4.296296296296297.

Experiment

We can explore how the parameters affect the stability and the existence of a unique positive-definite solution.

1. Open the Discrete Lyapunov Equation module.
2. Set the parameter a to 0.5, b to 1, and d to 0.5.
3. Observe that the solution P is calculated and the system is classified as "Stable", with "Positive definite P" being true.
4. Change the parameter a to 2 and d to 2 with b = 0.
5. Notice that the asymptotic stability changes to "Not asymptotically stable", and "Positive definite P" becomes false, indicating the system is unstable and the quadratic energy form is no longer positive definite.
6. Try setting a to 1 and d to 0. The spectral radius becomes 1, and you will see that there is "No unique solution" because the matrix operator is singular.

Limitations

This module is specifically restricted to 2x2 upper-triangular matrices with discrete parameter choices. Additionally, this implementation uses closed-form algebraic elimination in double-precision floating-point arithmetic. When the operator is singular ((1 - a2)(1 - ad)(1 - d2) = 0), the module reports "No unique solution" and hides P-dependent values rather than claiming that every singular case is inconsistent.

Common Pitfalls

A common misconception is that finding a symmetric matrix P with V(Ax) - V(x) = -||x||2 < 0 automatically proves stability. As the case A = 2I with P = -(1/3)I shows, V(x) decreases only because it is negative definite and growing more negative as ||xk|| → ∞; positive definiteness of P is essential for a valid Lyapunov certificate.

References

The following reference provides the implementation details for solving the discrete Lyapunov equation computationally in standard scientific libraries. Note the transpose convention used in different mathematical software ecosystems.

• SciPy solve_discrete_lyapunov Documentation. SciPy uses A X AH - X + Q = 0. The proposed AT P A - P + I = 0 energy form corresponds to its argument a = A.T.