Available online at [Link].

com

ScienceDirect
Journal of Differential Equations 453 (2026) 113865
[Link]/locate/jde

Exponential stability preserving of two spatially

discretized port-Hamiltonian systems ✩
Fu Zheng a,∗ , Hongjian Yin a , Zhongjie Han b , Bao-Zhu Guo c
a School of Mathematics and Statistics, Hainan University, Haikou, Hainan, 570228, China
b School of Mathematics, Tianjin University, Tianjin, 300072, China
c Key Laboratory of Systems and Control, Academy of Mathematics and Systems Science, Academia Sinica, Beijing
100190, China
Received 12 September 2024; revised 31 August 2025; accepted 14 October 2025

Abstract
For an ideal transmission line described by the telegrapher’s equations, a mixed finite element method-an
extension of widely used spatially discretized approach-has been introduced. This numerical approximation
approach maintains both the Dirac structure and passivity, guaranteeing that the spatially discretized system
preserves its port-Hamiltonian characteristics. In this paper, we employ this method to spatially discretize
two infinite-dimensional port-Hamiltonian systems characterized by variable coefficients and boundary con­
trols. Subsequently, we explore the preservation of exponential stability in the resulting semi-discretized
systems, establishing their uniform exponential stability concerning discretization parameters. Through fre­
quency domain analysis, uniform exponential stability is demonstrated for both semi-discretized models.
Finally, numerical simulations confirm the efficacy of this semi-discrete approach.
© 2025 Elsevier Inc. All rights are reserved, including those for text and data mining, AI training, and
similar technologies.

✩
This work was supported by the National Natural Science Foundation of China under grant Nos. 12371446,
U23B2033, 62473281; and Hainan Provincial National Natural Science Foundation of Hainan under grant
No. 123MS004; and Key Laboratory of Engineering Modeling and Statistical Computation of Hainan Province; and
the Specific Research Fund of the Innovation Platform for Academicians of Hainan Province; and Scientific Research
Initiation Fund of Hainan University (RZ2200001240).
* Corresponding author.
E-mail address: fuzheng@[Link] (F. Zheng).

[Link]
0022-0396/© 2025 Elsevier Inc. All rights are reserved, including those for text and data mining, AI training, and similar
technologies.
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Keywords: Port-Hamiltonian system; Exponential stabilization; Mixed finite element; Semi-discretization; Frequency
domain

1. Introduction

It is universally acknowledged that control systems modeled by partial differential equations

(PDEs) are intrinsically infinite-dimensional. A fundamental challenge in simulating and con­
trolling these systems is achieving finite-dimensional approximations. To tackle this, various
numerical approximation methods have been developed within the realm of PDE numerical com­
putation, with the goal of preserving system structures. A common initial approach involves
spatially discretizing the PDEs while keeping the time variable continuous, a process known as
semi-discretization. This approach offers several benefits. Firstly, the resulting semi-discretized
models are finite-dimensional, enabling their analysis and design by control engineers. Secondly,
if these models retain many of the physical properties of the original PDEs, they can be consid­
ered physical modes rather than mere approximations of the infinite-dimensional system. As a
result, substantial research over the past two decades has been dedicated to preserving geomet­
ric structure [3], [14], [6], passivity [14,8], [6], controllability [11], [15], observability [24], and
exponential stability [18,8,12,23,17,22,7].
The pursuit of maintaining exponential stability in numerical approximations can be traced
back to the 1980s, when Gibson et al. [5,4] explored linear quadratic optimal control. Con­
currently, Banks et al. [2] delved into numerical approximations for linear quadratic regulator
problems and solutions to operator Riccati equations for distributed parameter systems with un­
bounded input and boundary control, emphasizing the crucial role of exponential stability. In
1990, Banks et al. [1] initiated research on exponentially stable approximations for wave equa­
tions with boundary dampers, observing that widely used discretization techniques such as finite
differences or finite elements often fail to preserve a uniform decay rate in semi-discretized sys­
tems. Subsequently, Infante et al. [9] provided a mathematical explanation for this issue, showing
that spurious high-frequency modes generated during spatial semi-discretization can undermine
this preservation.
An early significant effort involved researchers modifying the finite difference scheme to
recover uniform exponential stability. For instance, Tebou and Zuazua introduced an artificial
numerical viscosity term into the finite difference semi-discretization of the 1-D wave equation
with boundary control [18], thereby restoring uniform exponential stability in the new discrete
systems. Similar strategies have been employed to achieve uniform controllability [11] and
uniformly exponentially stable approximations for a class of abstract second-order evolution
equations [16]. More recently, Liu and Guo developed a novel semi-discretized order reduction
finite difference scheme to explore uniform approximation of the 1-D wave equation [12]. They
reformulated the wave equation as a first-order port-Hamiltonian system and applied the central
finite difference scheme. These semi-discretized systems demonstrated uniform exponential sta­
bility without additional measures. Using this semi-discretization scheme, they investigated the
uniform exponential stability of the wave equation with local viscosity and observer-based con­
trol [23], the Schrödinger equation [7], the Timoshenko beam [20], and the heat-wave coupled
system [22].
However, while the finite difference scheme offers numerous benefits for handling PDEs with
constant coefficients, examining the uniform exponential stability of its semi-discretization for
PDEs with variable coefficients presents challenges, often making it difficult to validate key

2
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

conclusions [21]. Therefore, it is crucial to explore alternative semi-discretization schemes for

PDEs with variable coefficients that are suitable for studying uniform exponential stability.
To the best of our knowledge, the mixed finite element method may be the most appropriate
approximation scheme for the ideal transmission line described by the telegrapher’s equations
[6]. This is supported by three primary reasons. Firstly, the spatially discretized system retains
its port-Hamiltonian nature, and the geometric structure of the spatially discretized system on
a partition interval mirrors that of the entire spatial domain, greatly simplifying the analysis
of discrete dynamics. Secondly, the discretized scheme derived from this mixed finite element
method shares similarities with the order reduction finite difference method proposed in [12].
Specifically, the temporal derivative component of the spatially discretized system corresponds
to a weighted average of two adjacent node functions (see Remark 2.1 below), which contrasts
with the temporal derivative component of the order reduction finite difference method. This
similarity allows for the application of extensive research on uniform exponential stability ap­
proximations for the order reduction finite difference method to mixed finite elements. Finally,
despite the significant challenge of identifying an appropriate energy multiplier, the frequency
domain characterization of uniform exponential stability for a family of contractive semigroups
[13][20][7] can be applied to this scheme.
In this paper, we conduct a study of both the ideal transmission line and the Timoshenko
beam. The paper is organized as follows. Section 2 revisits the spatial discretization method based
on mixed finite elements for the ideal transmission line. Section 3 explores the preservation of
exponential stability in the resulting semi-discretized system for the ideal transmission line, with
a focus on ensuring its uniform exponential stability with respect to discretization parameters.
Section 4 introduces a novel semi-discretized scheme for the Timoshenko beam and examines its
exponential stability preservation through frequency domain analysis. Finally, Section 5 presents
several numerical simulations to validate the effectiveness of our theoretical findings.

2. Revisiting the results

2.1. Continuous PDE system

The telegrapher’s equations for the transmission line with a boundary resistor are given by
PDEs with variable coefficients:
⎧ ∂q(t, z) ∂I (t, z) ∂ϕ(t, z) ∂V (t, z)
⎪
⎨ =− , =− ,
∂t ∂z ∂t ∂z (2.1)
⎩ V (t, 0) = 0, V (t, S) = RI (t, S),
⎪
q(0, z) = q0 (z), ϕ(0, z) = ϕ0 (z),

where q(·, ·) denotes the charge density, ϕ(·, ·) is the flux density, and the current I (·, ·) and
voltage V (·, ·) are defined as

ϕ(t, z) q(t, z)
I (t, z) = , V (t, z) = (2.2)
L(z) C(z)

with C(z) and L(z) representing the distributed capacitance and distributed inductance of the
line, respectively. The spatial variable z belongs to the domain ℐ := [0, S] with S > 0, and the
voltage at z = 0 is set to zero, while a resistor R > 0 is placed at the other end, resulting in
V (t, S) = RI (t, S).

3
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

The energy of the system (2.1) is expressed as

∫︂
q(t, z)2 ϕ(t, z)2
E(t) = + dz. (2.3)
C(z) L(z)
ℐ

To guarantee the exponential decay of energy along the solution to equation (2.1), we require
the following assumption.

Assumption 2.1. Assume that C(·), L(·) ∈ L∞ ∞

p (ℐ), which means that C(·), L(·) ∈ L (ℐ) and
there are positive constants m and M such that the inequalities m ≤ C(z), L(z) ≤ M hold.

Theorem 2.1. Given Assumption 2.1 and the condition that q0 , ϕ0 ∈ L2 (ℐ), the system (2.1) has
a unique solution in the state space X := L2 (ℐ; R2 ). Moreover, there are two positive constants
M0 and ω0 such that the energy along the solution to (2.1) satisfies the following inequality:

E(t) ≤ M0 e−ω0 t E(0), (2.4)

Proof. The proof for Theorem 2.1 closely parallels the approach taken in the proof of Exam­
ple 9.2.1 presented in [10]. Due to this strong similarity, we choose to omit the detailed proof
here. □

Drawing inspiration from the concepts and notations in [6], the telegrapher’s equations (2.1)
can be reformulated in a geometric version as

∂q(t, z) ∂ϕ(t, z)
= −deϕ (t, z), = −deq (t, z), (2.5)
∂t ∂t

where q(t, z), ϕ(t, z) are interpreted as one-forms, and eϕ (t, z), eq (t, z) are zero-forms given by
∫︂ [︃ ]︃
∗q(z) ∗ϕ(z)
eq = δq H, eϕ = δϕ H, H = q(z) + ϕ(z) , (2.6)
2C(z) 2L(z)
ℐ

with δ representing the variational derivative.

Here, we give some explanations on the terminologies zero-form, one-form, the exterior
derivative d and the Hodge star operator ∗. In the context of a one-dimensional spatial domain,
a function is a zero-form, while the one-form can be integrated over every sub-interval of the
interval. If we consider a spatial coordinate z for some interval, then a function f is simply given
by the values f (z) ∈ R for every coordinate value z in the interval, while a one-form is given as
g(z)dz for a certain density function g.
Given a coordinate z for the spatial domain we obtain by spatial differentiation of a function
f (z) the one-form g := dfdz (z)dz. In coordinate-free language this is denoted as g = df , where d
is called the exterior derivative, converting zero-forms into one-forms.
The Hodge star operator ∗ generally converts any k-form ρ on a n-dimensional spatial do­
main to an (n − k)-form ∗ρ for 0 ≤ k ≤ n. The definition of the Hodge star operator relies on
the assumption of a Riemannian metric on the spatial domain. However, in the present paper
this Riemannian metric will simply be the Euclidean inner product corresponding to a choice

4
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

of local coordinates. Thus on the one-dimensional spatial domain with spatial coordinate z we
simply have ∗f (z) = f (z)dz and ∗(g(z)dz) = g(z). The introduction of these notations is both
convenient for the readers and intended to lay the groundwork for future discussions on high­
dimensional systems.

