Learning Lyapunov Operators for Nonlinear Systems
arXiv:2609.18894v1 [math.AP] 16 Sep 2026
Amartya Mukherjee, Maxwell Fitzsimmons, David C. Del Rey Fernández, and Jun Liu
Abstract— Constructing Lyapunov functions for nonlinear dynamical systems is a central problem in stability analysis, yet remains challenging. Lyapunov functions are commonly characterized as solutions to first-order partial differential equations (PDEs), but these solutions are typically obtained for single systems, limiting their reuse across systems. In this paper, we study the Lyapunov solution operator that maps a vector field to the corresponding Lyapunov function defined by a dissipation-based Lyapunov PDE. We establish that, on compact subsets of the domain of attraction and under exponential stability assumptions, this operator is well-defined, unique, and continuous with respect to perturbations of both the vector field and the dissipation function. These results provide a theoretical foundation for approximating Lyapunov functions uniformly over families of nonlinear systems. Building on these theoretical foundations, we employ Fourier Neural Operators (FNOs) as a data-driven approximation of the Lyapunov solution operator. Numerical experiments demonstrate that a single trained operator can accurately approximate the numerical Lyapunov functions across parameterized families of dynamics. This illustrates the potential of neural operators for approximating Lyapunov functions.
I. INTRODUCTION One of the longstanding challenges in nonlinear systems and control is the construction of Lyapunov functions [4]. While Lyapunov functions can essentially be characterized by solutions to partial differential equations (PDEs) and neural network solutions to such PDEs can effectively provide approximations to Lyapunov functions [9], solving PDEs for each system can still be time-consuming. Additionally, prior works that use neural networks to learn Lyapunov functions only aid in the verification of a single system [11], thus failing to generalize into systems with slightly different dynamics. This could pose difficulties in real-life systems. Recently, operator learning has emerged as a paradigm for approximating mappings between function spaces, such as those defined by PDEs. Unlike conventional deep learning architectures that learn finite-dimensional mappings, neural operators generalize across function spaces, allowing them to learn solution operators of PDEs from data. This perspective is particularly attractive for Lyapunov analysis: the correspondence between vector fields and their associated Lyapunov functions can be viewed as an operator defined by a PDE constraint. Approximating this operator directly offers This work was supported in part by the Natural Sciences and Engineering Research Council of Canada and the Canada Research Chairs Program. Amartya Mukherjee, Maxwell Fitzsimmons, David C. Del Rey Fernández, and Jun Liu are with the Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1 (email: [email protected] (Jun Liu)).
the promise of learning a single model that can compute Lyapunov functions for a large class of systems. The Fourier Neural Operator (FNO) is an example of a neural operator [5], [7]. It lifts the input functions to a higher-dimensional feature space using a linear layer, followed by repeated applications of Fourier convolution layers. Each layer performs a fast Fourier transform (FFT), applies learnable filters in the frequency domain, and then inverts the transform to return to the spatial domain. This global convolution mechanism enables FNO to efficiently capture long-range dependencies in the input. In this paper, we bridge the perspectives of Lyapunov stability and operator learning. We begin by formulating the Lyapunov stability condition as a PDE, thereby recasting the problem into the operator learning framework. We then establish regularity and continuity results for this PDE, showing that the assumptions required for the FNO universality theorem [6] hold in our setting. This provides theoretical justification for approximating Lyapunov operators using FNOs. We finally train an FNO on nonlinear systems and demonstrate that it can generate Lyapunov functions with small error with respect to the true Lyapunov function. II. P ROBLEM F ORMULATION We consider autonomous nonlinear dynamical systems of the form (1) ẋ(t) = f (x(t)), x ∈ Rn , where the vector field f : Rn → Rn is continuously differentiable and satisfies f (0) = 0. Let φ f (t, x) denote the flow of the system initialized at x at time t = 0. A. Admissible class of dynamics Define the function space F0 := { f ∈ C1 (Rn ; Rn ) | f (0) = 0}, and the subclass of locally exponentially stable vector fields ∂f (0) is Hurwitz . Fst := f ∈ F0 ∂x For f ∈ Fst , classical results guarantee that the origin is a locally exponentially stable equilibrium and that solutions exist and are unique in a neighborhood of the origin. The domain of attraction of the origin is defined as DOA( f ) := {x ∈ Rn | lim φ f (t, x) = 0}. t→∞
(2)
A set Ω is called positively invariant if φ f (t, x) ∈ Ω for all t ≥ 0 and x ∈ Ω. B. Dissipation functions and Lyapunov PDE
∇ω(0) = 0,
∇2 ω ≻ 0,
(3)
where ∇2 ω denotes the Hessian matrix. Define the admissible class W>0 := ω ∈ C2 (Rn , R) ω satisfies conditions (3) . Given f ∈ Fst and ω ∈ W>0 , we consider Lyapunov functions defined as solutions of the first-order PDE ∇V (x) · f (x) = −ω(x),
V (0) = 0.
(4)
Equations of this form arise naturally in converse Lyapunov theory and characterize dissipation-based Lyapunov functions. In applications, local Lyapunov functions are found to prove local asymptotic/exponential stability, but they are also used to estimate the size of DOA( f ). This is because if V is a Lyapunov function for f on K = {x ∈ Rn : V (x) < r} for a fixed r, then K is positively invariant and K ⊂ DOA( f ). Moreover, there is a Lyapunov function V1 for f such that {x ∈ Rn : V1 (x) < r} = DOA( f ) [15], [16]. C. Integral representation and well-posedness Let K ⊂ DOA( f ) be compact and positively invariant. Under this assumption, the Lyapunov PDE (4) admits the integral representation Z ∞ V f ,ω (x) := ω φ f (t, x) dt, x ∈ K, (5) 0
which is well defined due to the exponential decay of trajectories and the regularity of ω. Proposition 1 (Existence and uniqueness): Let f ∈ Fst , ω ∈ W>0 , and let K ⊂ DOA( f ) be compact and positively invariant. Then the following are equivalent: 1) There exists a unique function V ∈ C1 (K) satisfying (4) 2) This solution is given by the integral formula (5). Proof: If (1) holds, then we see that d (V f ,ω ◦ φ f (t, x)) = −ω ◦ φ f (t, x) dt for t ≥ 0 and x ∈ K. Integrating from 0 to T , we find V f ,ω ◦ φ f (T, x) −V f ,ω (x) =
Z T 0
R
V f ,ω ◦ φ f (T, x) −V f ,ω (x) = −T ω ◦ φ f (s′ , x)
Let ω : Rn → R be a twice continuously differentiable function satisfying ω(0) = 0,
If (2) holds then V f ,ω ◦ φ f (T, x) − V f ,ω (x) = 0T −ω ◦ φ f (s, x)ds. By the mean value theorem for integrals, we obtain
−ω ◦ φ f (s, x)ds.
We note that, as ω is Lipschitz, the integrated solution is unique. Taking the limit as T → ∞ yields the result, after noting that φ f (T, x) → 0 if x ∈ K and V (0) = 0.
for s′ ∈ [0, T ]. Dividing both sides by T and taking the limit T → 0 yields the result. Finally, V is continuously differentiable on K as proved by [9]. As a consequence, V f ,ω is a Lyapunov function for f on K, and its sublevel sets define invariant subsets of the domain of attraction. D. The Lyapunov solution operator Rather than constructing Lyapunov functions on a per-system basis, we adopt an operator-theoretic viewpoint. Define the Lyapunov solution operator G : Fst × W>0 ,
G ( f , ω) := V f ,ω ,
mapping a vector field and dissipation function to the corresponding Lyapunov function. The central objective of this work is to establish continuity of this operator on compact subsets of the domain of attraction. E. Sobolev regularity of the data and solution The universality results for FNOs are formulated for operators acting between Sobolev spaces, while Proposition 1 defines solutions to the Lyapunov PDE in the C1 and C2 space. To place the Lyapunov solution operator within this framework, it is necessary to ensure that both the input and output functions admit sufficient Sobolev regularity. The Sobolev regularity results here are standard in the analysis of PDEs and play a central role in neural operator theory [1]. They ensure that the Lyapunov PDE (4) is well defined pointwise while simultaneously allowing us to work in function spaces compatible with neural operator approximation. The next two theorems ensure that the solutions to the Lyapunov PDE can be embedded in appropriate Sobolev spaces. Theorem 1 (Sobolev embedding [13]): Let Ω ⊂ Rn be a bounded open domain with Lipschitz boundary, and let k ≥ 0 be an integer. If s > n/2 + k, then the Sobolev space H s (Ω) is continuously and compactly embedded in Ck (Ω), i.e., H s (Ω) ,→ Ck (Ω). Moreover, there exists a constant C > 0 such that ∥u∥Ck (Ω) ≤ C∥u∥H s (Ω) ,
∀u ∈ H s (Ω).
Corollary 1 (Sobolev embedding for vector-valued functions): Let Ω ⊂ Rn be a bounded open domain with Lipschitz boundary, and let k ≥ 0 be an integer. If s > n/2 + k, then the Sobolev space H s (Ω; Rn ) consisting of Rn -valued functions with s-Sobolev regularity is continuously and compactly embedded in Ck (Ω; Rn ), i.e., H s (Ω; Rn ) ,→ Ck (Ω; Rn ). Moreover, there exists a constant C > 0 such that ∥u∥Ck (Ω;Rn ) ≤ C∥u∥H s (Ω;Rn ) ,
∀u ∈ H s (Ω; Rn ).
The proof of Corollary 1 follows directly from the scalar case by applying the component-wise argument. III. C ONTINUITY AND U NIVERSALITY OF THE LYAPUNOV S OLUTION O PERATOR
A. FNO universality theorem The work of [6] provides the key conditions imposed on the input functions and output function of an operator so that the FNO universality theorem applies. In our setting, the operator is the Lyapunov operator, which parametrizes the Lyapunov PDE over the vector field f and the dissipation function ω. To apply the universality result, we must verify the following conditions: 1) the input vector fields belong to a Sobolev space H s (Ω; Rn ) and the dissipation functions belong to a Sobolev space H r (Ω) with sufficiently high regularity; 2) the corresponding Lyapunov functions belong to a ′ Sobolev space H s (Ω); 3) the Lyapunov solution operator is continuous on compact subsets of the product of the input spaces. Under these assumptions, the modified universality theorem of [6] for the Lyapunov PDE yields convergent approximation by FNOs. We restrict attention to compact subsets of the admissible classes Fst and W>0 , viewed as subsets of appropriate Sobolev spaces via Sobolev embedding. Theorem 2 (Modification of Theorem 9 by [6]): Let Ω ⊂ Rn be a bounded domain with Lipschitz boundary such that Ω ⊂ (0, 2π)n . Let s, r, s′ ≥ 0, and let K f ⊂ Fst ∩ H (Ω; R ), n
Kω ⊂ W>0 ∩ H (Ω) r
be compact sets of admissible vector fields and dissipation functions, respectively. Assume that for each ( f , ω) ∈ K f × Kω , the Lyapunov PDE ∇V (x) · f (x) = −ω(x),
x ∈ Ω,
(6)
s′
admits a unique solution V f ,ω ∈ H (Ω), and that the associated Lyapunov solution operator ′
G : K f × Kω → H s (Ω),
G ( f , ω) = V f ,ω ,
is continuous with respect to the product topology induced by ′ the H s (Ω; Rn ) and H r (Ω) norms on the input and the H s (Ω) norm on the output. Let Tn denote the n-dimensional torus. Then, for every ε > 0, there exist 1) continuous linear extension operators E f : H s (Ω; Rn ) → H s (Tn ; Rn ),
′
Nε : H s (Tn ; Rn ) × H r (Tn ) → H s (Tn ), such that
In this section, we establish the main theoretical result of the paper: continuity of the Lyapunov solution operator with respect to perturbations of the vector field and the dissipation function. This property is essential for approximating Lyapunov functions uniformly over families of nonlinear systems.
s
2) a Fourier Neural Operator
Eω : H r (Ω) → H r (Tn ),
sup ( f ,ω)∈K f ×Kω
G ( f , ω) − Nε (E f f , Eω ω) Ω H s′ (Ω) < ε.
In particular, the Lyapunov solution operator can be approximated arbitrarily well on compact subsets of admissible vector fields and dissipation functions by a suitable FNO. The regularity conditions required for this theorem are satisfied in our setting. Specifically, from the Sobolev embedding results in Section II-E, for s > n/2 +1 we have H s (Ω; Rn ) ,→ C1 (Ω; Rn ), ensuring that vector fields are continuously differentiable. For the dissipation functions, we require r > n/2+2 to guarantee H r (Ω) ,→ C2 (Ω), which is the natural regularity for ω in the Lyapunov PDE. The continuity of the Lyapunov solution operator G is established in Theorem 3, providing the final ingredient needed to apply the universality result. B. Regularity of the Lyapunov PDE Before turning to continuity, we briefly comment on regularity. Under Section II-E, the Lyapunov PDE admits a unique solution V f ,ω ∈ C1 (J) on compact invariant sets J ⊂ DOA( f ). In numerical settings, vector fields and Lyapunov functions are often represented in Sobolev spaces. For sufficiently large s > n/2+1, classical Sobolev embedding results ensure continuous and compact embeddings H s (J; Rn ) ,→ C1 (J; Rn ),
H s (J) ,→ C1 (J),
which guarantee that the Lyapunov PDE is well defined in the function spaces used by neural operator architectures. C. Main continuity result We begin by stating the central theorem. Theorem 3 (Continuity of the Lyapunov solution operator): Let f ∈ Fst and ω ∈ W>0 . Let K ⊂ DOA( f ) be compact. Then there exists a compact set, K ⊂ J ⊂ DOA( f ), such that the Lyapunov solution operator G : Fst × W>0 → C1 (J),
G (g, ψ) := Vg,ψ ,
is continuous at ( f , ω) when Fst × W>0 is endowed with the norm ∥g∥C1 (J) + ∥ψ∥C2 (J) , and C1 (J) is endowed with the uniform norm ∥ · ∥∞,J . Remark 1: Theorem 3 formalizes the intuition that small perturbations of the system dynamics and dissipation function induce small changes in the associated Lyapunov function on compact invariant sets. This robustness property
provides the theoretical foundation for learning Lyapunov functions uniformly over families of nonlinear systems. The proof of Theorem 3 proceeds in three steps. First, we establish robustness of exponential stability under perturbations of the vector field. Second, we show uniform decay of the Lyapunov function outside small sublevel sets. Finally, we prove uniform integrability of the Lyapunov integral representation. D. Robust exponential stability We now establish robustness of local exponential stability with respect to perturbations of the vector field. Lemma 1 (Robust exponential stability): Let f ∈ Fst , ω ∈ W>0 . There exist constants r > 0, δ > 0, M > 0, and c > 0 such that for any g ∈ F0
∥g − f ∥C1 (J) < r,
with
the origin remains locally exponentially stable for ẋ = g(x) and ∥φg (t, x)∥ ≤ Me−ct ∥x∥, ∀t ≥ 0, ∥x∥ ≤ δ . Proof: Define the ball B(x, η) := {y ∈ Rn : ∥x−y∥ < η}, and define B(x, η) as its closure. Choose η > 0 such that B(0, η) ⊂ DOA( f ). Let J = K ∪ B(0, η). Let A := ∂∂ xf (0), Q = ∇2 ω(0). Since f ∈ Fst , A is Hurwitz, and since ω ∈ W>0 , Q ≻ 0. Let P ∈ Rn×n be the positive definite matrix that satisfies [4] PA + AT P = −Q. Let λmin (·) and λmax (·) refer to the smallest and largest eigenvalues respectively. Define W (x) := 12 xT Px, then W is a local Lyapunov function for f . Let λmin (Q) c1 := . 6λmax (P) We will show that, for sufficiently small δ > 0 and r > 0, the function W is a strict local Lyapunov function for the system ẋ = g(x) whenever ∥g − f ∥C1 (J) < r.
Since xT PAx = 12 xT (PA + AT P)x = − 12 xT Qx, we obtain d W (φg (t, x)) dt t=0 1 ∂g ∂ f ≤ − λmin (Q)∥x∥2 + ∥P∥ ∥x∥2 − 2 ∂x ∂x ∞ ∂f + ∥P∥ − A ∥x∥2 ∂x ∞ 1 λmin (Q) 2 ≤ − λmin (Q)∥x∥ + ∥P∥r∥x∥2 + ∥P∥ ∥x∥2 2 6∥P∥ 1 1 1 ≤ − λmin (Q)∥x∥2 + λmin (Q)∥x∥2 + λmin (Q)∥x∥2 2 6 6 1 2 = − λmin (Q)∥x∥ . 6 Using 2W (x) ≤ λmax (P)∥x∥2 , it follows that λmin (Q) d ≤− W (φg (t, x)) W (x) = −2cW (x) dt 3λmax (P) t=0 for all ∥x∥ ≤ δ1 . Therefore, along any trajectory of ẋ = g(x) that remains in B(0, δ1 ), d W (φg (t, x)) ≤ −2cW (φg (t, x)). dt By Grönwall’s inequality, W (φg (t, x)) ≤ e−2ct W (x),
In particular, since W is decreasing, the sublevel set Ωδ := {x ∈ Rn : W (x) ≤ 12 λmin (P)δ 2 } is positively invariant whenever δ ≤ δ1 . Choose δ = δ1 . Then every trajectory with ∥x∥ ≤ δ remains in B(0, δ ) ⊂ B(0, η) ⊂ J for all t ≥ 0, so the above estimate is valid globally in time. Finally, using the bounds relating W and ∥x∥2 , 2 W (φg (t, x)) λmin (P) 2 ≤ e−2ct W (x) λmin (P) λmax (P) −2ct ≤ e ∥x∥2 . λmin (P)
∥φg (t, x)∥2 ≤
Since f ∈ C1 , there exists δ1 > 0 such that δ1 ≤ η and ∂f (ξ ) − A ≤ c1 . ∥ξ ∥ ≤ δ1 =⇒ ∂x Choose r := c1 . Now let g ∈ F0 satisfy ∥g − f ∥C1 (J) < r, and let ∥x∥ ≤ δ1 . Since g(0) = 0, Taylor’s theorem gives
Hence, s
Z 1
g(x) =
∂g (tx)dt x 0 ∂x
Hence, d W (φg (t, x)) = xT Pg(x) dt t=0 Z 1 ∂g = xT P (tx)dt x 0 ∂x
t ≥ 0.
∥φg (t, x)∥ ≤
λmax (P) −ct e ∥x∥. λmin (P)
Setting s M :=
λmax (P) , λmin (P)
we obtain ∂ f ∂ g = xT PAx + xT P (tx) − (tx) dt x ∥φg (t, x)∥ ≤ Me−ct ∥x∥, ∀t ≥ 0, ∥x∥ ≤ δ . ∂x ∂x 0 Z 1 ∂f This proves local exponential stability of the origin for ẋ = T +x P (tx) − A dt x. g(x). ∂x 0 Z 1
E. Uniform decay of the Lyapunov derivative We next show that the Lyapunov function V f ,ω decreases uniformly along trajectories outside small sublevel sets. Lemma 2 (Uniform Lyapunov decay): Let f ∈ Fst , ω ∈ W>0 , and let J ⊂ DOA( f ) be compact. Define U(δ ) := {x ∈ Rn : V f ,ω (x) < δ }. Then there exists δ0 > 0 such that for every δ ∈ (0, δ0 ), there exist constants cδ > 0 and rδ > 0 such that for all g ∈ F0
with
Next, we establish a global lower bound on ω outside U(δ ). To separate our analysis into large and small level sets, we decompose the set K2 \U(δ ) = K2 \ B(0, δ1 ) ∪ B(0, δ1 ) \U(δ ) . Since ω is continuous and strictly positive, the minimum away from the origin, m1 :=
min
x∈K2 \B(0,δ1 )
ω(x),
is strictly positive. Combining this with (8), we obtain ( ) λmin (Q) δ2 min ω(x) ≥ min m1 , . 2 x∈K2 \U(δ ) ∥∇V f ,ω ∥2∞,K2
∥g − f ∥C1 (J) < rδ ,
the inequality ∇V f ,ω (x) · g(x) ≤ −cδ
Define
holds for all x ∈ J \U(δ ).
( ) 1 λmin (Q) δ2 cδ := min m1 , > 0. 2 2 ∥∇V f ,ω ∥2∞,K2
Let Q := ∇2 ω(0) ≻ 0, and denote λmin (Q) > 0
Proof: its smallest eigenvalue.
Then,
Let R := max V f ,ω (x), x∈J
min
K2 := {x ∈ Rn : V f ,ω (x) ≤ R}.
x∈K2 \U(δ )
Then K2 is compact, positively invariant under f , and satisfies J ⊆ K2 ⊆ DOA( f ). We first establish a lower bound on ω near the origin. Since ω ∈ C2 and ∇2 ω(0) = Q ≻ 0, there exists δ1 > 0 such that for all ∥x∥ ≤ δ1 , ∥∇2 ω(x) − Q∥ ≤
Let g ∈ F0 satisfy ∥g − f ∥C1 (J) < rδ . Then, for all x ∈ K2 \ U(δ ), ∇V f ,ω (x) · g(x) = −ω(x) + ∇V f ,ω (x) · (g(x) − f (x)) ≤ −2cδ + ∥∇V f ,ω ∥∞,K2 ∥g − f ∥∞,K2
λmin (Q) . 2
< −2cδ + ∥∇V f ,ω ∥∞,K2 rδ = −cδ .
1 ω(x) = xT ∇2 ω(ξ )x 2
This proves the result.
for some ξ on the segment between 0 and x. Hence, ω(x) = xT Qx + xT ∇2 ω(ξ ) − Q x
F. Uniform integrability of the Lyapunov integral We now establish convergence and robustness of the integral representation.
≥ λmin (Q)∥x∥2 − ∥∇2 ω(ξ ) − Q∥∥x∥2 λmin (Q) ∥x∥2 . ≥ 2
(7)
Since V f ,ω is continuous and V f ,ω (0) = 0, there exists δ0 > 0 such that for all δ ∈ (0, δ0 ), U(δ ) ⊂ B(0, δ1 ).
Lemma 3 (Uniform integrability): Let f ∈ Fst and ω ∈ W>0 . There exist constants δ > 0, r > 0, and C > 0 such that for all g ∈ F0 , ψ ∈ W>0
with
the integral
Fix such a δ , then for x ∈ B(0, δ1 ) \U(δ ), we have V f ,ω (x) ≥ δ . By Taylor’s theorem, followed by Cauchy-Schwarz inequality, ∥x∥ ≥
δ ∥∇V f ,ω ∥∞,K2
.
x ∈ B(0, δ1 ) \U(δ ).
0
ψ(φg (t, x)) dt
is well defined for all x ∈ J, and T
|ψ(φg (t, x))| dt ≤ Cδ
for all x ∈ J and all sufficiently large T .
Combining with (7), we obtain λmin (Q) δ2 ω(x) ≥ , 2 ∥∇V f ,ω ∥2∞,K2
∥g − f ∥C1 (J) + ∥ψ − ω∥C2 (J) < r, Z ∞
Vg,ψ (x) =
Z ∞
=⇒
(9)
Finally, we prove the main perturbation argument. Let cδ . rδ := ∥∇V f ,ω ∥∞,K2
By Taylor’s theorem, for such x,
δ ≤ ∥∇V f ,ω ∥∞,K2 ∥x∥
ω(x) ≥ 2cδ .
(8)
Proof: Fix δ > 0 sufficiently small as in Lemma 2, and define U(δ ) := {x ∈ Rn : V f ,ω (x) < δ }.
Let x ∈ K2 and define the hitting time
Fix ε > 0. We will show that for (g, ψ) sufficiently close to ( f , ω) in C1 (J) ×C2 (J),
T (δ , x, g) := inf{t ≥ 0 : V f ,ω (φg (t, x)) ≤ δ }.
∥Vg,ψ −V f ,ω ∥∞,J < ε.
From Lemma 2, we have d V f ,ω (φg (t, x)) ≤ −cδ dt
For x ∈ J, we use the integral representation Z ∞ V f ,ω (x) −Vg,ψ (x) = ω(φ f (t, x)) − ψ(φg (t, x)) dt.
for t ∈ [0, T (δ , x, g)].
Integrating and using V f ,ω (φg (T (δ , x, g), x)) = δ , we obtain T (δ , x, g) ≤
V f ,ω (x) − δ . cδ
|ω(φ f (t, x)) − ψ(φg (t, x))| ≤ |ω(φ f (t, x)) − ω(φg (t, x))|
Since V f ,ω is bounded on K2 , there exists a constant Tδ > 0 such that T (δ , x, g) ≤ Tδ for all x ∈ K2 . By construction, U(δ ) ⊂ B(0, δ1 ) for sufficiently small δ . From Lemma 1, there exist constants M1 > 0 and c1 > 0 such that for all t ≥ T (δ , x, g), −c1 (t−T (δ ,x,g))
∥φg (t, x)∥ ≤ M1 e
0
We decompose the integrand: + |ω(φg (t, x)) − ψ(φg (t, x))|. Fix δ > 0 and define Tδ := supx∈J T (δ , x, g), which is finite by Lemma 3. For t ∈ [0, Tδ ], the flows φ f and φg remain in J. By Grönwall inequality, there exists C1 > 0 such that ∥φ f (t, x) − φg (t, x)∥ ≤ C1 ∥ f − g∥C1 (J) .
sup x∈J,t∈[0,Tδ ]
∥φg (T (δ , x, g), x)∥.
Since ω ∈ C1 (J), it is Lipschitz on J, so there exists Lω > 0 such that
Since φg (T (δ , x, g), x) ∈ U(δ ), we have ∥φg (T (δ , x, g), x)∥ ≤ δ ,
|ω(φ f (t, x)) − ω(φg (t, x))| ≤ Lω ∥φ f (t, x) − φg (t, x)∥.
and therefore
Hence,
∥φg (t, x)∥ ≤ M1 δ e−c1 (t−T (δ ,x,g)) ,
t ≥ T (δ , x, g).
(10)
Z T
δ
0
lies in a bounded subset of C2 (K2 ), its gradient
Since ψ is uniformly bounded on K2 . Hence, there exists a constant L > 0 such that
for some constant C2 > 0. Similarly, Z T
|ψ(y)| ≤ L∥y∥,
δ
∀ y ∈ K2 .
|ω(φ f (t, x)) − ω(φg (t, x))| dt ≤ C2 ∥ f − g∥C1 (J)
0
|ω(φg (t, x)) − ψ(φg (t, x))| dt ≤ Tδ ∥ω − ψ∥∞,J .
For t ≥ T (δ , x, g), Lemma 3 yields
Using (10), for t ≥ T (δ , x, g),
Z ∞
|ψ(φg (t, x))| ≤ L∥φg (t, x)∥ ≤ LM1 δ e−c1 (t−T (δ ,x,g)) . Therefore, Z ∞ T (δ ,x,g)
T (δ ,x,g)
|ψ(φg (t, x))| dt ≤ Cδ .
Applying the same argument to ( f , ω), |ψ(φg (t, x))| dt ≤ LM1 δ
Z ∞
e−c1 s ds =
0
LM1 δ. c1
1 Define C := LM c1 . Then for all x ∈ K2 ,
Z ∞ T (δ ,x,g)
|ψ(φg (t, x))| dt ≤ Cδ .
Since T (δ , x, g) is uniformly bounded over K2 , the integral defining Vg,ψ (x) is finite for all x ∈ J.
Proof: Let K ⊂ DOA( f ) be compact. Define x∈K
T (δ ,x, f )
|ω(φ f (t, x))| dt ≤ Cδ .
Thus, the tail contribution satisfies Z ∞ Tδ
|ω(φ f (t, x)) − ψ(φg (t, x))| dt ≤ 2Cδ .
Combining the estimates, we obtain ∥V f ,ω −Vg,ψ ∥∞,J ≤ C2 ∥ f − g∥C1 (J) + Tδ ∥ω − ψ∥∞,J + 2Cδ . First choose δ > 0 such that 2Cδ < ε/3. Then choose (g, ψ) sufficiently close to ( f , ω) so that
G. Proof of Theorem 3
R := max V f ,ω (x),
Z ∞
J := {x ∈ Rn : V f ,ω (x) ≤ R}.
Then J is compact, positively invariant under f , and satisfies K ⊂ J ⊂ DOA( f ).
C2 ∥ f − g∥C1 (J) < ε/3,
Tδ ∥ω − ψ∥∞,J < ε/3.
It follows that ∥Vg,ψ −V f ,ω ∥∞,J < ε. Therefore, G is continuous at ( f , ω).
IV. N UMERICAL R ESULTS
the predicted and the true Lyapunov function and divides by the L1 norm of the true Lyapunov function. This is the default evaluation metric used in [7], [2].
A. Dataset Generation We constructed datasets for three representative nonlinear dynamical systems: the damped Duffing oscillator, the inverted pendulum, and the Van der Pol oscillator. For each system, 1,000 parameterized instances were generated, with parameters sampled from uniform distributions as specified below. Lyapunov functions were obtained by numerically solving the PDE V̇ = −xT Qx using the LyZNet toolbox [8]. Q is a randomly sampled positive definite matrix. Each dataset was divided into 800 training, 100 validation, and 100 test samples. Each system was projected onto the grid [−1, 1]2 divided into a 64x64 grid during training of the FNO.
Table I summarizes the L1 errors on the test set. FNO achieved a lower error (0.0182) compared to DeepONet (0.6483), demonstrating its improved ability to approximate Lyapunov functions. FNO 0.0182
DeepONet 0.6483
TABLE I: L1 errors of FNO and DeepONet on learning the Lyapunov PDE
D. Visualization of Learned Functions •
Damped Duffing Equation: ẋ1 = x2 ,
(11)
ẋ2 = −δ x2 − αx1 − β x13 ,
(12)
with parameters sampled as α ∼ U(1, 10), β ∼ U(0.1, 2.0), δ ∼ U(0.1, 1.0). • Inverted Pendulum: θ̇1 = θ2 ,
(13)
Figure 1 illustrates test cases across the systems. For each case, we plot the learned Lyapunov function against the ground-truth solution. The FNO reconstructions closely follow the true solutions, capturing its structure. V. C ONCLUSION
B. Model Implementation
In this work, we introduced a framework for learning Lyapunov functions of nonlinear dynamical systems using FNOs By formulating the Lyapunov condition as a PDE and leveraging the universality of FNOs on Sobolev spaces, we established theoretical guarantees ensuring that Lyapunov operators can be approximated with high fidelity. Our analysis demonstrated regularity and continuity properties that justify the application of neural operator learning in this setting. Through numerical experiments on nonlinear systems, we showed that FNOs achieve substantially lower approximation error compared to DeepONets.
We use the implementation of FNO by [7]. For our experiments, we adopt an extended version of the standard FNO by integrating Adaptive Instance Normalization (AdaIN, [3]) to handle conditioning on the parameters Q. These parameters are flattened and passed through a small MLP to compute normalization constants that modulate the intermediate features at various stages of the network,
A key challenge is extending this framework to higherdimensional dynamical systems. In principle, the universality of FNOs extends naturally to Rn , but practical implementation quickly becomes computationally prohibitive. Even for 3D problems, Fast Fourier Transforms (FFTs) scale as O(N 3 log N) in time and O(N 3 ) in memory, where N is the resolution along each dimension [14]. This growth imposes severe memory and runtime bottlenecks during training.
c θ̇2 = −cθ2 − sin θ1 , l with parameters l ∼ U(0.1, 100), c ∼ U(0.1, 10) • Van Der Pol:
(14)
ẋ1 = −x2 ,
(15)
ẋ2 = x1 − µx2 (1 − x12 ),
(16)
with parameters µ ∼ U(0.01, 1.0).
AdaIN(x, α(Q), β (Q)) = α(Q)
x − µ(x) + β (Q), σ (x)
(17)
where µ(x) and σ (x) are the mean and standard deviation of the layer, and α(Q), β (Q) are trainable MLPs that map the matrix Q to scalars that determine scaling and shifting in the normalized layer. C. Performance Comparison We compared the FNO with DeepONet [10] in learning solutions to the Lyapunov PDE. We assessed the performance of these models by evaluating them on the testing set using the relative L1 error, which computes the L1 error between
Finally, combining learned Lyapunov functions with controller design remains a promising avenue. By embedding Lyapunov certificates into feedback synthesis, one may learn controllers with formal guarantees. This would bridge datadriven stability analysis with practical control implementation, further motivating the development of scalable and interpretable neural operator architectures. An early step in this direction was accomplished using diffusion models [12]. R EFERENCES [1] L. C. Evans. Partial differential equations, volume 19. American mathematical society, 2022.
Fig. 1: Comparison of learned Lyapunov functions (output V ) against ground truth solutions (true V ) for representative test cases across the inverted pendulum (Row 1), damped Duffing oscillator (Row 2), and Van der Pol oscillator (Row 3). Each subplot shows the input vector field components ( f1 , f2 ) and the corresponding Lyapunov function values. The FNO reconstructions closely match the ground truth across all systems, closely matching the numerical reference solutions.
[2] M. Herde, B. Raonic, T. Rohner, R. Käppeli, R. Molinaro, E. de Bézenac, and S. Mishra. Poseidon: Efficient foundation models for pdes. Advances in Neural Information Processing Systems, 37:72525–72624, 2024. [3] X. Huang and S. Belongie. Arbitrary style transfer in real-time with adaptive instance normalization. In Proceedings of the IEEE international conference on computer vision, pages 1501–1510, 2017. [4] H. K. Khalil. Nonlinear Systems, Third Edition. Pearson Education, 2002. [5] J. Kossaifi, N. Kovachki, Z. Li, D. Pitt, M. Liu-Schiaffini, R. J. George, B. Bonev, K. Azizzadenesheli, J. Berner, and A. Anandkumar. A library for learning neural operators, 2024. [6] N. Kovachki, S. Lanthaler, and S. Mishra. On universal approximation and error bounds for fourier neural operators. Journal of Machine Learning Research, 22(290):1–76, 2021. [7] N. B. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. M. Stuart, and A. Anandkumar. Neural operator: Learning maps between function spaces. CoRR, abs/2108.08481, 2021. [8] J. Liu, Y. Meng, M. Fitzsimmons, and R. Zhou. Tool lyznet: A lightweight python tool for learning and verifying neural lyapunov functions and regions of attraction. In Proceedings of the 27th ACM International Conference on Hybrid Systems: Computation and Control, pages 1–8, 2024. [9] J. Liu, Y. Meng, M. Fitzsimmons, and R. Zhou. Physics-informed
neural network lyapunov functions: Pde characterization, learning, and verification. Automatica, 175:112193, 2025. [10] L. Lu, P. Jin, and G. E. Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv:1910.03193, 2019. [11] Y. Meng, R. Zhou, A. Mukherjee, M. Fitzsimmons, C. Song, and J. Liu. Physics-informed neural network policy iteration: Algorithms, convergence, and verification. In Forty-first International Conference on Machine Learning, 2024. [12] A. Mukherjee, T. Quartz, and J. Liu. Manifold-guided stabilization of nonlinear dynamical systems with diffusion models. In 2025 American Control Conference (ACC), pages 3850–3855. IEEE, 2025. [13] S. L. Sobolev. Sur un théorème d’analyse fonctionnelle. Recueil Mathématique (Nouvelle série), 4(46):471–497, 1938. [14] C. Van Loan. Computational frameworks for the fast Fourier transform. SIAM, 1992. [15] A. Vannelli and M. Vidyasagar. Maximal Lyapunov functions and domains of attraction for autonomous nonlinear systems. Automatica, 21(1):69–80, 1985. [16] V. I. Zubov. Methods of A. M. Lyapunov and Their Application. Noordhoff, 1964.