Eigenvalues: −1 ± i (complex, Re < 0)
Trajectories spiral inward toward the origin. The rate of convergence is set by Re(λ) = −1; the frequency of rotation is Im(λ) = 1.
Lyapunov: Any positive definite quadratic V = xᵀPx works. The isotropic V = x²+y² gives dV/dt = −2(x²+y²) < 0 everywhere: a strict global Lyapunov function.
Key insight: The symmetric part of A is ½(A+Aᵀ) = −I, which is negative definite — a sufficient condition for V=xᵀx to be a Lyapunov function with no calculation required.
Eigenvalues: −2, −1 (real, distinct, negative)
Trajectories approach the origin tangent to the slow eigenvector (y-axis, eigenvalue −1). Unlike a focus, there is no rotation.
Lyapunov: V = x²+y² works but gives non-circular level curves relative to the flow. The "natural" candidate from solving AᵀP+PA = −I is an ellipse aligned with the eigenvectors.
Key insight: The general quadratic with a nonzero cross-term b can give a tighter ROA estimate than the isotropic one because it better matches the geometry of the flow.
Eigenvalues: 1 ± i (complex, Re > 0)
Trajectories spiral outward. The time-reverse of the stable spiral.
Lyapunov: No Lyapunov function can certify stability. dV/dt will be positive somewhere near the origin for any positive definite V.
Eigenvalues: +1, −2 (real, opposite sign)
The y-axis is the stable manifold; the x-axis is the unstable manifold. Trajectories on the stable manifold converge to the origin; all others eventually escape.
Lyapunov: Not certifiable. Even though some trajectories approach the origin, the unstable manifold always provides an escape route.
Jacobian at origin: zero matrix — linearization predicts nothing. Stability must be established by nonlinear Lyapunov analysis.
Lyapunov: V = x²+y² gives dV/dt = −2(x⁴+y⁴) < 0 everywhere except the origin. The quartic candidate V = x⁴+y⁴ also works and gives dV/dt = −4(x⁶+y⁶). Convergence is slower than exponential (polynomial rate).
Key insight: This system illustrates why Lyapunov theory extends beyond linearization — it is the primary tool when linearization gives no information.
Eigenvalues: ±i (purely imaginary — center)
A conservative system. Trajectories are perfect circles; V = x²+y² = r² is constant along every orbit. The origin is Lyapunov stable (trajectories stay close) but not asymptotically stable (they never converge).
Connection to Hamiltonian mechanics: H = ½(x²+y²) = ½p² + ½q² is the total energy (T = ½p², U = ½q²). dH/dt = 0 identically — energy is exactly conserved.
LaSalle: Cannot strengthen the conclusion here because {dV/dt = 0} = entire phase plane. The largest invariant set is everything, not just the origin.
Eigenvalues of linearization at origin: μ/2 ± i√(1−μ²/4) — stable focus for μ < 0, unstable focus for μ > 0, center for μ = 0.
The damping term μ(1−x²)y is negative (dissipative) for |x| < 1 and positive (pumping energy in) for |x| > 1. This nonlinear balance creates a stable limit cycle for μ > 0.
μ < 0: Stable focus. This is the reverse-time Van der Pol oscillator, the standard region-of-attraction benchmark. The isotropic V = x²+y² gives dV/dt = 2μ(1−x²)y², negative inside the strip |x| < 1, and so it certifies the unit disk (c* = 1). That is only the estimate from this particular V: the true ROA is the interior of the unstable limit cycle (near amplitude 2), and optimizing the general quadratic certifies a noticeably larger region.
μ > 0: Unstable focus; energy pumped into small oscillations. Limit cycle at r ≈ 2. No Lyapunov certificate for the origin.
Historical note: Originally derived by Balthazar van der Pol (1926) to model oscillations in vacuum tube circuits.
In polar: ṙ = r(μ − r²), θ̇ = 1. Eigenvalues of linearization at origin: μ ± i — stable focus for μ < 0, unstable for μ > 0.
The canonical form for a supercritical Hopf bifurcation. As μ increases through zero, the stable focus loses stability and a stable limit cycle is born at r = √μ.
Lyapunov: V = x²+y² = r² gives dV/dt = 2r²(μ−r²). For μ < 0: dV/dt < 0 for all r > 0, and so the origin is globally asymptotically stable and V certifies it with no conservatism. For μ > 0: dV/dt > 0 for r < √μ, and so the origin is unstable (trajectories spiral outward to the stable limit cycle at r = √μ) and no Lyapunov function certifies it.
Key insight: the μ < 0 case is one of the rare instances where a plain quadratic certifies global stability exactly — dV/dt is strictly negative everywhere off the origin, and so the certified region is the entire plane with no gap to the true basin.
Equilibria: stable at x = 2kπ (eigenvalues (−c ± √(c²−4))/2), saddles at x = (2k+1)π. c < 2: stable focus; c = 2: critically damped node; c > 2: overdamped node.
A mechanical system (unit mass, unit length pendulum) with viscous damping c. State: x = angle, y = angular velocity.
Hamiltonian: H = (1−cos x) + ½y² = potential energy + kinetic energy. With damping: dH/dt = −cy² ≤ 0. The cross-terms cancel exactly due to the antisymmetry of Hamilton's equations — energy dissipation is the only survivor.
Lyapunov: The candidate is the energy itself, V = H = (1−cos x) + ½y², proportional to the Hamiltonian. It is positive definite about the hanging equilibrium, and dV/dt = −cy² ≤ 0, which makes V a Lyapunov function whenever c > 0. The decrease is nonstrict (it vanishes on the whole line y = 0), and so asymptotic stability follows from LaSalle rather than from strict decrease.
ROA via LaSalle: with the energy function, c* = H(±π, 0) = 2, certifying the sublevel set {H < 2} bounded by the undamped separatrix. The true basin (bounded by the stable manifolds of the saddles at ±π) is somewhat larger, but this energy estimate captures essentially all of it.
Potential: U(x) = −½x² + ¼x⁴, with wells at x = ±1 and a local maximum at x = 0.
Eigenvalues of linearization at origin: (−c ± √(c²+4))/2 — always real with opposite signs, confirming a saddle for all c ≥ 0.
Eigenvalues of linearization at (±1, 0): (−c ± √(c²−8))/2 — stable focus for c < 2√2 ≈ 2.83, stable node for c > 2√2.
A damped particle in a double-well potential. Used in structural mechanics to model snap-through buckling.
Hamiltonian: H = ¼(1−x²)² + ½y². Level curves of H are closed around each well for H < ¼ (inside the homoclinic orbit through the saddle at origin).
Lyapunov: About either well the candidate is again the energy, V = H = ¼(1−x²)² + ½y², proportional to the Hamiltonian and normalized so that it vanishes at the well. It is positive definite there, and dV/dt = −cy² ≤ 0, and so LaSalle gives asymptotic stability of that well for c > 0. The estimate is local: the relevant sublevel set is the part of {H < ¼} surrounding one well, bounded by the homoclinic orbit through the saddle at the origin.
Multiple basins: initial conditions may converge to different wells depending on which side of the stable/unstable manifold of the origin they start on.
Lotka–Volterra prey–predator dynamics, translated so the coexistence equilibrium is at the origin (x = prey deviation, y = predator deviation). Eigenvalues of linearization at origin: ±i — a center in the linear approximation.
Conserved quantity: V = (x − ln(1+x)) + (y − ln(1+y)) satisfies dV/dt = 0 identically — it is a first integral, playing the same role as the Hamiltonian in a mechanical system. Level curves of V are the closed orbits (population cycles).
Interpretation: This is a non-mechanical example of a Hamiltonian-like conserved quantity arising from ecology. The "energy" here is an information-theoretic quantity (related to Kullback-Leibler divergence from equilibrium populations).
Adds proportional mortality c (constant per-capita death rate) to both species. Eigenvalues of linearization at origin: −c ± i — stable focus for all c > 0.
Lyapunov: The same V = (x−ln(1+x)) + (y−ln(1+y)) now gives dV/dt = −c[x²/(1+x) + y²/(1+y)] ≤ 0. This mirrors the mechanical case: adding dissipation to a conservative system turns a first integral into a Lyapunov function.
LaSalle: {dV/dt = 0} only at origin → origin is asymptotically stable. The certified region via this V is bounded by where V becomes undefined (x = −1 or y = −1, i.e. population extinction).
This single family contains every quadratic Lyapunov candidate, written as V = xᵀPx with P = [[1, b/2], [b/2, c]]. The naming follows the usual quadratic convention ax² + bxy + cy² with the leading coefficient normalized to a = 1. It is positive definite exactly when c > 0 and b² < 4c (the determinant condition det P = c − b²/4 > 0); outside that region V is not a valid candidate and the explorer hatches the plane as "V undefined".
Why fixing the x² coefficient to 1 loses nothing: scaling V by a positive constant leaves the certified ROA identical (if V₂ = kV₁ then dV₂/dt = k·dV₁/dt has the same sign everywhere, and {V₂ ≤ kc} = {V₁ ≤ c} as sets). So every p·x² + q·xy + r·y² normalizes to this form, and the diagonal special cases are just b = 0: c = 1 is the isotropic V = x²+y² (circles), c > 1 weights y, c < 1 weights x. The cross-term b rotates the ellipse off the coordinate axes.
Connection to the Lyapunov matrix equation: for a linear system ẋ = Ax, solving AᵀP + PA = −Q for positive definite Q yields a P that is generally not diagonal — the cross-term b is exactly that off-diagonal coupling. A tilted or stretched ellipse routinely certifies a larger ROA than the circle for systems whose flow is not aligned with the axes (the stable spiral, the damped pendulum).
Presets in the Candidate Parameters panel jump to the named special cases (isotropic, y-weighted, x-weighted, cross-term); the sliders then move freely from there.
Higher-degree candidate. Necessary for systems like the Cubic Drain where the Jacobian at the origin is zero and quadratic candidates may give uninformative dV/dt = 0 at the origin.
Trade-off: The gradient vanishes to second order at the origin, making the level curves "flat" near zero. This means dV/dt is very small near the origin even when negative — the color map may appear nearly uniform there.
The total mechanical energy of the pendulum: potential U = 1−cos x (= mgl(1−cosθ) in physical units) plus kinetic T = ½y² (= ½mℓ²ω²). The ½ coefficient on y² is not arbitrary — it is the unique coefficient that produces exact cancellation of cross-terms.
ROA boundary: V(±π, 0) = 2. The level curve V = 2 is the separatrix — the boundary between the basin of attraction of the origin and the "over-the-top" trajectories that converge to x = ±2π.
Potential: U = ¼(1−x²)² — double well with minima at x = ±1.
Total energy for the Duffing oscillator. The potential U has minima at x = ±1 (the stable equilibria) and a local maximum at x = 0 (the saddle). V = 0 at (±1, 0), not the origin — so V is not positive definite at the origin, but is centered on each well.
Note: Pairing this candidate with the Duffing system correctly certifies the ROA around each well. Pairing it with the origin-centered view reveals why no certificate exists there.
The natural Lyapunov function for the centered Lotka–Volterra system. Note that t − ln(1+t) ≥ 0 for all t > −1 with equality only at t = 0 — so V is positive definite near the origin.
Undamped: dV/dt = 0 everywhere — V is a first integral (conserved quantity). Level curves = population cycles.
With damping c: dV/dt = −c[x²/(1+x) + y²/(1+y)] ≤ 0. Exactly mirrors the mechanical case: damping extracts the ecologically conserved quantity, driving populations to equilibrium.
Information-theoretic interpretation: V is related to the Kullback-Leibler divergence between the current population state and the equilibrium — a measure of "surprise" or distance from equilibrium in an information sense.
Certifying that a polynomial is non-negative everywhere is NP-hard in general, but certifying that it is a sum of squares — a sufficient condition for non-negativity — is convex. A polynomial p(x) is SOS exactly when p(x) = z(x)ᵀ Q z(x) for a positive semidefinite matrix Q, with z(x) a vector of monomials; searching for such a Q is a semidefinite program (SDP), solvable to global optimality in polynomial time. Lyapunov's conditions — V positive definite, dV/dt non-increasing — are exactly the polynomial non-negativity statements this turns into an SDP.
Parameterize V as a polynomial with unknown coefficients and require V and −dV/dt each to dominate a fixed positive definite function (each of V − ℓ₁ and −dV/dt − ℓ₂ is SOS). When feasible, the program returns a V whose existence certifies asymptotic stability, searching a whole family of candidates at once rather than testing them one by one. The one choice that matters in practice is how strictly to enforce the decrease: a strict margin ℓ₂ certifies convergence outright but is infeasible for degenerate cases such as the Cubic Drain, and so this tool drops it and finishes convergence with a short per-system argument. That trade-off, and the one-line change that restores the strict form, is developed on the SOS Synthesis tab.
The "Optimize parameters for largest ROA" button on the Explorer tab is a low-dimensional stand-in for the same idea: it searches the two free entries of the quadratic's matrix P by brute-force grid evaluation, where SOS searches the whole cone of positive definite matrices by convex optimization. The crucial difference is what comes out — the button's region is a numerical estimate from sampling the grid, whereas SOS produces an actual certificate that holds at every point. Enlarging the certified region, and scaling beyond the quadratic to higher-degree V or non-polynomial dynamics, is where the dedicated software on the SOS Code tab takes over.
An equilibrium x* is Lyapunov stable if trajectories that start close to x* stay close for all future time. Formally: for every ε > 0 there exists δ > 0 such that ‖x(0)−x*‖ < δ ⟹ ‖x(t)−x*‖ < ε for all t ≥ 0.
This is a boundedness condition — it rules out trajectories that drift away, but does not require convergence. A center (e.g. the SHO ★) is Lyapunov stable but not asymptotically stable.
An equilibrium x* is asymptotically stable if it is Lyapunov stable and trajectories starting sufficiently close to x* converge to x* as t → ∞. The set of initial conditions from which convergence occurs is the region of attraction (ROA).
Asymptotic stability is strictly stronger than Lyapunov stability. All the stable spirals, nodes, and foci in the Explorer ★ are asymptotically stable.
An equilibrium is globally asymptotically stable if it is asymptotically stable and its region of attraction is the entire state space. The Cubic Drain ★ and Stable Spiral ★ are GAS; the Damped Pendulum ★ is only locally asymptotically stable (the ROA is bounded by the separatrix).
For a GAS equilibrium, a strict Lyapunov function with dV/dt < 0 everywhere certifies global convergence directly, without needing to compute the ROA boundary.
An equilibrium x* is exponentially stable if ‖x(t)−x*‖ ≤ ke−αt‖x(0)−x*‖ for positive constants k, α. This is the strongest form: it specifies an exponential rate of convergence.
Asymptotic stability does not imply exponential stability in general — the Cubic Drain ★ is GAS but converges at a polynomial rate (not exponential) because its Jacobian at the origin is zero.
An equilibrium that is not Lyapunov stable. Trajectories starting arbitrarily close can escape to any distance. Unstable foci ★, saddle points ★, and the origin of the Duffing system ★ are all unstable.
Note the distinction between a saddle (some trajectories approach, some escape — along the stable/unstable manifolds) and an unstable focus (all nearby trajectories spiral away).
Stability is local if it is guaranteed only for initial conditions within some neighborhood of the equilibrium; global if it holds from the entire state space. A Lyapunov function certifies local stability when its sublevel set is bounded; global stability when dV/dt < 0 everywhere with V → ∞ as ‖x‖ → ∞ (a radially unbounded V).
A point x* where f(x*) = 0 — the system is stationary there. Linearization around x* via the Jacobian J = ∂f/∂x|x* classifies the local behavior by the eigenvalues of J.
Classification of a 2D equilibrium by the eigenvalues λ₁, λ₂ of its Jacobian:
Node: both eigenvalues real, same sign. Stable node (Re < 0 ★), unstable node (Re > 0). Trajectories approach/leave tangent to the slow eigenvector.
Focus (Spiral): complex conjugate eigenvalues (λ = α ± βi). Stable focus (α < 0 ★), unstable focus (α > 0 ★). Trajectories spiral in or out.
Center: purely imaginary eigenvalues (λ = ±βi). Closed orbits, Lyapunov stable only. SHO ★, undamped predator–prey ★.
Saddle: real eigenvalues of opposite sign ★. Has a stable manifold (trajectories converging to it) and an unstable manifold (trajectories leaving). Always unstable.
An isolated closed orbit — one that nearby trajectories spiral toward (stable limit cycle) or away from (unstable). Unlike the closed orbits of a center, a limit cycle is an isolated trajectory that attracts (or repels) nearby ones.
The Van der Pol oscillator ★ (μ > 0) and Hopf Normal Form ★ (μ > 0) both have stable limit cycles. No quadratic Lyapunov function can certify stability of the limit cycle itself, though one can certify the ROA around the unstable origin.
A bifurcation in which a stable equilibrium loses stability as a parameter crosses a critical value, with a limit cycle being born simultaneously. In a supercritical Hopf bifurcation the born limit cycle is stable; in a subcritical one it is unstable.
The Hopf Normal Form ★ is the canonical example: at μ = 0 the origin changes from stable focus to unstable focus, and a stable limit cycle grows from it for μ > 0 with radius r = √μ.
A set S is positively invariant if every trajectory that starts in S stays in S for all t ≥ 0. It is invariant (two-sided) if this holds for t ∈ ℝ. A sublevel set {V ≤ c} where dV/dt ≤ 0 on its boundary is positively invariant — this is the key mechanism behind LaSalle's principle ★.
For a saddle point x*, the stable manifold Ws is the set of all initial conditions whose trajectories converge to x* as t → +∞; the unstable manifold Wu is the set converging to x* as t → −∞ (i.e. diverging forward in time).
For a linear saddle, Ws and Wu are the eigenspaces of the negative and positive eigenvalues respectively. Nonlinearly, they are smooth curves tangent to those eigenspaces at x*.
The stable manifold of the saddle at (±π, 0) in the Damped Pendulum ★ forms the separatrix bounding the ROA of the origin.
A curve (or surface) in phase space that divides regions of qualitatively different long-term behavior. In 2D, a separatrix is typically the stable manifold of a saddle point.
In the Damped Pendulum ★, the separatrix separates initial conditions that converge to the origin (pendulum coming to rest upright) from those that wrap around to adjacent equilibria (pendulum going over the top). The Lyapunov analysis correctly identifies V = 2 as the energy level of this separatrix.
A homoclinic orbit is a trajectory that connects a saddle point to itself — it lies in both the stable and unstable manifold of the same equilibrium. It forms a loop enclosing one basin of attraction.
A heteroclinic orbit (or heteroclinic connection) connects two different saddle points — it lies in the unstable manifold of one and the stable manifold of the other.
The Duffing Double-Well ★ has a pair of homoclinic orbits through the origin (the figure-eight shaped separatrix visible at c = 0). The undamped pendulum ★ has homoclinic orbits through each saddle at x = ±π.
With damping (c > 0), these become the separatrices bounding the ROA of each well: trajectories starting inside converge to the well, those outside escape to another equilibrium.
A function V : ℝⁿ → ℝ is positive definite (PD) on a neighborhood of the origin if V(0) = 0 and V(x) > 0 for all x ≠ 0 in that neighborhood. It is positive semi-definite (PSD) if V(x) ≥ 0 (allowing V to be zero away from the origin).
For a Lyapunov function, PD is required for V. The condition on dV/dt allows PSD (dV/dt ≤ 0) — this is the setting for LaSalle's principle. Strict negativity (dV/dt < 0 except at origin) directly gives asymptotic stability without LaSalle.
A continuously differentiable function V : ℝⁿ → ℝ that is positive definite and has dV/dt = ∇V · f(x) ≤ 0 along trajectories. If such a V exists, x* = 0 is at minimum Lyapunov stable.
Lyapunov's theorem: if V is PD and dV/dt is negative definite (ND), then x* = 0 is asymptotically stable. If additionally V is radially unbounded (V → ∞ as ‖x‖ → ∞), then x* = 0 is GAS.
There is no systematic method for finding V — construction remains an art, though sum-of-squares (SOS) programming and neural-network methods automate the search for certain problem classes.
A function V with dV/dt = 0 identically along all trajectories. Level curves of V are exactly the system's orbits. A first integral proves Lyapunov stability but cannot prove asymptotic stability — the largest invariant set in {dV/dt = 0} is the entire phase space.
The Hamiltonian H of an undamped mechanical system is a first integral ★ (SHO, undamped pendulum, undamped predator–prey). Adding dissipation breaks the conservation: dV/dt = −D ≤ 0, turning a first integral into a Lyapunov function.
The set of all initial conditions x(0) from which the trajectory converges to the equilibrium: ROA = {x : x(t) → x* as t → ∞}. For a GAS equilibrium the ROA is all of ℝⁿ; for a locally stable one it is some proper subset.
Lyapunov analysis gives certified inner approximations of the ROA: any sublevel set {V ≤ c} that lies entirely in {dV/dt < 0} is positively invariant and contained in the true ROA. The orange dashed boundary in the Explorer ★ marks the largest such set certifiable with the chosen V.
Extends asymptotic stability conclusions to cases where dV/dt ≤ 0 (semi-definite) rather than dV/dt < 0 (definite).
Statement: Let Ω be a compact positively invariant set. Let E = {x ∈ Ω : dV/dt(x) = 0}. Let M be the largest invariant set contained in E. Then every trajectory starting in Ω converges to M as t → ∞.
How it's applied: if M = {0} (the only invariant trajectory in {dV/dt = 0} is the equilibrium itself), then the origin is asymptotically stable for all initial conditions in Ω. This handles, e.g., the Damped Pendulum ★ where dV/dt = −cy² is zero on the entire y = 0 axis, but no trajectory other than the equilibrium can stay on that axis.
LaSalle's principle strictly generalizes Lyapunov's theorem: whenever dV/dt is ND, the set E = {0} trivially and LaSalle recovers the same conclusion.
For a linear system ẋ = Ax, the quadratic V = xᵀPx (P symmetric positive definite) has dV/dt = xᵀ(AᵀP + PA)x. For dV/dt to be ND we need AᵀP + PA = −Q for some PD matrix Q — this is the Lyapunov matrix equation.
The equation always has a unique PD solution P for any PD Q if and only if A is Hurwitz (all eigenvalues have negative real part). Solving it for Q = I gives the "natural" Lyapunov function for the linear system. For nonlinear systems, this solution at the linearization gives a local Lyapunov function valid in a neighborhood of the origin.
The general quadratic candidate V = x² + b·xy + c·y² in the Explorer ★ represents exactly this structure: it is the general 2×2 symmetric positive definite quadratic xᵀPx with P = [[1, b/2], [b/2, c]], the leading coefficient normalized to 1. Optimizing the two sliders is searching for the P that best certifies the region of attraction.
A dynamical system is Hamiltonian if it can be written as ẋ = ∂H/∂p, ṗ = −∂H/∂x for a scalar function H(x, p) called the Hamiltonian (usually total energy: H = T + U = kinetic + potential).
The defining property is exact conservation: dH/dt = (∂H/∂x)ẋ + (∂H/∂p)ṗ = (∂H/∂x)(∂H/∂p) + (∂H/∂p)(−∂H/∂x) = 0. The two terms cancel identically — a consequence of the antisymmetric (symplectic) structure of Hamilton's equations.
This means H is always a first integral ★. Adding dissipation breaks the cancellation and leaves dH/dt = −D ≤ 0 (where D ≥ 0 is the dissipation), turning the Hamiltonian into a Lyapunov function for the damped system.
Examples in the Explorer ★: SHO (H = ½x² + ½y²), undamped pendulum (H = (1−cos x) + ½y²), undamped predator–prey (H = (x−ln(1+x)) + (y−ln(1+y))), Duffing (H = ¼(1−x²)² + ½y²). Each becomes a Lyapunov function when its corresponding damped version is selected.
A system with input u and output y is passive if there exists a non-negative storage function V(x) ≥ 0 such that the energy balance inequality holds: dV/dt ≤ uᵀy. In words: the rate of change of stored energy is at most the externally supplied power uᵀy — the system cannot generate energy internally.
When u = 0 (autonomous), this reduces to dV/dt ≤ 0 — the storage function V is a Lyapunov function. All the Hamiltonian-based Lyapunov functions in the Explorer ★ are storage functions for passive systems.
Passivity and Lyapunov stability are deeply linked: a passive system with a positive definite storage function and zero-state observability (the output y = 0 implies x → 0) is asymptotically stable. This is essentially LaSalle's principle in input-output language.
Key benefit: passivity is preserved under negative feedback interconnection — connecting two passive systems in feedback gives a passive composite system. This makes passivity a powerful framework for stability analysis of complex networked systems (power grids, multi-robot systems, biomechanical models).
Passivity-based control designs a controller by shaping the system's energy function rather than canceling its natural dynamics. The two key steps are:
1. Energy shaping: find a control law that makes the closed-loop system behave as if its Hamiltonian were a modified Hd(x) with a minimum at the desired equilibrium x*. This modifies the potential energy landscape so the system "wants" to rest at x*.
2. Damping injection: add a dissipative term so that dHd/dt ≤ 0 with equality only at x*. This ensures convergence via LaSalle's principle.
Concrete example (pendulum ★): The uncontrolled pendulum with damping c already satisfies dH/dt = −cy² ≤ 0, with H = (1−cos x) + ½y². LaSalle certifies that the only invariant set in {y = 0} is the origin x = 0 (the hanging equilibrium). To instead stabilize the inverted position (x = π), one would reshape Hd to have a minimum there — this is the "swing-up and balance" problem in robotics.
Why PBC instead of feedback linearization? Feedback linearization cancels the system's nonlinearities exactly, requiring precise knowledge of the model and often demanding large control effort. PBC works with the system's natural dynamics; it is inherently robust and often requires less energy.
The textbook SOS encoding of asymptotic stability mirrors Lyapunov's theorem term for term. Parameterize V as a polynomial with unknown coefficients and impose two SOS constraints, each holding a quantity above a fixed positive definite function ℓ:
V − ℓ₁(x) is SOS (V is positive definite)
−dV/dt − ℓ₂(x) is SOS (dV/dt is negative definite)
with ℓ₁, ℓ₂ positive definite, most often ℓ = ε‖x‖². The second constraint forces dV/dt ≤ −ℓ₂(x) < 0 for every x ≠ 0, which is strict decrease, and so asymptotic stability — in fact exponential stability, for a quadratic ℓ₂ — follows directly, with no further argument. This is the default in the SOS-synthesis literature: Parrilo (2000), Papachristodoulou and Prajna (2002), and the region-of-attraction work of Jarvis-Wloszek et al. (2003), Tan and Packard (2008), and Topcu, Packard, and Seiler (2008) all carry a positive definite margin on the derivative, as does the Lyapunov example in the SOSTOOLS user's guide. The attraction of this form is a self-contained certificate: the solver proves convergence, not merely non-increase.
A quadratic margin ℓ₂ = ε‖x‖² demands that V decrease at least quadratically near the equilibrium. That holds whenever the linearization is non-degenerate, a hyperbolic equilibrium, which covers most systems. It fails at a non-hyperbolic equilibrium whose decrease is inherently higher order. The Cubic Drain (ẋ = −x³) is the clean example: with V = ‖x‖², its derivative is dV/dt = −2(x₁⁴ + x₂⁴), strictly negative for x ≠ 0 yet bounded below by no multiple of ‖x‖². Requiring −dV/dt − ε‖x‖² to be SOS is then infeasible, even though the origin is perfectly asymptotically stable.
The conventional fix is not to abandon the margin but to match its degree to the decay order: use ℓ₂ ∝ ‖x‖⁴ for quartic decay, and a positive definite form of the right degree in general. That keeps the certificate self-contained, at the cost of choosing the margin degree system by system — which a single uniform template cannot do.
To run one template across the whole catalog, the Cubic Drain included, without per-system margin tuning, the generated code drops ℓ₂ from the derivative constraint and asks only for
V − ε‖x‖² is SOS, λ is SOS, −dV/dt − λ·(R² − ‖x‖²) is SOS.
V stays positive definite — a quadratic floor is fine on the V side because the degenerate decay lives in dV/dt, not in V — while the derivative is held only to dV/dt ≤ 0 on the disk ‖x‖ ≤ R. The multiplier λ is the S-procedure: instead of demanding −dV/dt ≥ 0 everywhere, the code demands −dV/dt − λ·(R² − ‖x‖²) ≥ 0 with λ itself SOS, and so the bound is forced only where R² − ‖x‖² ≥ 0 — that is, inside the disk — and may be relaxed outside it. A non-strict bound certifies Lyapunov stability and forward-invariance of the sublevel set, but not convergence on its own. Convergence is recovered by one short, explicit argument per system, and supplying that argument is the deliberate teaching point rather than a gap. It also matches how dissipative and energy methods are taught: for those systems dV/dt ≤ 0 is the natural statement and an invariance argument does the rest.
The trade is worth naming plainly. The conventional form moves the entire proof into the SDP; this form keeps the SDP uniform and leaves a one-line finishing step. Neither is more correct — they are two standard camps, and the dissipative systems below genuinely belong to the second.
The default template on the SOS Code tab proves only dV/dt ≤ 0 for the polynomial V it finds, which gives stability and an invariant region but not convergence. Two routes supply the rest, and which one applies is a property of the Lyapunov function, not of the system: every system here already has a region of attraction worked out by energy or LaSalle arguments on the Explorer and Reference tabs, with no SOS at all.
Tighten the SDP. Re-add a strict margin (section 1). Every equilibrium in the catalog except the Cubic Drain is hyperbolic, and a quadratic margin then certifies convergence outright. The Cubic Drain is non-hyperbolic, and its quartic decay needs a quartic margin to match (section 2); the SDP still certifies it. The non-strict default exists to dodge this per-system tuning, not because the SDP cannot certify these systems.
Argue by hand. Leave the certificate non-strict and establish convergence separately. For a V whose dV/dt is already negative away from the origin, that is a one-line observation. For a V whose dV/dt vanishes on more than the origin — the textbook case being the energy function on the Damped Pendulum and the Duffing wells, where dV/dt = −c·y² is zero along y = 0 — LaSalle's invariance principle applies, and the largest invariant set inside {dV/dt = 0} reduces to the origin. That LaSalle step belongs to the energy function, not to the pendulum: handed to the SDP with a margin, the same system admits a different V that needs no such step.
The SOS program certifies the largest level set of V that fits inside a disk of radius R, and the disk caps how big the certified region can be. Two routes push past a single hand-picked disk, and they trade off differently.
Alternating (V–s) iteration. The textbook route searches over the region itself, enlarging the certified set and re-fitting V to it in alternation: fix V and grow the multipliers and the level, then fix those and re-solve for V, and repeat. Each half is a convex SDP, but the product of V and the multiplier makes the joint problem bilinear and non-convex, and the iteration can stall in a local optimum and depends on where it starts. This is the standard approach of Tan and Packard (2008) and Topcu, Packard, and Seiler (2008).
Bisection on the region radius. A simpler alternative keeps every step convex and searches a single scalar instead. Hold the formulation fixed and bisect on R: solve the convex SDP at a trial radius, grow R while it stays feasible, and shrink it when it does not, converging to the largest disk the certificate can cover. Each evaluation re-finds V from scratch – the Lyapunov function is never iterated against the multiplier – and feasibility is monotone in R, which makes the bisection well posed because a disk that is certifiable stays certifiable when shrunk. Driving R toward the true basin boundary this way forces the SDP's V close to what converse Lyapunov theory guarantees exists: Massera's theorem (1956) says an asymptotically stable system admits a smooth Lyapunov function whose sublevel sets fill the basin, and the largest-disk certificate is the degree-limited SOS approximation of it. The construction is not literally unique because many converse integrands work, and V here is capped at degree d, but in practice the bisection recovers very nearly the full basin with far less tuning than the alternating iteration. It is the same move as the Explorer tab's "Optimize parameters for largest ROA" button, which searches the quadratic's two free parameters for the largest certified area; here, the scalar searched is the region radius rather than the shape of V.
The generated code takes the bisection route by default: its auto_R flag wraps the single
convex program in exactly this search, and turning it off certifies the single fixed disk R instead. Enlarging the
region beyond the largest disk, reshaping it to follow the flow rather than a circle, returns to the
non-convex territory of the alternating method.
SOS requires the system to have a polynomial vector field. For systems with non-polynomial terms — like sin(x₁) in the Damped Pendulum — those terms can be replaced by their degree-k Taylor polynomials before applying the SOS method. However, the certificate is then rigorous for that approximation, not the true system, and the two agree only near the origin. On its own the SDP therefore proposes a candidate, not a proof about the real dynamics — and for a hyperbolic equilibrium that adds little to the linearization.
Two routes restore rigor: bound the Taylor remainder and carry a strict margin large enough to absorb it, which keeps the decrease valid once the true f is substituted back (Han and Panagou, 2017); or recast the nonlinearity exactly, as in using trigonometric identities (e.g., s² + c² = 1 when s = sin x and c = cos x) to convert a non-polynomial system into a polynomial system with no approximation to bound. The generated code does neither; it stops at the Taylor step.
auto_R flag, and so the same script switches back to the single fixed disk by setting it false. The SOS Synthesis tab places this beside the alternating method it replaces.pip install SumOfSquares)