2.2. Semi-discretization of a part of the transmission line

This subsection is adapted from select parts of [6] to ensure the self-contained nature of the
paper. Initially, a semi-discretization process is applied to a segment of the transmission line
spanning between two points a and b (0 ≤ a < b ≤ S), with ℐab = [a, b].
Step 1: Approximations of one-forms on the domain ℐab .
On the interval ℐab , the energy variables q(t, z) and ϕ(t, z), as well as the infinitesimal charge
rate ∂q(t, z)/∂t and the infinitesimal flux rate ∂ϕ(t, z)/∂t, are approximated as
{︃ q
q(t, z) = Qab (t)ωab (z),
ϕ (2.7)
ϕ(t, z) = Φab (t)ωab (z),

and
⎧
⎪
⎨
∂q(t, z) q q
= fab (t)ωab (z),
∂t (2.8)
⎪
⎩
∂ϕ(t, z) ϕ ϕ
= fab (t)ωab (z),
∂t
q ϕ
where one-forms ωab (z) and ωab (z) satisfy
∫︂ ∫︂
q ϕ
ωab (z) = ωab (z) = 1 (2.9)
ℐab ℐab

and
dQab (t) q dΦab (t) ϕ
= fab (t), = fab (t). (2.10)
dt dt
Step 2: Approximations of zero-forms on the domain ℐab .
The co-energy variables eq (t, z) and eϕ (t, z) are approximated on ℐab according to the fol­
lowing expressions:
{︃ q q q q
eq (t, z) = ea (t)ωa (z) + eb (t)ωb (z),
ϕ ϕ ϕ ϕ (2.11)
eϕ (t, z) = ea (t)ωa (z) + eb (t)ωb (z),

q q ϕ ϕ
where zero-forms ωa (z), ωb (z) and ωa (z), ωb (z) satisfy the boundary value conditions:

q q q q
ωa (a) = 1, ωa (b) = 0, ωb (a) = 0, ωb (b) = 1,
ϕ ϕ ϕ ϕ
ωa (a) = 1, ωa (b) = 0, ωb (a) = 0, ωb (b) = 1,

as well as the compatibility of the forms, given by

5
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

ϕ q q q ϕ
−dωaϕ = dωb = ωab , −dωa = dωb = ωab . (2.12)

Substituting (2.8), (2.11) and (2.12) into (2.5) gives

{︄
q q q ϕ q
fab (t)ωab (z) = eaϕ (t)ωab (z) − eb (t)ωab (z),
ϕ ϕ q ϕ q ϕ (2.13)
fab (t)ωab (z) = ea (t)ωab (z) − eb (t)ωab (z).

By integrating these identities over the interval ℐab and utilizing expression (2.9), we obtain

q ϕ ϕ q q
fab (t) = eaϕ (t) − eb (t), fab (t) = ea (t) − eb (t). (2.14)

Step 3: Approximates of Hamiltonian and dynamics

On the one hand, by employing the terminology associated with Dirac structures, the power
supplied to or taken from the electrical part of the part of the transmission line dHdt
ab (t)
can be
q ϕ
expressed as the summation of the products between the flow variables fab (t), fab (t), and the
q ϕ
corresponding effort/co-energy variables eab (t), eab (t). By substituting expressions (2.8) and
(2.11) into the continuous power
∫︂ [︃ ]︃
dHab (t) ∂q(t, z) q ∂ϕ(t, z) ϕ
= ∗e (t, z) + ∗e (t, z)
dt ∂t ∂t
ℐab

and applying Proposition 1 from [6], we arrive at

dHab (t) q q q ϕ ϕ
= [αab ea (t) + (1 − αab )eb (t)]fab (t) + [(1 − αab )eaϕ (t) + αab eb (t)]fab (t)
dt

where
∫︂
q q
αab = wab (z)∗wa (z).
ℐab

Thus, we are able to define the co-energy variables as follows:

{︃ q q q
eab (t) := αab ea (t) + (1 − αab )eb (t),
ϕ ϕ ϕ (2.15)
eab (t) := (1 − αab )ea (t) + αab eb (t),

q q ϕ ϕ
such that dH ab (t)/dt = fab (t)eab (t) + fab (t)eab (t).
On the other hand, if we define

∫︂ q q ∫︂ ϕ ϕ
−1 wab (z) ∗ wab (z) wab (z) ∗ wab (z)
Cab = , L−1
ab = ,
C(z) L(z)
ℐab ℐab

then through (2.6)-(2.7) the energy of the considered segment can be approximated by

6
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Q2ab (t) Φ2ab (t)

Hab (t) = + .
2Cab 2Lab

By the variational derivative, the discrete co-energy variables are again derived through the ex­
pression of Hab (t):
⎧
⎪ q Qab (t)
⎨ eab (t) = δQab Hab (t) = ,
Cab (2.16)
⎪ Φ (t)
⎩ eϕ (t) = δΦab Hab (t) = ab .
ab
Lab
q ϕ
Finally, substituting and eliminating the co-energy variables eab (t) and eab (t), the discrete­
time dynamics of system (2.1) over the interval ℐab can be derived from equations (2.10),
(2.14)-(2.16). The resulting discrete dynamics are expressed as follows:
⎧
⎪ d
⎨ Cab [αab eaq (t) + (1 − αab )eq (t)] = eaϕ (t) − eϕ (t),
b b
dt (2.17)
⎪ d
⎩ Lab [(1 − αab )ea (t) + αab e (t)] = ea (t) − eq (t).
ϕ ϕ q
b b
dt

Remark 2.1. Setting hab = |ℐab | as the length of the interval ℐab and dividing both sides of the
identities in (2.17) by hab , we obtain
⎧ ϕ ϕ
⎪
⎪ C d e (t) − ea (t)
⎨ ab q q
[αab ea (t) + (1 − αab )eb (t)] = − b ,
hab dt hab
q q
⎪
⎪ ab
L d e (t) − ea (t)
⎩ ϕ
[(1 − αab )eaϕ (t) + αab eb (t)] = − b .
hab dt hab
q q
Cab [αab ea (t)+(1−αab )eb (t)]
In the above equation, the convex combinations of hab and
ϕ ϕ
Lab [(1−αab )ea (t)+αab eb (t)]
hab can be interpreted as approximations of q(t, z) and ϕ(t, z) on ℐab ,
q q ϕ ϕ
e (t)−e (t) e (t)−e (t)
respectively. Similarly, b hab a and b hab a represent approximations of the derivatives
de (t, z) and de (t, z) on ℐab . This approach is analogous to the finite difference scheme pro­
q ϕ

posed in [21]. For brevity, (2.17) will be retained for further discussion.

The discrete energy associated with the considered segment of the transmission line is given
by

1 [︂ q q ϕ
]︂
Hab (t) = Cab (αab ea (t) + (1 − αab )eb (t))2 + Lab ((1 − αab )eaϕ (t) + αab eb (t))2 , (2.18)
2
and it possesses the following noteworthy property.

Proposition 2.1. The discrete energy defined in (2.18) satisfies

dHab (t) q q ϕ
= ea (t)eaϕ (t) − eb (t)eb (t). (2.19)
dt
7
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Proof. Starting from (2.17) and (2.18), we have

dHab (t) q q ϕ
= [αab ea (t) + (1 − αab )eb (t)][eaϕ (t) − eb (t)]
dt
ϕ q q
+ [(1 − αab )eaϕ (t) + αab eb (t)][ea (t) − eb (t)].

By simple calculation, (2.19) is derived from the above identity. □

2.3. Semi-discretization of the whole transmission line

To achieve the semi-discretization of the entire transmission line, for any positive integer n, we
insert n − 1 points Sj (with j = 1, · · · , n − 1) into the spatial domain ℐ. This results in a partition
of ℐ into segments ℐSj −1 Sj for j = 1, · · · , n, with S0 = 0, Sn = S and ℐSj −1 Sj = [Sj −1 , Sj ], and
ℐSj −1 Sj = [Sj −1 , Sj ]. On each segment ℐSj −1 Sj , we select

q ∗C(z) ϕ ∗L(z)
ωSj −1 Sj = ∫︁ , ωSj −1 Sj = ∫︁
ℐSj −1 Sj C(z)dz ℐSj −1 Sj L(z)dz

as one-form finite elements, along with their corresponding compatible zero-form finite elements
q q ϕ ϕ
ωSj −1 , ωSj , ωSj −1 , and ωSj to carry out the semi-discretization process described in the previous
q q ϕ ϕ
subsection. For brevity, we set ej (t) = eSj (t), ej (t) = eSj (t), and define
∫︂ ∫︂ ∫︂
q q
Cj = C(z)dz, Lj = L(z)dz, αj = ωSj −1 Sj ∗ωSj −1
ℐSj −1 Sj ℐSj −1 Sj ℐSj −1 Sj

with ˆ︁
αj = 1 − αj for j = 0, 1, · · · , n. Subsequently, we obtain the semi-discretized approxima­
tion of (2.5) on ℐSj −1 Sj as follows:
⎧
⎪ d
⎪ q q ϕ ϕ
⎨ Cj dt [αj ej −1 (t) + ˆ︁
⎪ αj ej (t)] = ej −1 (t) − ej (t),
d ϕ ϕ q q (2.20)
⎪
⎪ Lj [ˆ︁ αj ej −1 (t) + αj ej (t)] = ej −1 (t) − ej (t),
⎪ q dt
⎩ q
e0 (t) = 0, en (t) = Renϕ (t), j = 0, 1, · · · , n.

The final two equations in (2.20) are derived directly from the boundary conditions specified in
(2.1). Furthermore, the discrete energy is given by

1 ∑︂ [︂ ]︂
n
q q ϕ ϕ
H△n (t) = Cj (αj ej −1 (t) + ˆ︁
αj ej (t))2 + Lj (ˆ︁
αj ej −1 (t) + αj ej (t))2 . (2.21)
2
j =1

Remark 2.2. The discrete scheme (2.20) represents a generalization of two well-known semi­
discretization approaches. Specifically, when C(z) and L(z) are set to constant values, and αj =
1, the scheme reduces to the finite difference method on staggered grids as presented in [19].
On the other hand, if C(z) = L(z) = 1 and αj is set to 12 , then (2.20) corresponds to the order­
reduced finite difference scheme described in [12,23].

8
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Now, let us define some notations and outline certain assumptions.

Assumption 2.2. Let hj := |ℐSj −1 Sj | = Sj − Sj −1 denote the length of each segment, and let
△n := max1≤j ≤n hj be the maximum segment length. We assume that △n < 1 for all n ∈ N +
and that △n = 𝒪(n−1 ), meaning there exists a positive constant C such that △n ≤ 𝒞n−1 for all
n ∈ N+.

Assumption 2.3. Under the Assumptions 2.1 and 2.2, we assume the existence of positive con­
stants c and C such that the coefficients αj satisfy

c ≤ min αj ≤ max αj ≤ C < 1

1≤j ≤n 1≤j ≤n

for all j = 1, 2, · · · , n and n ∈ N + .

q q
Remark 2.3. It is straightforward to verify that the function ωSj −1 (z), defined as ωSj −1 (z) :=
Sj −z q q
Sj −Sj −1 , satisfies the boundary conditions ωSj −1 (Sj −1 ) = 1 and ωSj −1 (Sj ) = 0. Noting that
ωSj −1 Sj = ∫︁ ∗C(z)
q
C(z)dz
, we proceed with the following derivation:
ℐ Sj −1 Sj

∫︂ ∫︂
q q M Sj − z
αj = ωSj −1 Sj ∗ ωSj −1 ≤ dz
m Sj − Sj −1
ℐSj −1 Sj ℐSj −1 Sj
M M
= (Sj − Sj −1 ) = hj = 𝒪(n−1 ).
2m 2m
This demonstrates that the inequality max1≤j ≤n αj ≤ C < 1 holds true. The other inequality can
be established through a similar line of reasoning.

Assumption 2.4. Assumptions 2.1 and 2.2 imply that max1≤j ≤n Cj = 𝒪(n−1 ), min1≤j ≤n Cj =
𝒪(n−1 ), max1≤j ≤n Lj = 𝒪(n−1 ) and min1≤j ≤n Lj = 𝒪(n−1 ).

3. Uniform exponential stability of (2.20)

Firstly, based on (2.18) and Proposition 2.1, the discrete energy H△n (t) defined in (2.21)
satisfies the following balance equation and dissipative property.

Proposition 3.1. The balance equation and dissipative property, given by

dH△n (t) q ϕ q
= e0 (t)e0 (t) − en (t)enϕ (t) = −R|enϕ (t)|2 (3.1)
dt
hold for all positive integers n.

Proof. Indeed, we can express the discrete energy H△n (t) as the sum of local contributions:

∑︂
n
H△n (t) = HSj −1 Sj (t)
j =1

9
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

where each local energy HSj −1 Sj (t) is given by:

1 [︂ q q ϕ ϕ
]︂
HSj −1 Sj (t) = Cj (αj eSj −1 (t) + ˆ︁
αj eSj (t))2 + Lj (ˆ︁
αj eSj −1 (t) + αj eSj (t))2 .
2

By applying Proposition 2.1 to each HSj −1 Sj (t), we obtain

dHSj −1 Sj (t) q ϕ q ϕ
= eSj −1 (t)eSj −1 (t) − eSj (t)eSj (t).
dt

Therefore, summing over all segments gives

dH△n (t) ∑︂ dHSj −1 Sj (t) ∑︂ [︂ q ]︂

n n
ϕ q ϕ q ϕ q
= = eSj −1 (t)eSj −1 (t) − eSj (t)eSj (t) = e0 (t)e0 (t) − en (t)enϕ (t).
dt dt
j =1 j =1

Finally, substituting the boundary conditions from (2.20) yields the dissipative property stated in
(3.1). This completes the proof of the proposition. □

Secondly, to analyze (2.20) from a frequency domain perspective, it is necessary to reformu­

late the equation into suitable state equations defined on a particular state space. To this end, we
define the state space 𝒳△n = C 2n endowed with the inner product:

⟨︁ ⟩︁ ∑︂
n
X△n , Y△n △n
= Cj (αj xj −1 + ˆ︁
αj xj )(αj yj −1 + ˆ︁
αj yj )
j =1

∑︂
n
+ αj xn+j + αj xn+j +1 )(ˆ︁
Lj (ˆ︁ αj yn+j + αj yn+j +1 ), (3.2)
j =1

where X△n = (x1 , · · · , x2n ) and Y△n = (y1 , · · · , y2n ) ∈ 𝒳△n are elements of 𝒳△n . Note that in
the definition of the inner product, additional terms like x0 , y0 , x2n+1 and y2n+1 appear. These
are assigned specific values: x0 = y0 = 0, x2n+1 = kxn and y2n+1 = kyn with k = R −1 > 0 to
unify the notations of the right hand-side of inner product (3.2). In the sequel, this treatment
is utilized several times similarly. However, given that x0 = 0 is a known value and x2n+1 can
be determined based on the relationship x2n+1 = kx2n once xn is established by the differential
equation (2.20), there is no need to incorporate x0 and x2n+1 into the state space 𝒳Δ .
Now, we recast (2.20) into a vectorial form for clarity and convenience. To do so, we introduce
the vectors W△n (t) = (e1 (t), · · · , en (t))⊤ and V△n (t) = (e0 (t), · · · , en−1 (t))⊤ as the unknown
q q ϕ ϕ

variables of (2.20). Additionally, we define the n × n matrices D△n , D ˆ︁△n and M△n as follows:

⎛ ⎞ ⎛ ⎞
ˆ︁
α1 ˆ︁
α1 α1
⎜ α2 ⎟ ⎜ .. ⎟
⎜ ˆ︁
α2 ⎟ ⎜ ˆ︁
α2 . ⎟
D△n =⎜ ⎟, ˆ︁△n = ⎜
D ⎟,
⎝ .. .. ⎠ ⎜ .. ⎟
. . ⎝ . αn−1 ⎠
αn ˆ︁
αn ˆ︁
αn

10
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

⎛ ⎞
1 −1
⎜ .. ⎟
⎜ 1 . ⎟
M△ n =⎜
⎜
⎟.
⎟ (3.3)
⎝ ..
. −1 ⎠
1

We further introduce the n × n diagonal matrices

𝒞△n = diag{C1 , C2 , · · · , Cn }, ℒ△n = diag{L1 , L2 , · · · , Ln }

and ℰ△n = diag{0, 0, · · · , 1} to facilitate the vectorial representation of (2.20). Consequently, the
equation (2.20) can be rewritten in the concise form:

′ ′
ℋ△n Φ△n X△ n
(t) = Ψ△n X△n (t) ⇔ X△ n
(t) = A△n X△n (t), (3.4)

where A△n = Φ−1 −1

△n ℋ△n Ψ△n , and the vector X△n (t) combines both W△n (t) and V△n (t) as

(︃ )︃
W△n (t)
X△n (t) = .
V△n (t)

The matrices involved in this representation are defined as follows:

(︃ )︃ (︃ )︃ (︃ )︃
𝒞△n D△n 0 −kℰ△n M△ n
ℋ△n = , Φ△n = ˆ︁△n , Ψ △n = ⊤ .
ℒ △n kαn ℰ△n D −M△ n
0

It is noteworthy that Φ△n serves as a weighting operator, reflecting the influence of αj and ˆ︁αj .
If we substitute the weighting operator Φ△n with the identity operator I , then the equation (3.4)
simplifies to the classical finite difference scheme of (2.1).
Furthermore, the inner product defined in (3.2) within the space 𝒳△n can be elegantly ex­
pressed as

⟨︁ ⟩︁ ⟨︁ ⟩︁
X△n , Y△n 𝒳 = Φ△n X△n , H△n Φ△n Y△n , ∀X△n , Y△n ∈ 𝒳△n ,
△n

where X△n = (x1 , · · · , x2n ), Y△n = (y1 , · · · , y2n ) ∈ 𝒳△n . Here ⟨·, ·⟩ represents the standard inner
product in C 2n .

Lemma 3.1. The operator A△n is dissipative and generates a C0 -semigroup of contractions
T△n (t) on the state space 𝒳△n .

Proof. By setting A△n X△n = Y△n with the additional conditions x0 = y0 = 0, x2n+1 = kxn and
y2n+1 = kyn , for j = 1, · · · , n and any X△n ∈ 𝒳△n , we obtain

Cj (αj yj −1 + ˆ︁
αj yj ) = xn+j − xn+j +1 , Lj (ˆ︁
αj yn+j + αj yn+j +1 ) = xj −1 − xj .

11
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Subsequently, utilizing the definition of the inner product ⟨·, ·⟩△n , we derive:

2Re⟨A△n X△n , X△n ⟩△n = ⟨A△n X△n , X△n ⟩△n + ⟨X△n , A△n X△n ⟩△n
∑︂
n
= 2Re Cj (αj xj −1 + ˆ︁
αj xj )(αj yj −1 + ˆ︁
αj yj )
j =1

∑︂
n
+2Re αj xn+j + αj xn+j +1 )(ˆ︁
Lj (ˆ︁ αj yn+j + αj yn+j +1 )
j =1

∑︂
n ∑︂
n
= 2Re (αj xj −1 + ˆ︁
αj xj )(xn+j − xn+j +1 ) + 2Re αj xn+j + αj xn+j +1 )(xj −1 − xj )
(ˆ︁
j =1 j =1

∑︂
n ∑︂
n
= 2Re [xn+j xj −1 − xj xn+j +1 ] = 2Re [xn+j xj −1 − xn+j +1 xj ]
j =1 j =1

= 2Re[xn+1 x0 − x2n+1 xn ] = −2k|xn | . 2

(3.5)

Equation (3.5) directly implies that the operator A△n is dissipative. The second statement is self­
evident, thus completing the proof of this lemma. □

Remark 3.1. In the reviewing process, one reviewer proposed an enlightening proof of the
Lemma 3.1. More precisely, if setting Λ△n to be the diagonal matrix consisting of α1 , · · · , αn ,
then it is easy to see that

D △ n = I − Λ △ n M△
T ˆ︁△n = I − Λ△n MΔn .
,D
n

Using this we find that

(︄ )︄
T
Λ △ n M△ 0
Φ△n = I − n .
−kαn ℰ Λ△n MΛ△n

Thus, it is easy to see that

⟨A△n X△n , X△n ⟩△n + ⟨X△n , A△n X△n ⟩△n = ⟨Ψ△n X△n , Φ△n X△n ⟩ + ⟨Φ△n X△n , Ψ△n X△n ⟩,

and some basic calculations give

ΨT△n Φ△n + ΦT△n Ψ△n

(︃ )︃ [︃ (︃ )︃]︃
−kℰ△n −M△n T
Λ △ n M△ 0
= T I − n
M△ n
0 −kαn ℰ△n Λ△n M△n
[︄ (︃ )︃T ]︄ (︃ )︃T
T
Λ △ n M△ 0 −kℰ△n −M△n
+ I− n
T
−kαn ℰ△n Λ△n M△n M△ n
0

12
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

(︃ )︃
−2R −1 ℰ 0
= .
0 0

These show that the operator A△n is dissipative.

Lemma 3.2. The intersection of the spectral set σ (A△n ) of the operator A△n with the imaginary
axis is empty.

Proof. Given X△n ∈ 𝒳△n and s ∈ R, by setting A△n X△n = isX△n with the additional conditions
x0 = 0 and x2n+1 = kxn , we obtain from (3.4)

isCj (αj xj −1 + ˆ︁
αj xj ) = xn+j − xn+j +1 , isLj (ˆ︁
αj xn+j + αj xn+j +1 ) = xj −1 − xj . (3.6)

From (3.5) and the fact that A△n X△n = isX△n , we have

0 = Re⟨A△n X△n , X△n ⟩△n = −k|xn |2 . (3.7)

Setting j = n in (3.6) and using (3.7) along with x2n+1 = kxn = 0, we get

isCj αj xn−1 − x2n = 0, −xn−1 + isLj ˆ︁

αj x2n = 0.

This implies that xn−1 = x2n = 0 since the determinant of the coefficient matrix
⃓ ⃓
⃓ isCj αj −1 ⃓⃓
⃓ = −s 2 αj ˆ︁
αj Cj Lj − 1
⃓ −1 αj ⃓
isLj ˆ︁

is nonzero.
Similarly, by setting j = n − 1 in (3.6) and using the fact that xn−1 = x2n = 0, we can deduce
that xn−2 = x2n−1 = 0. By induction on j , we can show that xj = 0 for all j = 1, · · · , 2n.
This implies that X△n = 0, which contradicts the assumption that X△n is a nonzero eigenvector
corresponding to an eigenvalue is. Therefore, is does not belong to the spectral set σ (A△n ) of
the operator A△n . Thus, we have completed the proof of Lemma 3.2. □

Now, we are ready to present the uniform stability standard, as outlined in [13], which will be
instrumental in proving Theorem 3.2.

Theorem 3.1. Let h∗ > 0 and consider a family of semigroups of contractions (Sh (t)) on the
Hilbert space (X ˜︁h ). Let (A
˜︁h ) denote the corresponding infinitesimal generators. The family
(Sh (t)) is uniformly exponentially stable if and only if the following two conditions are met:
˜︁h ) of (A
(i) For all h ∈ (0, h∗ ), iR is contained in the resolvent set ρ(A ˜︁h ).
˜︁ −1
(ii) suph∈(0,h∗ ),β∈R ∥(iβI − Ah ) ∥L(X ˜︁ h ) < ∞.

By using the preceding results, we can now arrive at the main result of this section.

Theorem 3.2. Let h∗ = maxn≥1 △n . Then the semigroups T△n (t) generated by the operators A△n
are uniformly exponentially stable, i.e., there exist two positive constants M0 and ω0 independent
of △n and t such that

13
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

sup ∥T△n (t)∥△n ≤ M0 e−ω0 t .

n∈N +

Proof. Since (i) of Theorem 3.1 has already been established in Lemma 3.2, our focus in this
proof is to demonstrate that condition (ii) of Theorem 3.1 also holds. To do so, we employ a
proof by contradiction.
Suppose, to the contrary, that condition (ii) is false. Then, for any n ∈ N, there exist βn ∈ R,
mn ∈ N, and X△mn ∈ 𝒳△mn such that ∥X△mn ∥△mn = 1 and

∥Y△mn ∥△mn ≤ n−2 , (3.8)

where Y△mn = (iβn − A△mn )X△mn , with the property that mn → ∞ as n → ∞.

Firstly, the equation Y△mn = (iβn − A△mn )X△mn is equivalent to the following system
{︃
Cj [αj (iβn xj −1 − yj −1 ) + ˆ︁
αj (iβn xj − yj )] = xmn +j − xmn +j +1 ,
Lj [ˆ︁
αj (iβn xmn +j − ymn +j ) + αj (iβn xmn +j +1 − ymn +j +1 )] = xj −1 − xj ,

where j = 1, 2 · · · , 2mn , with the additional conditions x0 = y0 = 0, x2mn +1 = kxmn , and

y2mn +1 = kymn . By regrouping the terms in the above equations, we obtain
{︃
iβn αj Cj xj −1 − xmn +j = Cj (αj yj −1 + ˆ︁
αj yj ) − iβnˆ︁
αj Cj xj − xmn +j +1 ,
(3.9)
−xj −1 + iβnˆ︁ αj Lj xmn +j = Lj (ˆ︁
αj ymn +j + αj ymn +j +1 ) − iβn αj Lj xmn +j +1 − xj .

Since the determinant of the coefficient matrix

⃓ ⃓
⃓ iβn Cj αj −1 ⃓⃓
Dj := ⃓⃓ = −βn2 αj ˆ︁
αj Cj Lj − 1
−1 αj ⃓
iβn Lj ˆ︁

is always nonzero, we can uniquely solve the system (3.9) to obtain

⎧ α L [︁
iβ ˆ︁ ]︁
⎪
⎪ xj −1 = nDjj j Cj (αj yj −1 + ˆ︁
αj yj ) − iβnˆ︁
αj Cj xj − xmn +j +1
⎪
⎪ [︁ ]︁
⎨ αj ymn +j + αj ymn +j +1 ) − iβn αj Lj xmn +j +1 − xj ,
+ D1j Lj (ˆ︁
[︁ ]︁ (3.10)
⎪
⎪ xmn +j = D1j Cj (αj yj −1 + ˆ︁
αj yj ) − iβnˆ︁
αj Cj xj − xmn +j +1
⎪
⎪
⎩ iβ α C [︁ ]︁
+ nDjj j Lj (ˆ︁αj ymn +j + αj ymn +j +1 ) − iβn αj Lj xmn +j +1 − xj .

For any βn ∈ R, we derive the following relations: |Dj | ≥ 1, |Cj Lj αj ˆ︁

αj βn2 /Dj | < 1, and
⃓√︁ ⃓
⃓ ⃓
⃓ Cj [αj yj −1 + ˆ︁ αj yj ]⃓ = 𝒪(n−2 ), (3.11)
⃓√︁ ⃓
⃓ ⃓
αj ymn +j + αj ymn +j +1 ]⃓ = 𝒪(n−2 ),
⃓ Lj [ˆ︁ (3.12)
⃓ ⃓
⃓ Cj βn ⃓
Ij := ⃓
1 ⃓ ⃓ = 𝒪(1), (3.13)
Dj ⃓
⃓ ⃓
⃓ Lj βn ⃓
Ij2 := ⃓⃓ ⃓ = 𝒪(1). (3.14)
Dj ⃓

14
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

The relations (3.11) and (3.12) are straightforward consequences of ∥Y△mn ∥△mn = 𝒪(n−2 ). For
(3.13), we have
⃓ ⃓
⃓ Cj βn ⃓
Ij1 ⃓
=⃓ ⃓ ≤ √︁ Cj ,
αj Cj Lj βn ⃓
1 + αj ˆ︁ 2 αj ˆ︁
αj Cj Lj

where we utilized the inequality a 2 + b2 ≥ 2ab. By Assumptions 2.2-2.4, it follows that Ij1 =
𝒪(1). The inequality (3.14) can be proven similarly.
Secondly, (3.8) and (3.5) imply that

⃓ ⃓
k|xmn |2 = ⃓Re⟨(iβn − A△mn )X△mn , X△mn ⟩△mn ⃓
⃓ ⃓
= ⃓Re⟨Y△mn , X△mn ⟩△mn ⃓ ≤ ∥Y△mn ∥△mn = 𝒪(n−2 ).

Thus, by considering x2mn +1 = kxmn , we derive

{︃
|xmn |2 = 𝒪(n−2 ),
(3.15)
|x2mn +1 |2 = 𝒪(n−2 )

Substituting (3.15) into (3.10) with (3.10) with j = mn , we obtain

{︃
|xmn −1 |2 = 𝒪(n−2 ),
(3.16)
|x2mn |2 = 𝒪(n−2 ),

since the coefficients of Cj αj (yj −1 +ˆ︁ αj ymn +j + αj ymn +j +1 ), xmn and x2mn in (3.10)
αj yj ), Lj (ˆ︁
are all of the forms of 1/Dj , Ij1 , Ij2 and |Cj Lj αj ˆ︁ αj βn2 /Dj |, respectively. Similarly, for j =
1, 2, · · · , mn − 1, it can be readily proven by induction and using an analogous approach as from
(3.15) to (3.16) that

{︃
|xmn −j |2 = 𝒪(n−2 ),
(3.17)
|x2mn +1−j |2 = 𝒪(n−2 ).

Finally, by the definition of the norm ∥ · ∥△mn and the Assumptions 2.2 and 2.4, we have

∑︂
mn ∑︂
mn
∥X△mn ∥2△mn = Cj |αj xj −1 + ˆ︁
αj xj | + 2
Lj |ˆ︁
αj xmn +j + αj xmn +j +1 |2
j =1 j =1
mn (︂
∑︂ )︂
≤4C△mn |xj |2 + |xmn +j |2 + C△mn |x2mn +1 |2
j =1
−2
=𝒪(n ),

which contradicts the assumption that ∥X△mn ∥△mn = 1. □

15
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

4. Uniform exponential stability of Timoshenko beam

The Timoshenko beam equation, which describes the dynamics of a beam under transverse
loading and rotational effects, is given by
{︃
ρ ẅ(t, z) = K(w ′ (t, z) − ϕ(t, z))′ ,
(4.1)
Iρ ϕ̈(t, z) = EI ϕ ′′ (t, z) + K(w ′ (t, z) − ϕ(t, z)),

where

• w(t, z) represents the transverse displacement of the beam at time t and spatial position
z ∈ ℐ;
• ϕ(t, z) is the rotation angle of a filament of the beam at time t and spatial position z;
• ρ(z), Iρ (z), EI (z) and K(z) are the material properties of the beam: mass per unit length,
rotary moment of inertia of a cross-section, the product of Young’s modulus of elasticity and
the moment of inertia of a cross-section, and the shear modulus, respectively.

Here, the dot ˙ and the prime ′ denote derivatives with respect to time and spatial variables, respec­
tively. Additionally, we assume that the beam is clamped at the left-hand side (z = 0), meaning
that both the transverse displacement and the rotation angle are zero there. At the right-hand side
(z = S), we apply a damping force proportional to the velocity of the transverse displacement.
Therefore, the boundary conditions are

ẇ(t, 0) = 0, ϕ̇(t, 0) = 0, (4.2)

EI (S)ϕ ′ (t, S) = −k1 ϕ̇(t, S), (4.3)
′
K(S)(w (t, S) + ϕ(t, S)) = −k2 ẇ(t, S). (4.4)

4.1. PDE Timoshenko beam

To formulate the Timoshenko beam model as a port-Hamiltonian system, we introduce the

following physical notations that represent key dynamical quantities:

x1 (t, z) = wz (t, z) − ϕ(t, z),

which represents the shear displacement, capturing the difference between the transverse dis­
placement gradient and the rotation angle;

x2 (t, z) = ρwt (t, z),

which is the momentum associated with the transverse displacement of the beam;

x3 (t, z) = ϕz (t, z),

denoting the angular displacement gradient, or the rate of change of the rotation angle with
respect to the spatial position;

16
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

x4 (t, z) = Iρ ϕt (t, z),

which represents the angular momentum associated with the rotation of the beam’s filaments.
By utilizing equations (4.1) through (4.4), we derive the time derivatives of the variables
xi (t, z) for i = 1, 2, 3, 4, which are given by
⎧
⎪
⎪ ẋ1 (t, z) = [β2 x2 (t, z)]′ − β4 x4 (t, z),
⎪
⎪
⎪
⎪ ẋ2 (t, z) = [β1 x1 (t, z)]′ ,
⎨
ẋ3 (t, z) = [β4 x4 (t, z)]′ ,
(4.5)
⎪
⎪ ẋ4 (t, z) = [β3 x3 (t, z)]′ + β1 x1 (t, z),
⎪
⎪
⎪
⎪ x1 (t, S) = −𝒦2 x2 (t, S), x2 (t, 0) = 0,
⎩
x3 (t, S) = −𝒦1 x4 (t, S), x4 (t, 0) = 0,

where β1 = K, β2 = ρ −1 , β3 = EI and β4 = Iρ−1 , 𝒦1 = k1 β2 (S)/β1 (S) and 𝒦2 = k2 β4 (S)/

β3 (S).
Let us define the vector x(t, z) = (x1 (t, z), x2 (t, z), x3 (t, z), x4 (t, z))⊤ and introduce the di­
agonal matrix ℒ(z) = diag{β1 (z), β2 (z), β3 (z), β4 (z)}. Additionally, we define the matrices P1
and P0 as follows:
(︃ ∑︁ )︃ (︃ )︃
0 ∑︁ 0 1
P1 := 2 ∑︁ with 2 =
0 2 1 0

and
⎛ ⎞
0 0 0 −1
⎜0 0 0 0 ⎟
⎜
P0 := ⎝ ⎟.
0 0 0 0 ⎠
1 0 0 0

Utilizing these definitions and the time derivative expressions derived in (4.5), we can rewrite the
Timoshenko beam equations without boundary conditions into the form

ẋ(t, z) = P1 [ℒ(z)x(t, z)]′ + P0 [ℒ(z)x(t, z)]. (4.6)

The Hamiltonian/energy associated with either (4.6) or (4.5) is defined as

∫︂S
1
ℋ(t) = x ⊤ (t, z)ℒ(z)x(t, z)dz. (4.7)
2
0

According to Exercise 9.2 in [10], it is shown that if both k1 and k2 are positive, then the system
described by (4.5) is exponentially stable with respect to the energy H (t). For a more comprehen­
sive understanding of the well-posedness and exponential stability of (4.5), we refer to theorem
2.1 of [20].

17
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Theorem 4.1. Assuming βi (·) ∈ L∞ p (ℐ) for 1 ≤ i ≤ 4 and the initial value xi (0, z) := x0 ∈
i

L2 (ℐ), the system (4.5) admits a unique solution in the state space X := L2 (ℐ; R4 ) for the feed­
back gains k1 ≥ 0 and k2 ≥ 0. Furthermore, if k1 > 0 and k2 > 0, then there exist two positive
constants M1 and ω1 such that the energy defined by (4.7) along the solution to (4.5) satisfies

ℋ(t) ≤ M1 e−ω1 t ℋ(0). (4.8)

4.2. Semi-discretization of Timoshenko beam

In this subsection, we apply the concepts introduced in subsections 2.2 and 2.3 to discretize
the spatial variable of the Timoshenko beam in its port-Hamiltonian form (4.6). For this purpose,
we treat the components of the flow variables x(t, z) and the matrix L(z) as one-forms, whereas
the components e1 (t, z), e2 (t, z), e3 (t, z), e4 (t, z) of effort variables e(t, z) := ℒ(z)x(t, z) are
considered as zero-forms. The partitioning of the interval [0, S] into ℐ = ∪nj=1 ℐSj −1 Sj remains
valid in this section.
To approximate x1 (t, z) and x3 (t, z) on each interval [Sj −1 , Sj ], we utilize the one-form finite
q
elements ωSj −1 Sj . Specifically,

q
x1 (t, z) = X1,j (t)ωSj −1 Sj (z),
q
x3 (t, z) = X3,j (t)ωSj −1 Sj (z),

ϕ
where Ẋ1,j (t) = x1,j (t) and Ẋ3,j (t) = x3,j (t). The one-form finite elements ωSj −1 Sj are em­
ployed to approximate x2 (t, z) and x4 (t, z) within the interval [Sj −1 , Sj ]. Specifically, the ap­
proximations take the form

ϕ
x2 (t, z) = X2,j (t)ωSj −1 Sj (z),
ϕ
x4 (t, z) = X4,j (t)ωSj −1 Sj (z),

q
where Ẋ2,j (t) = x2,j (t) and Ẋ4,j (t) = x4,j (t). However, the zero-form finite elements ωSj −1 ,
q ϕ ϕ
ωSj , and ωSj −1 , ωSj are utilized to approximate the effort variables e1 (t, z), e3 (t, z), and e2 (t, z),
e4 (t, z) respectively, within the interval [Sj −1 , Sj ]. More precisely, we set:

q q
ei (t, z) = ei,j −1 (t)ωSj −1 (z) + ei,j (t)ωSj (z), i = 1, 3,
ϕ ϕ
ei (t, z) = ei,j −1 (t)ωSj −1 (z) + ei,j (t)ωSj (z), i = 2, 4,

where ei,j (t) will be determined later for i = 1, · · · , 4 and j = 0, 1, · · · , n. By utilizing the
concepts outlined in Remark 2.1, we derive the following relationships on the interval [Sj −1 , Sj ]:

ei,j (t) − ei,j −1 (t)

ei′ (t, z) ≈ , j = 1, 2, · · · , n,
hj
βi,j
ei (t, z) ≈ [αj ei,j −1 (t) + ˆ︁
αj ei,j (t)], i = 1, 3,
hj

18
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

βi,j
ei (t, z) ≈ [ˆ︁
αj ei,j −1 (t) + αj ei,j (t)], i = 2, 4,
hj

where hj = Sj − Sj −1 and
∫︂
βi,j = βi−1 (z)dz, i = 1, · · · , 4.
ℐSj −1 Sj

Consequently, we formulate the semi-discretization scheme for (4.5) as follows:

β1,j [αj ė1,j −1 (t) + ˆ︁

αj ė1,j (t)] = e2,j (t) − e2,j −1 (t) − hj [ˆ︁
αj e4,j −1 (t) + αj e4,j (t)], (4.9)
β2,j [ˆ︁
αj ė2,j −1 (t) + αj ė2,j (t)] = e1,j (t) − e1,j −1 (t), (4.10)
β3,j [αj ė3,j −1 (t) + ˆ︁
αj ė3,j (t)] = e4,j (t) − e4,j −1 (t), (4.11)
β4,j [ˆ︁
αj ė4,j −1 (t) + αj ė4,j (t)] = e3,j (t) − e3,j −1 (t) + hj [αj e1,j −1 (t) + ˆ︁
αj e1,j (t)], (4.12)
e1,n (t) = −𝒦1 e2,n (t), e2,0 (t) = 0, (4.13)
e3,n (t) = −𝒦2 e4,n (t), e4,0 (t) = 0, (4.14)

where j = 1, 2, · · · , n. Furthermore, the discrete energy ℋ△n is defined as follows:

∑︂
n
ℋ△n (t) = β1,j [αj e1,j −1 (t) + ˆ︁
αj e1,j (t)]2 + β2,j [ˆ︁
αj e2,j −1 (t) + αj e2,j (t)]2
j =1

∑︂
n
+ β3,j [αj e3,j −1 (t) + ˆ︁
αj e3,j (t)]2 + β4,j [ˆ︁
αj e4,j −1 (t) + αj e4,j (t)]2 . (4.15)
j =1

The state space associated with equations (4.9) through (4.14) resides in the Hilbert space, de­
noted as
⎧ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎫
⎪
⎪ Z1,△n zl,1 zk,0 ⎪
⎪
⎪
⎨ ⎜ zl,2 ⎟ ⎜ zk,1 ⎟ zl,j , zk,j ∈ C, ⎪
⎬
⎜ Z2,△ ⎟ ⎜ ⎟ ⎜ ⎟
X△n = Z△n = ⎝ ⎜ n⎟
∈ C 4n
: Z = ⎜ . ⎟ , Z = ⎜ . ⎟ , ,
⎪
⎪ Z3,△n ⎠ l,△ n
⎝ .. ⎠
k,△ n
⎝ .. ⎠ k = 1, 3, l = 2, 4, ⎪ ⎪
⎪
⎩ ⎪
⎭
Z4,△n zl,n zk,n−1
(4.16)

˜︁△n ⟩△n is defined as

where the inner product ⟨Z△n , Z

˜︁△n ⟩△n
⟨Z△n , Z
∑︂
n
= β1,j [αj z1,j −1 + ˆ︁
αj z1,j ][αj˜︁
z1,j −1 + ˆ︁
αj˜︁
z1,j ]
j =1

∑︂
n
+ β2,j [ˆ︁
αj z2,j −1 + αj z2,j ][ˆ︁
αj˜︁
z2,j −1 + αj˜︁
z2,j ]
j =1

19
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

∑︂
n
+ β3,j [αj z3,j −1 + ˆ︁
αj z3,j ][αj˜︁
z3,j −1 + ˆ︁
αj˜︁
z3,j ]
j =1

∑︂
n
+ β4,j [ˆ︁
αj z4,j −1 + αj z4,j ][ˆ︁
αj˜︁
z4,j −1 + αj˜︁
z4,j ],
j =1

˜︁△n ∈ X△n , subject to the additional constraints z2,0 = z4,0 = 0, z1,n = −𝒦1 z2,n ,
where Z△n , Z
and z3,n = −𝒦2 z4,n .
To express equations (4.9) through (4.14) in a vectorial form, we introduce several n × n
matrices. Specifically, let M△n and ℰ△n be defined as outlined in section 3. Furthermore, define
C△n = diag{h1 , h2 , · · · , hn },

⎛ ⎞ ⎛ ⎞
α1 α1 ˆ︁
α1
⎜ ˆ︁ ⎟ ⎜ .. ⎟
⎜ α2 α2 ⎟ ˆ︁ ⎜ α2 . ⎟
B△n =⎜ .. .. ⎟ , B△n = ⎜
⎜
⎟,
⎟
⎝ ⎠ ⎝ ..
. . αn−1 ⎠
. ˆ︁
ˆ︁
αn αn αn
ℒ△n = diag(L1 , L2 , L3 , L4 ), with Li = diag{βi,1 , βi,2 , · · · , βi,n }, i = 1, . . . , 4,
⎛ ⎞
ˆ︁△n −ˆ︁
B αn 𝒦1 ℰ△n 0 0
⎜ 0 B△n 0 0 ⎟
Φ△n = ⎜ ⎝ 0
⎟,
0 ˆ︁
B△n −ˆ︁ αn 𝒦2 ℰ△n ⎠
0 0 0 B△n
⎛ ⊤ ˆ︁⊤ ⎞
0 M△ n
0 −C△n B △n
⎜ −M△ 0 0 0 ⎟
Ψ △n = ⎝ ⎜ n ⎟,
0 0 0 M△ ⊤ ⎠
n
C△n B△n 0 −M△n 0
⎛ ⎞
0 0 0 0
⎜0 −𝒦1 ℰ△n 0 0 ⎟
Ω △n = ⎜ ⎝0
⎟.
⎠
0 0 0
0 −hnˆ︁ αn 𝒦1 ℰ△n 0 −𝒦2 ℰ△n

The state variables of equations (4.9) through (4.14) are consolidated into the vector Y△n (t),
defined as:
(︂ )︂⊤
⊤ ⊤ ⊤ ⊤
Y△n (t) = y1,△ n
(t), y 2,△ n
(t), y 3,△ n
(t), y 4,△ n
(t) ,

where each component vector yi,△n (t) is specified as follows

y1,△n (t) = (e1,0 (t), · · · , e1,n−1 (t))⊤ ,

y2,△n (t) = (e2,1 (t), · · · , e2,n (t))⊤ ,
y3,△n (t) = (e3,0 (t), · · · , e3,n−1 (t))⊤ ,
y4,△n (t) = (e4,1 (t), · · · , e4,n (t))⊤ .

20
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Subsequently, the system of equations (4.9) through (4.14) can be equivalently expressed in the
compact form
{︃
Ẏ△n (t) = 𝒜△n Y△n (t),
(4.17)
Y△n (0) = (yh1 , yh2 , yh3 , yh4 )⊤ ∈ X△n ,

where the matrix 𝒜△n = Φ−1 −1

△n ℒ△n (Ψ△n + Ω△n ), assuming that both Φ△n and ℒ△n are invertible,
as stated.
˜︁△n ∈ X△n , the inner product in the space X△n can be reformulated
Furthermore, for any Z△n , Z
as
⟨︁ ⟩︁
˜︁△n ⟩△n = Φ△n Z△n , ℒ△n Φ△n Z
⟨Z△n , Z ˜︁△n , (4.18)

where ⟨·, ·⟩ denotes the standard inner product in C 4n .

In the context of this section, Assumptions 2.2 and 2.3 continue to hold. However, Assump­
tion 2.4 is superseded by the following:

Assumption 4.1. Assume that βi (·) ∈ L∞p (ℐ) for 1 ≤ i ≤ 4. Furthermore, when combined with
Assumption 2.2, it implies that both max1≤j ≤n βi,j = 𝒪(n−1 ) and min1≤j ≤n βi,j = 𝒪(n−1 ) are
valid.

Furthermore, we require an additional assumption regarding hj , which is crucial for our anal­
ysis.

Assumption 4.2. Assume that hj satisfies the following inequality for all j = 1, 2, · · · , n:
[︃ ]︃
1 β1,j β2,j
hj ≤ + . (4.19)
2αj ˆ︁
αj β4,j β3,j

Given these assumptions, we can now present the following result on dissipativity.

Lemma 4.1. The operator 𝒜△n is dissipative on the space X△n for all n ∈ N + . Consequently,
𝒜△n generates a family of semigroups of contractions 𝒯△n (t) that ensures the existence and
uniqueness of the solution Y△n (t) to the system of equations (4.9)-(4.14). This solution satisfies
the balance equations in their discrete version, which can be expressed as

ℋ̇△n (t) = −𝒦1 |e2,n (t)|2 − 𝒦2 |e4,n (t)|2 . (4.20)

Proof. Given Z△n ∈ X△n defined by (4.16), we can derive the following identity by considering
the inner product involving the operator 𝒜△n :

⟨Z△n , 𝒜△n Z△n ⟩△n + ⟨𝒜△n Z△n , Z△n ⟩△n

⟨︁ ⟩︁ ⟨︁ ⟩︁
= Φ△n Z△n , (Ψ△n + Ω△n )Z△n + (Ψ△n + Ω△n )Z△n , Φ△n Z△n . (4.21)

From the definition of Φ△n , Ψ△n and Ω△n , we have

21
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

⟨︁ ⟩︁
Φ△n Z△n , (Ψ△n + Ω△n )Z△n

∑︂
n ∑︂
n
= β1,j [αj z1,j −1 + ˆ︁
αj z1,j ][z2,j − z2,j −1 ] + β2,j [ˆ︁
αj z2,j −1 + αj z2,j ][z1,j − z1,j −1 ]
j =1 j =1

∑︂
n ∑︂
n
+ β3,j [αj z3,j −1 + ˆ︁
αj z3,j ][z4,j − z4,j −1 ] + β4,j [ˆ︁
αj z4,j −1 + αj z4,j ][z3,j − z3,j −1 ]
j =1 j =1

∑︂
n
− hj [αj z1,j −1 + ˆ︁
αj z1,j ][ˆ︁
αj z4,j −1 + αj z4,j ]
j =1

∑︂
n
+ hj [ˆ︁
αj z4,j −1 + αj z4,j ][αj z1,j −1 + ˆ︁
αj z1,j ], (4.22)
j =1

and

⟨︁ ⟩︁
(Ψ△n + Ω△n )Z△n , Φ△n Z△n

∑︂
n ∑︂
n
= β1,j [z2,j − z2,j −1 ][αj z1,j −1 + ˆ︁
αj z1,j ] + β2,j [z1,j − z1,j −1 ][ˆ︁
αj z2,j −1 + αj z2,j ]
j =1 j =1

∑︂
n ∑︂
n
+ β3,j [z4,j − z4,j −1 ][αj z3,j −1 + ˆ︁
αj z3,j ] + β4,j [z3,j − z3,j −1 ][ˆ︁
αj z4,j −1 + αj z4,j ]
j =1 j =1

∑︂
n
− hj [ˆ︁
αj z4,j −1 + αj z4,j ][αj z1,j −1 + ˆ︁
αj z1,j ]
j =1

∑︂
n
+ hj [αj z1,j −1 + ˆ︁
αj z1,j ][ˆ︁
αj z4,j −1 + αj z4,j ]. (4.23)
j =1

By substituting equations (4.22) and (4.23) into (4.21), and leveraging the same calculation
methodology as in (3.5), we arrive at the following result

Re⟨𝒜△n Z△n , Z△n ⟩△n = −𝒦1 |z2,n |2 − 𝒦2 |z4,n |2 , (4.24)

where we have utilized the identities z1,n = −𝒦1 z2,n and z3,n = −𝒦2 z4,n . This confirms that the
balance equation (4.20) holds true, as can be seen by considering the energy functional ℋ△n (t) =
⟨Y△n (t), Y△n (t)⟩△n and utilizing equation (4.17). □

To provide a precise frequency domain analysis, we establish a crucial Lemma 4.2 following.

22
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Lemma 4.2. Define Γj by

⃓ ⃓
⃓ iαj β1,j β 1 0 αj ⃓⃓
hj ˆ︁
⃓
⃓ 1 iˆ︁
αj β2,j β 0 0 ⃓
Γj : = ⃓⃓ ⃓
⃓
⃓ 0 0 iαj β3,j β 1 ⃓
⃓ −hj αj 0 1 αj β4,j β ⃓
iˆ︁
= 1 + (β1,j β2,j + β3,j β4,j − h2j αj ˆ︁
αj β2,j β3,j )αj ˆ︁
αj β 2 + β1,j β2,j β3,j β4,j αj2ˆ︁
αj2 β 4 (4.25)

for any β ∈ R and 0 < hj ≤ Δn . Then,

(β1,j β2,j + β3,j β4,j )αj ˆ︁

αj β 2
Γj ≥ 1 + + β1,j β2,j β3,j β4,j αj2ˆ︁
αj2 β 4 ≥ 1. (4.26)
2

Furthermore, if we define Ii,j = (β∗,j )i β i Γ−1j for i = 0, 1, 2, 3, 4, where (β∗,j ) is the product
i

of i terms selected from {β1,j , β2,j , β3,j , β4,j }, then Ii,j are uniformly bounded, and the upper
bounds are all independent of β. Specifically, the following estimates hold:

|Ii,j | = 𝒪(1), i = 0, 1, 2, 3, 4. (4.27)

Proof. The inequality (4.26) holds due to Assumption 4.2, which implies that the middle term
of the right-hand side of (4.25) satisfies

1
(β1,j β2,j + β3,j β4,j − h2j αj ˆ︁ αj β 2 ≥ [β1,j β2,j + β3,j β4,j ]αj ˆ︁
αj β2,j β3,j )αj ˆ︁ αj β 2 .
2
This ensures the inequality |I0,j | ≤ 1. Furthermore, by (4.26), Assumption 4.1, and Young’s
inequality, we can derive bounds for |I1,j | for i = 1, 2, 3, 4. Specifically,

2β∗,j |β| β∗,j

|I1,j | ≤ ≤ √︁ = 𝒪(1),
2 + (β1,j β2,j + β3,j β4,j )αj ˆ︁
αj β 2 2(β1,j β2,j + β3,j β4,j )αj ˆ︁
αj
(β∗,j )2 β 2 1 (β∗,j )2
|I2,j | ≤ ≤ √︁ = 𝒪(1),
1 + β1,j β2,j β3,j β4,j αj2ˆ︁
αj2 β 4 2αj ˆ︁
αj β3,j β1,j β2,j β4,j

2(β∗,j )3 |β|
|I3,j | ≤
(β1,j β2,j + β3,j β4,j )αj ˆ︁
αj + 2β1,j β2,j β3,j β4,j αj2ˆ︁
αj2 β 2
(β∗,j )3
≤ √︁ = 𝒪(1),
αj ˆ︁
αj 2(β1,j β2,j + β3,j β4,j )β1,j β2,j β3,j β4,j αj ˆ︁
αj
2(β∗,j )4 β 4
|I4,j | ≤
2 + (β1,j β2,j + β3,j β4,j )αj ˆ︁
αj β 2 + 2β1,j β2,j β3,j β4,j αj2ˆ︁
αj2 β 4
(β∗,j )4
≤ = 𝒪(1),
β1,j β2,j β3,j β4,j αj2ˆ︁
αj2

which implies that (4.27) hold for i = 1, 2, 3, 4. □

23
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

The dissipativity of 𝒜△n guarantees that the spectral set σ (𝒜△n ) of 𝒜△n lies strictly within the
open left half-plane of the complex plane. This is a strengthening of the basic fact that σ (𝒜△n )
lies within the closed left half-plane, as dissipativity implies that the real part of every eigenvalue
of 𝒜△n is strictly negative, ensuring that no eigenvalues lie on the imaginary axis. Thus, for any
n ∈ N, the spectral set σ (𝒜△n ) is contained entirely within the open left half-plane of C.

Lemma 4.3. For every n ∈ N, iR ⊂ ρ(𝒜△n ).

Proof. If there exist β ∈ R and nonzero Z△n ∈ X△n such that iβZ△n = 𝒜△n Z△n , then it follows
from the definition of 𝒜△n that 𝒜△n Z△n = iβZ△n is equivalent to

iββ1,j [αj z1,j −1 + ˆ︁

αj z1,j ] = z2,j − z2,j −1 − hj [ˆ︁
αj z4,j −1 + αj z4,j ], (4.28)
iββ2,j [ˆ︁
αj z2,j −1 + αj z2,j ] = z1,j − z1,j −1 , (4.29)
iββ3,j [αj z3,j −1 + ˆ︁
αj z3,j ] = z4,j − z4,j −1 , (4.30)
iββ4,j [ˆ︁
αj z4,j −1 + αj z4,j ] = z3,j − z3,j −1 + hj [αj z1,j −1 + ˆ︁
αj z1,j ]. (4.31)

On the other hand, from (4.24), we derive the following:

⟨︁ ⟩︁ ⟨︁ ⟩︁
0 = Re iβZ△n , Z△n △ = Re 𝒜△n Z△n , Z△n △ = −𝒦1 |z2,n |2 − 𝒦2 |z4,n |2 .
n n

Given that z1,n = −𝒦1 z2,n and z3,n = −𝒦2 z4,n , the only way for the above equation to hold true
is if

z1,n = z2,n = z3,n = z4,n = 0. (4.32)

Substituting these values into (4.28)-(4.31) with j = n, we obtain

iββ1,n αn z1,n−1 + z2,n−1 + hnˆ︁

αn z4,n−1 = 0, (4.33)
z1,n−1 + iββ2,nˆ︁
αn z2,n−1 = 0, (4.34)
iββ3,n αn z3,n−1 + z4,n−1 = 0, (4.35)
−hn αn z1,n−1 + z3,n−1 + iββ4,nˆ︁
αn z4,n−1 = 0. (4.36)

The coefficients determinant of above equations is just Γn which is defined in Lemma 4.2 with
j = n. Since the Lemma 4.2 shows that the coefficients determinant Γj is positive, and by using
Cramer’s rule, Γj ̸= 0 is equivalent to the unique solution to the equations (4.34)-(4.35) are zeros,
i.e.,

z1,n−1 = z2,n−1 = z3,n−1 = z4,n−1 = 0. (4.37)

By induction and employing the same reasoning from (4.32) to (4.37), we can conclude that
zi,j = 0 for all i = 1, 2, 3, 4 and j = 0, 1, · · · , n. Consequently, Z△n = 0, which contradicts the
initial assumption that Z△n was nonzero. □

24
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Now, we are poised to present the main result of this section. Accordingly, we define h∗ as
{︃ [︃ ]︃}︃
1 β1,j β2,j
h∗ = max △n , max + .
n∈N 1≤j ≤n 2αj ˆ︁
αj β4,j β3,j

Theorem 4.2. For the matrices 𝒜△n defined by (4.17), the associated family of C0 -semigroups
𝒯△n (t) generated by 𝒜△n is uniformly exponentially stable. Specifically, there exist two constants
M > 0 and ω > 0, both independent of t and △n ∈ (0, h∗ ), such that for all △n in this range,

∥𝒯△n (t)∥ ≤ Me−ωt , ∀t ≥ 0. (4.38)

Proof. Building upon Lemma 4.1, we know that for every △n ∈ (0, h∗ ), the semigroup 𝒯△n (t)
is a C0 -semigroup of contractions. Moreover, Lemma 4.3 has already verified that 𝒜△n satisfies
the first condition of Theorem 3.1.
Now, following a similar contradiction argument as in the proof of Theorem 3.2, we suppose
that the second condition of Theorem 3.1 is false. If this supposition holds, then there exists a
sequence {(βk , △nk , Z△k )}
n k∈N + with nk + 1 = [1/△nk ] (where [a] denotes the largest integer
k
less than or equal to the real number a), βk ∈ R, △nk ∈ (0, h∗ ), and Z△
k
n
∈ X△nk such that
k

⎧ k
⎨ ∥Z△nk ∥X△nk = 1,
⎪
U△k n := (iβk I△nk − 𝒜△nk )Z△
k ,
(4.39)
⎪
⎩ ∥U kk ∥ −2
nk

△nk X△nk ≤k .

The proof is structured into the following three steps for clarity and precision.
Step 1: U△k n = (iβk I△nk − 𝒜△nk )Z△ k
n
is equivalent to
k k

β1,j [iβk (αj z1,j −1 + ˆ︁

αj z1,j ) − (αj u1,j −1 + ˆ︁
αj u1,j )] = z2,j − z2,j −1 − hj (ˆ︁
αj z4,j −1 + αj z4,j ),
(4.40)
β2,j [iβk (ˆ︁
αj z2,j −1 + αj z2,j ) − (ˆ︁
αj u2,j −1 + αj u2,j )] = z1,j − z1,j −1 , (4.41)
β3,j [iβk (αj z3,j −1 + ˆ︁
αj z3,j ) − (αj u3,j −1 + ˆ︁
αj u3,j )] = z4,j − z4,j −1 , (4.42)
β4,j [iβk (ˆ︁
αj z4,j −1 + αj z4,j ) − (ˆ︁
αj u4,j −1 + αj u4,j )] = z3,j − z3,j −1 + hj (αj z1,j −1 + ˆ︁
αj z1,j ),
(4.43)

for j = 1, · · · , nk . To streamline the above formulas, we introduce the simplifications z2,0 =

z4,0 = 0, z1,nk = −𝒦1 z2,nk , and z3,nk = −𝒦2 z4,nk for Z△k . Similarly, we apply these settings to
n k
U△k n .
k
Rearranging equations (4.40) through (4.43) results in

iβk β1,j αj z1,j −1 + z2,j −1 + hj ˆ︁

αj z4,j −1 =b1,j , (4.44)
z1,j −1 + iβk β2,j ˆ︁
αj z2,j −1 =b2,j , (4.45)
iβk β3,j αj z3,j −1 + z4,j −1 =b3,j , (4.46)
−hj αj z1,j −1 + z3,j −1 + iβk β4,j ˆ︁
αj z4,j −1 =b4,j , (4.47)

25
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

where

b1,j = β1,j (αj u1,j −1 + ˆ︁

αj u1,j ) + z2,j − iβk β1,j ˆ︁
αj z1,j − hj αj z4,j , (4.48)
b2,j = β2,j (ˆ︁
αj u2,j −1 + αj u2,j ) − iβk β2,j αj z2,j + z1,j , (4.49)
b3,j = β3,j (αj u3,j −1 + ˆ︁
αj u3,j ) − iβk β3,j ˆ︁
αj z3,j + z4,j , (4.50)
b4,j = β4,j (ˆ︁
αj u4,j −1 + αj u4,j ) + hj ˆ︁
αj z1,j + z3,j − iβk β4,j αj z4,j . (4.51)

It is straightforward to observe that the determinant of the coefficient matrix associated with
equations (4.44) through (4.47) is precisely Γj , as defined in Lemma 4.2, where β = βk . Ac­
cording to Lemma 4.2, Γj is guaranteed to be nonzero, implying that the system of equations
(4.44)-(4.47) possesses a unique solution for (z1,j −1 , z2,j −1 , z3,j −1 , z4,j −1 )⊤ , given by
⎛ ⎞
z1,j −1
⎜ z2,j −1 ⎟
⎜ ⎟
⎝ z3,j −1 ⎠
z4,j −1
⎛ ⎞
a2,j Bj Bj hj ˆ︁
αj a2,j −hj ˆ︁
αj a2,j a3,j
⎜ −Bj a1,j Bj + ˆ︁
αj αj h2j a3,j −hj ˆ︁
αj αj a3,j ⎟
hj ˆ︁
= Γ−1 ⎜ ⎟
j ⎝ −h α a hj αj a4,j Aj + αj ˆ︁
αj h2j a2,j −Aj ⎠
j j 2,j
−hj ˆ︁
αj a2,j a3,j −hj αj a3,j −Aj a3j Aj
⎛ ⎞
b1,j
⎜ b2,j ⎟
×⎜ ⎟
⎝ b3,j ⎠ (4.52)
b4,j

where a1,j = iβk β1,j αj , a2,j = iβk β2,j ˆ︁

αj , a3,j = iβk β3,j αj , a4,j = iβk β4,j ˆ︁
αj , Aj = a1,j a2,j −
1, and Bj = a3,j a4,j − 1.
Step 2: It is easy to see that ∥U△k n ∥X△n ≤ k −2 directly implies that
k k

αj u1,j )| = 𝒪(k −2 ), β2,j |(ˆ︁

β1,j |(αj u1,j −1 + ˆ︁ αj u2,j −1 + αj u2,j )| = 𝒪(k −2 ), (4.53)
αj u3,j )| = 𝒪(k −2 ), β4,j |(ˆ︁
β3,j |(αj u3,j −1 + ˆ︁ αj u4,j −1 + αj u4,j )| = 𝒪(k −2 ), (4.54)

for all j = 1, 2, · · · , nk . From (4.39) and (4.24), we deduce that

𝒦1 |z2,nk |2 + 𝒦2 |z4,nk |2 = Re⟨(iβk − 𝒜△nk )Z△

k
n
k
, Z△ ⟩
n X△n
= Re⟨U△k n , Z△
k
⟩
n X△n
≤ k −2 .
k k k k k k

Combining this with the given conditions z1,nk = −𝒦1 z2,nk and z3,nk = −𝒦2 z4,nk , we obtain

zi,nk = 𝒪(k −1 ), i = 1, 2, 3, 4. (4.55)

Step 3: We assert the claim that

zi,j = 𝒪(k −1 ), i = 1, 2, 3, 4, j = 1, 2, · · · , nk . (4.56)

26
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

In fact, from (4.48)-(4.51) and (4.52) with j = nk , we derive

| z1,nk −1 |
⃓
= | Γ−1 ⃓
nk | a2,nk Bnk (β1,nk (αnk u1,nk −1 + ˆ︁
αnk u1,nk ) + z2,nk − iβk β1,nk ˆ︁
αnk z1,nk − hnk αnk z4,nk )
+ Bnk (β2,nk (ˆ︁
αnk u2,nk −1 + αnk u2,nk ) − iβk β2,nk αnk z2,nk + z1,nk )
+ hnk ˆ︁
αnk a2,nk (β3,nk (αnk u3,nk −1 + ˆ︁
αnk u3,nk ) − iβk β3,nk ˆ︁
αj z3,nk + z4,nk )
⃓
−hnk ˆ︁ αnk z1,nk + z3,nk − iβk β4,nk αnk z4,nk )⃓
αnk u4,nk −1 + αnk u4,nk ) + hnk ˆ︁
αnk a2,nk a3,nk (β4,nk (ˆ︁
[︁
≤ | Γ−1
nk | |a2,nk Bnk |(β1,nk |(αnk u1,nk −1 + ˆ︁
αnk u1,nk )| + |z2,nk | + β1,nk |βk ||z1,nk | + |z4,nk |)
+ |Bnk |(β2,nk |ˆ︁
αnk u2,nk −1 + αnk u2,nk | + β2,nk |βk ∥z2,nk | + |z1,nk |)
+ |a2,nk |(β3,nk |αnk u3,nk −1 + ˆ︁
αnk u3,nk | + β3,nk |βk ∥z3,nk | + |z4,nk |)
]︁
+|a2,nk a3,nk |(β4,nk |ˆ︁
αnk u4,nk −1 + αnk u4,nk | + |z1,nk | + |z3,nk | + β4,nk |βk ∥z4,nk |) .

The coefficients of β1,nk |(αnk u1,nk −1 + ˆ︁ αnk u1,nk )|, β2,nk |ˆ︁
αnk u2,nk −1 + αnk u2,nk |,
β3,nk |αnk u3,nk −1 + ˆ︁
αnk u3,nk |, β4,nk |ˆ︁
αnk u4,nk −1 + αnk u4,nk | along with |zi,nk |(i = 1, 2, 3, 4) all
share the form of |Ii,nk | as defined in Lemma 4.2 with j = nk and β = βk . According to
Lemma 4.2 and equations (4.53)-(4.55), it follows that z1,nk −1 = 𝒪(k −1 ). Similarly, we can
deduce that |zi,nk −1 | = 𝒪(k −1 ) for i = 2, 3, 4. This establishes that (4.55) implies (4.56) for
j = nk − 1 with the aid of (4.52)-(4.54). By extending this reasoning through induction, we can
prove that (4.56) holds for all j = 1, 2, · · · , nk .
Finally, utilizing (4.56) alongside Assumptions 2.2 and 4.1, we obtain

∑︂
nk ∑︂
nk
∥ Z△
k
n
∥2X△ = β1,j |αj z1,j −1 + ˆ︁
αj z1,j |2 + β2,j |ˆ︁
αj z2,j −1 + αj z2,j |2
k nk
j =1 j =1

∑︂
nk ∑︂
nk
+ β3,j |αj z3,j −1 + ˆ︁
αj z3,j |2 + β4,j |ˆ︁
αj z4,j −1 + αj z4,j |2
j =1 j =1
4 ∑︂
∑︂ nk
≤2C△nk |zi,j |2 = 𝒪(k −2 ),
i=1 j =0

which contradicts the fact that ∥ Z△

k
n
∥X△n = 1. This completes the proof of the theorem. □
k k

5. Numerical simulations

In this section, we present numerical simulations to demonstrate the validity of our theoretical
findings. These simulations are performed under the assumption of uniform mesh size hj = h =
S/n, where S = 1. The coefficients Cj and Lj in (2.20) are approximated as

C(j h)h ≈ Cj , L(j h)h ≈ Lj , j = 0, · · · , n.

We present four figures to highlight the significance of the discrete schemes (2.20) and (3.4).
In Fig. 1, the blue points represent the maximal real parts of the eigenvalues of A△n obtained

27
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Fig. 1. Maximal real parts of eigenvalues of A△N with αj = 3/4. (For interpretation of the colors in the figure(s), the
reader is referred to the web version of this article.)

Fig. 2. Eigenvalues distribution of A△n with n = 100 and αj = 3/4.

from (3.4) with αj = 3/4 and the inclusion of the weighting operator Φ△n . In contrast, the red
points represent the maximal real parts of the eigenvalues of A△n using the same αj but without
the weighting operator Φ△n , which essentially reduces to the classical finite difference scheme.
Notably, the maximal real parts of the red points approach zero, indicating that the classical finite
difference scheme fails to uniformly preserve the exponential stability of (2.1). This observation
aligns with the conclusions presented in [9]. However, for the semidiscrete scheme (2.20) with
the same step size, the maximal real parts of the eigenvalues approach a negative number, which
is consistent with the statement of Theorem 3.1.
With αj = 3/4, Fig. 2 illustrates the distribution of the eigenvalues of A△n both with and
without the weighting operator Φ△n . This figure reinforces the same conclusions drawn from
Fig. 1. Similarly, Figs. 3 and 4 mirror Figs. 1 and 2, respectively, but with αj = 1/2. For αj =
1/2, similar numerical simulations were previously reported in [21]. In these figures, we use
C(z) = ln(1 + z), L(z) = exp(z), R = 5.

28
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Fig. 3. Maximal real parts of eigenvalues of A△N with αj = 1/2.

Fig. 4. Eigenvalues distribution of A△n with n = 100 and αj = 1/2.

Furthermore, numerical experiments indicate that the discrete schemes (2.20) and (3.4) with
αj = 1/2 and identical parameters exhibit an optimal decay rate. Specifically, the maximal real
parts observed in Figs. 3 and 4 are notably smaller than those in Figs. 1 and 2, respectively.
These findings are consistent across various numerical simulations conducted with αj ̸= 1/2.
However, we are only drawing conclusions from numerical simulation results, which deserve
rigorous theoretical verification in the future.
For the Timoshenko beam model, we have generated four additional figures, numbered Fig. 5
through Fig. 8, to showcase the numerical results obtained from the discrete schemes (4.9)-(4.12)
or (4.17). For these simulations, we set 𝒦1 = 𝒦2 = 1, and choose β1 = exp(z), β2 = 1 + 2z,
β3 = 2 + sin z, and β4 = 1 + z2 . Analogous to the previous cases, βi,j is approximated as

∫︂Sj
βi,j = βi−1 (z)dz ≈ hβi−1 (j h),
Sj −1

29
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Fig. 5. Maximal real parts of eigenvalues of 𝒜△N with αj = 3/4.

Fig. 6. Eigenvalues distribution of 𝒜△n with n = 100 and αj = 3/4.

where Sj = j h for i = 1, 2, 3, 4 and j = 0, · · · , n. Figs. 5 through 8 demonstrate the same effec­

tiveness and conclusions as Figs. 1 through 4.

6. Concluding remarks and further researches

In this paper, a generalization of spatially discretization scheme by mixed finite element

method is firstly proposed for the transmission line with a boundary resistor. The resulting fi­
nite dimensional approximation scheme preserves exponential stability of the original infinite
dimensional system. Notably, the main contribution of this paper is verifying the uniform expo­
nential stability through the frequency domain characterization by Liu and Zheng in [13], which
gives an equivalent relation between this property and the uniform order of the corresponding re­
solvent operator on the imaginary axis. To estimate this order, a contradiction argument is used.
Secondly, the proposed spatially discretization scheme is easily translated to the Timoshenko
beam with boundary damping. The same result is derived by the same method though the fre­

30
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Fig. 7. Maximal real parts of eigenvalues of 𝒜△N with αj = 1/2.

Fig. 8. Eigenvalues distribution of 𝒜△n with n = 100 and αj = 1/2.

quency domain characterization for the continuous model is absent. Our results can be potentially
applied to the related LQR problem and numerical approximating of state reconstruction in the
further research. Moreover, it also deserves to studying the preservations of some control prop­
erties of more complex infinite dimensional systems using the spatially discretization scheme of
this paper and the idea of Remark 3.1.

Acknowledgments

The authors wish to express their sincere gratitude to the anonymous referees for their meticu­
lous review of the manuscript, along with their insightful comments and constructive suggestions
that have significantly enhanced its quality. Particular appreciation is given to one reviewer for
providing an alternative proof, which is detailed in Remark 3.1.

31
F. Zheng, H. Yin, Z. Han et al. Journal of Differential Equations 453 (2026) 113865

Data availability

Data will be made available on request.

References

[1] H.T. Banks, K. Ito, C. Wang, Exponentially Stable Approximations of Weakly Damped Wave Equations, Internat.
Ser. Numer. Math., vol. 100, Birkhäuser Verlag, Basel, 1991, pp. 1--33.
[2] H.T. Banks, C. Wang, Optimal feedback control of infinite-dimensional parabolic evolution systems: approximation
techniques, SIAM J. Control Optim. 27 (1989) 1182--1219.
[3] A. Brugnoli, R. Rashad, S. Stramigioli, Dual field structure-preserving discretization of port-Hamiltonian systems
using finite element exterior calculus, J. Comput. Phys. 471 (2022) 111601.
[4] J.S. Gibson, Linear-quadratic optimal control of hereditary differential systems: infinite dimensional Riccati equa­
tions and numerical approximations, SIAM J. Control Optim. 21 (1983) 95--139.
[5] J.S. Gibson, I.G. Rosen, G. Tao, Approximation in control of thermoelastic systems, SIAM J. Control Optim. 30
(1992) 1163--1189.
[6] G. Golo, V. Talasila, A.J. van der Schaft, B. Maschke, Hamiltonian discretization of boundary control systems,
Automatica 40 (2004) 757--771.
[7] B.Z. Guo, F. Zheng, Uniform exponential stability for a Schrödinger equation and its semi-discrete approximation,
IEEE Trans. Autom. Control 69 (2024) 8900--8907.
[8] C. Harkort, J. Deutscher, Stability and passivity preserving Petrov-Galerkin approximation of linear infinite­
dimensional systems, Automatica 48 (2012) 1347--1352.
[9] J.A. Infante, E. Zuazua, Boundary observability for the space semi-discretizations of the 1-d wave equation, M2AN
Math. Model. Numer. Anal. 33 (1999) 407--438.
[10] B. Jacob, H. Zwart, Linear Port-Hamiltonian System on Infinite-Dimensional Space, Springer, Basel, 2012.
[11] L. León, E. Zuazua, Boundary controllability of the finite-difference space semi-discretizations of the beam equa­
tion, ESAIM Control Optim. Calc. Var. 8 (2002) 827--862.
[12] J. Liu, B.Z. Guo, A new semi-discretized order reduction finite difference scheme for uniform approximation of
1-D wave equation, SIAM J. Control Optim. 58 (2020) 2256--2287.
[13] Z.Y. Liu, S.M. Zheng, Uniform exponential stability and approximation in control of a thermoelastic system, SIAM
J. Control Optim. 32 (1994) 1226--1246.
[14] A. Macchelli, Energy shaping of distributed parameter port-Hamiltonian systems based on finite element approxi­
mation, Syst. Control Lett. 60 (2011) 579--589.
[15] S. Micu, C. Castro, Boundary controllability of a linear semi-discrete 1-D wave equation derived from a mixed finite
element method, Numer. Math. 102 (2006) 413--462.
[16] K. Ramdani, T. Takahashi, M. Tucsnak, Uniformly exponentially stable approximations for a class of second order
evolution equations-application to LQR problems, ESAIM Control Optim. Calc. Var. 13 (2007) 503--527.
[17] H.J. Ren, B.Z. Guo, Uniform exponential stability of semi-discrete scheme for observer-based control of 1-D wave
equation, Syst. Control Lett. 168 (2022) 105346.
[18] L.T. Tebou, E. Zuazua, Uniform boundary stabilization of the finite difference space discretization of the 1-d wave
equation, Adv. Comput. Math. 26 (2007) 337--365.
[19] V. Trenchant, H. Ramirez, Y. Le Gorrec, P. Kotyczka, Finite differences on staggered grids preserving the port­
Hamiltonian structure with application to an acoustic duct, J. Comput. Phys. 373 (2018) 673--697.
[20] X.F. Wang, W.L. Xue, Y. He, F. Zheng, Uniformly exponentially stable approximations for Timoshenko beams,
Appl. Math. Comput. 451 (2023) 128028.
[21] B.F. Zhang, F. Zheng, Y. He, Uniformly exponentially stable approximation for the transmission line with variable
coefficients and its application, J. Appl. Anal. Comput. 14 (2024) 2228--2256.
[22] F. Zheng, S. Zhang, H. Wang, B.Z. Guo, The exponential stabilization of a heat-wave coupled system and its
approximation, J. Math. Anal. Appl. 521 (2023) 126927.
[23] F. Zheng, H. Zhou, State reconstruction of the wave equation with general viscosity and non-collocated observation
and control, J. Math. Anal. Appl. 502 (2020) 125257.
[24] E. Zuazua, Propagation, observation, and control of waves approximated by finite difference methods, SIAM Rev.
47 (2005) 197--243.

32