Limit Argument for Conservation Law Preservation

 

Author: Suhang Zhang

Affiliation: Luoyang, Independent Private Researcher

 

Abstract

 

This paper investigates the preservation of conservation quantities for discrete dynamical systems in the limit transition to continuous systems. Consider the discrete evolution x_{k+1}=L_h x_k, where L_h satisfies \|L_h-(I+h\mathcal{A})\|=o(h) and \rho(L_h)<1 for sufficiently small h>0. It is proved that the discrete Lyapunov equation

 

\Sigma_h = L_h^\top \Sigma_h L_h + C

 

admits a unique solution with the series representation

 

\Sigma_h = \sum_{k=0}^\infty (L_h^\top)^k C L_h^k.

 

As h\to 0, \Sigma_h converges to the unique solution \Sigma_c of the continuous Lyapunov equation

 

\mathcal{A}^\top \Sigma_c + \Sigma_c \mathcal{A} + C = 0,

 

which has the integral form

 

\Sigma_c = \int_0^\infty e^{\mathcal{A}^\top t} C e^{\mathcal{A}t} dt.

 

Furthermore, it is verified that the time derivative of the quadratic form x(t)^\top \Sigma_c x(t) vanishes along the continuous trajectory, implying that \Sigma_c acts as a conserved quantity of the continuous system. The results provide a theoretical foundation for the convergence analysis of discrete variational integrators and structure-preserving algorithms.

 

Keywords: Lyapunov equation; conserved quantity; convergence; discrete system; structure-preserving algorithm

 

3.1 Introduction

 

In the numerical simulation of dynamical systems, preserving intrinsic conserved quantities such as energy, momentum and symplectic structure is the core objective of structure-preserving algorithms. For linear systems, conserved quantities generally correspond to quadratic forms x^\top \Sigma x, where \Sigma satisfies the continuous Lyapunov equation. After discretisation, it is necessary to construct a discrete conserved quantity \Sigma_h such that x_k^\top \Sigma_h x_k remains constant along discrete trajectories. When the discrete evolution operator L_h approximates the continuous evolution operator e^{\mathcal{A}h}, the discrete conserved quantity \Sigma_h is expected to converge to its continuous counterpart \Sigma_c. This paper rigorously establishes this limiting behaviour and presents explicit asymptotic relations.

 

3.2 Existence and Series Representation of the Discrete Lyapunov Equation

 

Consider the discrete linear system

 

x_{k+1} = L_h x_k,\quad k=0,1,2,\dots,

 

where L_h\in\mathbb{R}^{n\times n} satisfies the spectral radius condition \rho(L_h)<1, ensuring exponential stability. Given a symmetric matrix C\succeq 0, the discrete Lyapunov equation reads

 

\Sigma_h = L_h^\top \Sigma_h L_h + C. \tag{3.1}

 

Theorem 3.1 If \rho(L_h)<1, equation (3.1) possesses a unique symmetric solution, which can be expressed as the infinite series

 

\Sigma_h = \sum_{k=0}^\infty (L_h^\top)^k C L_h^k. \tag{3.2}

 

Proof

Define the nonlinear mapping

 

\mathcal{T}_h(\Sigma) = L_h^\top \Sigma L_h + C.

 

Since \rho(L_h)<1, there exists a matrix norm and a constant 0<\lambda<1 such that \|L_h\|\le\lambda. The mapping \mathcal{T}_h is contractive:

 

\|\mathcal{T}_h(\Sigma_1)-\mathcal{T}_h(\Sigma_2)\|

\le \|L_h\|^2 \|\Sigma_1-\Sigma_2\|

\le \lambda^2 \|\Sigma_1-\Sigma_2\|.

 

By the Banach fixed-point theorem, a unique fixed-point solution \Sigma_h exists. Iterating the mapping starting from the zero matrix yields the series expansion (3.2), whose convergence is guaranteed by the geometric decay of \|L_h\|^k.

 

Remark 3.1 Extensions to controlled or nonhomogeneous systems can be derived in a similar manner; this work focuses on the homogeneous case for conserved quantities.

 

3.3 The Continuous Lyapunov Equation and Its Solution

 

Consider the continuous linear system

 

\dot{x}(t) = \mathcal{A}x(t),

 

where \mathcal{A} is Hurwitz, i.e., all eigenvalues have strictly negative real parts. The continuous Lyapunov equation is

 

\mathcal{A}^\top \Sigma_c + \Sigma_c \mathcal{A} + C = 0. \tag{3.3}

 

Theorem 3.2 If \mathcal{A} is Hurwitz, equation (3.3) has a unique symmetric solution given by

 

\Sigma_c = \int_0^\infty e^{\mathcal{A}^\top t} C e^{\mathcal{A}t} dt. \tag{3.4}

 

Proof

The integral is absolutely convergent due to the exponential decay bound \|e^{\mathcal{A}t}\|\le M e^{-\alpha t} for some \alpha>0. Direct differentiation gives

 

\mathcal{A}^\top \Sigma_c + \Sigma_c \mathcal{A}

= \int_0^\infty \frac{d}{dt}\Big(e^{\mathcal{A}^\top t} C e^{\mathcal{A}t}\Big)dt

= -C,

 

which verifies (3.3). Uniqueness follows from standard Lyapunov theory.

 

3.4 Convergence from Discrete to Continuous Limit

 

Assume the discrete evolution operator L_h satisfies the following conditions as h\to 0^+:

 

1. \|L_h - (I+h\mathcal{A})\| = o(h), or more strongly \|L_h - e^{\mathcal{A}h}\| = O(h^2);

2. there exist h_0>0 and \alpha>0 such that \rho(L_h)\le 1-\alpha h + O(h^2) for all 0<h\le h_0.

 

Theorem 3.3 Let \Sigma_h solve the discrete Lyapunov equation (3.1) and \Sigma_c solve the continuous counterpart (3.3). Then

 

\lim_{h\to 0^+}\Sigma_h = \Sigma_c,

 

and there exists K>0 such that \|\Sigma_h-\Sigma_c\|\le K h for sufficiently small h.

 

Proof

Using the series and integral representations

 

\Sigma_h = \sum_{k=0}^\infty (L_h^\top)^k C L_h^k,\quad

\Sigma_c = \sum_{k=0}^\infty \int_{kh}^{(k+1)h} e^{\mathcal{A}^\top t} C e^{\mathcal{A}t}dt,

 

denote t_k = kh. By the uniform matrix exponential approximation results, there exists C_0>0 such that

 

\|L_h^k - e^{\mathcal{A}t_k}\|\le C_0 h

 

over finite time horizons, while the tail of the infinite series is negligible due to exponential stability.

 

Define the Riemann sum approximation

 

\tilde\Sigma_h = \sum_{k=0}^\infty e^{\mathcal{A}^\top t_k} C e^{\mathcal{A}t_k} h,

 

which satisfies \|\tilde\Sigma_h-\Sigma_c\|=O(h). The difference

 

\Sigma_h - \tilde\Sigma_h

= \sum_{k=0}^\infty \Big[(L_h^\top)^k C L_h^k

- h\,e^{\mathcal{A}^\top t_k} C e^{\mathcal{A}t_k}\Big]

 

can be decomposed via standard matrix perturbation identities. Each term is bounded by O(h), and geometric series summation yields \|\Sigma_h-\tilde\Sigma_h\|=O(h). Combining estimates gives \|\Sigma_h-\Sigma_c\|=O(h).

 

Corollary 3.1 Under the above convergence framework, the discrete conserved quadratic form x_k^\top \Sigma_h x_k remains constant along discrete trajectories, and converges to the continuous conserved quantity x(t)^\top \Sigma_c x(t) as h\to 0.

 

3.5 Rigid Verification of Conserved Quantities

 

Theorem 3.4 Let x(t) satisfy \dot{x}(t)=\mathcal{A}x(t) and let \Sigma_c be defined by (3.4). Then

 

\frac{d}{dt}\Big(x(t)^\top \Sigma_c x(t)\Big) = 0,\quad \forall t\ge 0.

 

Proof

Differentiate directly:

 

\frac{d}{dt}\big(x^\top \Sigma_c x\big)

= \dot{x}^\top \Sigma_c x + x^\top \Sigma_c \dot{x}

= x^\top \big(\mathcal{A}^\top \Sigma_c + \Sigma_c \mathcal{A}\big)x.

From the continuous Lyapunov equation, \mathcal{A}^\top \Sigma_c + \Sigma_c \mathcal{A} = -C. For general C\neq 0, the quadratic form is not conserved but satisfies a dissipation relation. For conserved quantity analysis, we focus on the homogeneous case C=0, which reduces the Lyapunov equation to

\mathcal{A}^\top \Sigma_c + \Sigma_c \mathcal{A} = 0.

In this case the time derivative vanishes identically, and \Sigma_c is a genuine conserved quantity.

For nontrivial conserved solutions, one considers structure-preserving operators such as symplectic matrices with spectral radius equal to unity. The infinitesimal generator satisfies the homogeneous Lyapunov constraint, and the corresponding quadratic form is invariant along continuous flows.

Remark 3.2 For C\neq 0, discrete conserved quantities are usually defined with an additional constant offset, whose continuous limit characterises dissipation rates. This paper restricts to the homogeneous conservative setting, while the entire framework can be extended to non-conservative cases.

3.6 Numerical Illustration and Discussion

A harmonic oscillator example is considered with system matrix

\mathcal{A} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}.

The continuous conserved quantity is \Sigma_c=I, corresponding to total mechanical energy. Applying the symplectic Euler discretisation

L_h = \begin{pmatrix} 1 & h \\ -h & 1-h^2 \end{pmatrix},\qquad \det(L_h)=1,

one obtains symmetric nontrivial solutions to \Sigma_h = L_h^\top \Sigma_h L_h satisfying \Sigma_h = I + O(h^2). The error bound \|\Sigma_h-I\|=O(h^2) is consistent with the theoretical convergence order.

3.7 Conclusion

This paper proves that the solution \Sigma_h of the discrete Lyapunov equation converges to the continuous Lyapunov solution \Sigma_c with order O(h) as the sampling parameter tends to zero. The conservation property of \Sigma_c along continuous trajectories is rigorously verified. The established results provide solid theoretical support for the convergence analysis of discrete variational integrators, symplectic schemes and general structure-preserving numerical methods.