Kunjie Ji
School of Mathematics and Statistics, Shandong Normal University, Jinan, China.
DOI: 10.4236/eng.2026.189019   PDF    HTML   XML   6 Downloads   64 Views   Citations

Abstract

This paper considers the joint identification of the wave speed and non-generic initial states of a one-dimensional wave equation from finite-time boundary measurements. An on-off boundary input partitions the observation horizon into zero-input and controlled phases, whose complementary information guarantees the unique identifiability of both the unknown wave speed and the initial states. By exploiting the exponential-series representation of the boundary response, the inverse problem is reformulated as the recovery of spectral data and modal coefficients from a perturbed exponential sequence. A multi-stage identification framework is then constructed by combining matrix-pencil spectral estimation, controlled-residual fitting, and regularized least-squares reconstruction. Perturbation bounds are established for the recovered discrete poles and characteristic exponents, quantifying the effect of modal truncation on the spectral estimates. Numerical experiments demonstrate accurate identification of the wave speed and reconstruction of the initial displacement and velocity. Monte Carlo simulations under multiplicative measurement noise further confirm the robustness of the proposed framework.

Keywords

Wave Equation, Inverse Problem, Simultaneous Identification, Boundary Switching Control, Matrix Pencil Method

Share and Cite:

1. Introduction

In many dynamical systems, the identification of unknown parameters is essential for understanding system behavior and ensuring effective control [1]-[8]. This process becomes even more critical in distributed parameter systems, where parameters vary across space and time, and accurate identification is the key to reliable modeling and monitoring [9]-[12]. In such systems, especially for wave equations, intrinsic physical parameters (e.g., wave speed), initial states, and internal disturbances are generally not directly observable and must be inferred from indirect measurements [13]. Therefore, reconstructing these unknown parameters through boundary measurements has received extensive attention in inverse problems associated with wave-type distributed parameter systems. For hyperbolic systems, particularly wave equations, inverse problems play a crucial role in understanding wave propagation dynamics and enabling effective monitoring and control [14]. Typical inverse problems for wave equations include the identification of wave speed, source terms, initial conditions, and boundary disturbances from partial observations [2]-[4] [8]-[10].

Early studies in this field mainly focused on inverse problems involving a single unknown for wave systems. Over the past few decades, extensive investigations have been conducted on inverse problems for one-dimensional wave equations from multiple perspectives, including the estimation of medium parameters such as wave speed [15] [16], dispersion coefficients [17], and disturbance parameters [18] [19], as well as the reconstruction of initial state distributions [20] and spatial source terms [16] [21]. Provided that suitable observability conditions are satisfied, these works have obtained essential results concerning uniqueness and stability, while various analytical and computational approaches have been proposed, such as the boundary control method [22], spectral methods [23], and adaptive identification strategies [16] [24]. For a more comprehensive discussion of numerical techniques for inverse problems arising from partial differential equations, interested readers are referred to the monographs [25] [26].

In recent years, research attention has gradually shifted toward inverse problems involving multiple unknowns, also known as hybrid or coupled inverse problems. For wave equations, numerous works have investigated the simultaneous identification of multiple unknown quantities, either of the same type or different types, such as multiple boundary disturbances [27], two distinct source terms [28], combinations of wave speed and source terms [29], as well as measurement biases and initial states [30]. Nevertheless, most existing results rely on additional internal measurements or multiple boundary observations, which may be infeasible in practical engineering applications [19] [27]. Furthermore, the simultaneous reconstruction of the wave speed together with two initial states has been rarely addressed, especially under the constraint of a single boundary observation. The difficulties outlined above, together with the recent progress reported in [31] for inverse heat-conduction models, in which the diffusion coefficient and the initial state are recovered within the same framework, provide the main motivation for the present study. Our objective is to determine the wave speed and both unknown initial data in a unified manner by using measurements acquired at only one boundary.

This study investigates an inverse problem for wave propagation in a homogeneous bar of unit length. The unknown wave velocity, initial displacement, and initial velocity are uniquely identifiable simultaneously from the available measurement data. The dynamics of the system are governed by

$\{\begin{array}{ll}{u}_{tt}\left(x,t\right)=c^2{u}_{xx}\left(x,t\right),x\in \left(0,1\right),\hfill & t>0,\hfill \\ u_x\left(0,t\right)=U\left(t\right),\hfill & t\ge 0,\hfill \\ u_x\left(1,t\right)=-qu\left(1,t\right),\hfill & t\ge 0,q>0,\hfill \\ y\left(t\right)=u\left(0,t\right),\hfill & t\ge 0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),u_t\left(x,0\right)=u_1\left(x\right),\hfill & x\in \left[0,1\right].\hfill \end{array}$ (1.1)

Here, $x$ and $t$ are the spatial and temporal variables, respectively. The propagation velocity $c$ is an unknown positive constant satisfying $c\ge c_0>0$ , while $q>0$ denotes the elastic coefficient at the boundary. The two unknown initial profiles are given by $u_0\in H^1\left(0,1\right)$ and $u_1\in L^2\left(0,1\right)$ , corresponding to the displacement and velocity at $t=0$ , respectively. Apart from boundedness, no additional restriction is imposed on either initial datum. The signal $U\left(t\right)$ serves as a Neumann-type boundary input and describes the stress flux applied at the left endpoint of the bar. The measured output $y\left(t\right)$ is the displacement recorded at the boundary. To explicitly indicate the dependence of the solution on the control input and the initial data, the solution of (1.1) may be denoted by $u=u\left(x,t;U,u_0,u_1\right)$ .

Inverse Problem. Consider system (1.1), where the wave velocity $c$ and the initial data ${\left(u_0,u_1\right)}^{\top }$ are unknown. The objective is to construct a boundary input $U\in L^2\left(0,T\right)$ such that the observation $y\left(t\right)=u\left(0,t\right)$ , $t\in \left(0,T\right)$ , uniquely determines $c$ , $u_0\left(x\right)$ , and $u_1\left(x\right)$ .

The inverse problem considered here has two notable features. First, the wave velocity and both initial profiles must be determined simultaneously, whereas the available information consists solely of the displacement trace recorded at the left endpoint $x=0$ . Recovering several unknown quantities from such limited data is therefore particularly difficult. Second, an appropriately chosen boundary excitation $U\left(t\right)$ is essential for distinguishing the wave velocity from the two initial states. Its purpose is to enrich the system dynamics so that the resulting output $y\left(t\right)$ carries enough information to identify all three unknown quantities uniquely. In this respect, the boundary input serves a role analogous to that of persistent excitation in adaptive control. To demonstrate this fact, consider the following two different collections of unknown quantities, where ${\omega }_k$ and ${\omega }_m$ ($k\ne m$ ) denote two distinct positive roots of the characteristic equation $\omega \tan \omega =q$ of the Robin boundary-value problem. The corresponding eigenfunctions $\text{cos}\left({\omega }_{n}x\right)$ satisfy both the Neumann condition at $x=0$ and the Robin condition at $x=1$ :

• Case 1: Let $c_1>0$ be arbitrary, $u_{01}\left(x\right)=\text{cos}\left({\omega }_{k}x\right)$ , and $u_{11}\left(x\right)=0$ .

The zero-input solution takes the form $u_1\left(x,t\right)=\text{cos}\left({\omega }_k{c}_{1}t\right)\text{cos}\left({\omega }_{k}x\right)$ , which yields the boundary output $y\left(t\right)=\text{cos}\left({\omega }_k{c}_{1}t\right)$ at $x=0$ .

• Case 2: Let $c_2=\frac{{\omega }_k}{{\omega }_m}{c}_1$ , $u_{02}\left(x\right)=\cos \left({\omega }_{m}x\right)$ , and $u_{12}\left(x\right)=0$ .

The zero-input solution becomes

$u_2\left(x,t\right)=\text{cos}\left({\omega }_m{c}_{2}t\right)\text{cos}\left({\omega }_{m}x\right)=\text{cos}\left({\omega }_k{c}_{1}t\right)\text{cos}\left({\omega }_{m}x\right)$ , which produces the

identical boundary output $y\left(t\right)=\text{cos}\left({\omega }_k{c}_{1}t\right)$ at $x=0$ .

By construction, both initial profiles are eigenfunctions of the underlying Sturm-Liouville problem and thus satisfy the boundary conditions of system (1.1). Nevertheless, the two distinct sets of unknown parameters yield exactly the same boundary displacement trace. Hence, the measured trace cannot be used to discriminate between the two cases. Equivalently, when no boundary excitation is applied, the input-output map is not injective. Therefore, the observation $y\left(t\right)$ alone is insufficient to uniquely determine the wave velocity $c$ together with the initial data $u_0$ and $u_1$ .

The contributions of this work can be summarized in three points. First, motivated by the switched boundary-input strategy proposed in [31] for the joint identification of the diffusion coefficient and the initial profile of a one-dimensional heat equation, we adapt the underlying identification idea to a hyperbolic wave equation. Second, in view of the oscillatory nature of wave propagation, which is fundamentally different from the dissipative behavior of heat conduction, we establish a reconstruction framework for the simultaneous identification of the constant wave speed $c$ and two initial states $u_0\left(x\right)$ and $u_1\left(x\right)$ from a single boundary observation. Third, in contrast to many identifiability results based on inverse spectral theory, which assume that the initial data have a nonzero projection onto every eigenfunction, our method eliminates this requirement through the designed boundary control and therefore permits some modal coefficients of the initial states to vanish.

The remainder of this article is arranged as follows. Section 2 studies the joint identifiability of the wave speed $c$ and the two initial states by means of a Fourier-series argument. Section 3 develops an identification procedure that integrates the matrix pencil technique, controlled-residual fitting, and regularized least-squares reconstruction. Section 4 examines the errors associated with applying the matrix pencil technique to the corresponding infinite-dimensional spectral estimation problem. Finally, Section 5 reports numerical experiments that assess the performance of the procedure developed in Section 3.

2. Identifiability

System (1.1) is formulated on the Hilbert space $\mathscr{H}=H^1\left(0,1\right)\times L^2\left(0,1\right)$ . We are equipped $\mathscr{H}$ with the following inner product ${\langle \cdot ,\cdot \rangle }_{\mathscr{H}}$ defined by

${\langle \left(f_1,g_1\right),\left(f_2,g_2\right)\rangle }_{\mathscr{H}}={\displaystyle {\int }_0^1{c}^2{f^{\prime }}_1\overline{{f^{\prime }}_2}\text{d}x}+{\displaystyle {\int }_0^1{g}_1\overline{g_2}\text{d}x}+c^{2}q{f}_2\left(1\right)\overline{f_1\left(1\right)}$ (2.1)

and the norm induced by this inner product will be written as $\Vert \cdot \Vert$ .

Let the operator $\mathcal{A}:D\left(\mathcal{A}\right)\left(\subset \mathscr{H}\right)\to \mathscr{H}$ be specified by

$\mathcal{A}\left(\begin{array}{l}f\hfill \\ g\hfill \end{array}\right)=\left(\begin{array}{c}g\\ c^2{f}^{″}\end{array}\right),$ (2.2)

where its domain is given by

$D\left(\mathcal{A}\right)=\left\{{\left(f,g\right)}^{\top }\in H^2\left(0,1\right)\times H^1\left(0,1\right)|f^{\prime }\left(0\right)=0,f^{\prime }\left(1\right)=-qf\left(1\right),q>0\right\}.$ (2.3)

With respect to the inner product defined in (2.1), the operator $\mathcal{A}$ is skew-adjoint on the separable Hilbert space $\mathscr{H}$ . It follows that all its eigenvalues lie on the imaginary axis and appear in complex-conjugate pairs; namely,

${\lambda }_n^{\pm }=\pm \text{i}{\omega }_{n}c,n=1,2,3,\cdots ,$ (2.4)

Here, the sequence ${\left\{{\omega }_n\right\}}_{n=1}^{\infty }$ consists of all positive solutions to

$\omega \tan \omega =q,q>0.$ (2.5)

For each eigenvalue ${\lambda }_n^{\pm }$ , an associated eigenvector can be chosen as

${\varphi }_n^{\pm }\left(x\right)={\left(\begin{array}{l}\cos \left({\omega }_{n}x\right),\pm i{\omega }_{n}c\cos \left({\omega }_{n}x\right)\hfill \end{array}\right)}^{\top },n=1,2,3,\cdots .$ (2.6)

It is not immediate from the unboundedness and skew-adjointness of $\mathcal{A}$ that its eigenvectors form an orthogonal basis for $\mathscr{H}$ . This property can instead be established through the inverse operator. Indeed, the Rellich-Kondrachov compact embedding theorem ensures that ${\mathcal{A}}^{-1}$ is compact. Moreover, the skew-adjointness of $\mathcal{A}$ yields ${\left({\mathcal{A}}^{-1}\right)}^{\text{*}}=-{\mathcal{A}}^{-1}$ , and hence ${\mathcal{A}}^{-1}$ is also skew-adjoint. Since $\mathscr{H}$ is separable, the spectral theorem for compact normal operators guarantees that the eigenvectors associated with ${\mathcal{A}}^{-1}$ constitute a complete orthogonal system in $\mathscr{H}$ . Furthermore, $\mathcal{A}$ and ${\mathcal{A}}^{-1}$ possess the same eigenvectors, with their corresponding eigenvalues being reciprocal. It therefore follows that ${\left\{{\varphi }_n^{\pm }\left(x\right)\right\}}_{n=1}^{\infty }$ forms an orthogonal basis of $\mathscr{H}$ .

Normalizing these eigenvectors yields the orthonormal basis ${\left\{{\tilde{\varphi }}_n^{\pm }\left(x\right)\right\}}_{n=1}^{\infty }$ , i.e.,

${\tilde{\varphi }}_n^{\pm }\left(x\right)=\frac{1}{{\alpha }_n}{\varphi }_n^{\pm }\left(x\right),$ (2.7)

where the normalization constant is

${\alpha }_n=c\sqrt{{\omega }_n\left({\omega }_n+\frac{1}{2}\text{sin}\left(2{\omega }_n\right)\right)},$ (2.8)

satisfying ${\Vert {\tilde{\varphi }}_n^{\pm }\Vert }_{\mathscr{H}}=1$ .

Define the state vector

$z\left(t\right)={\left(u\left(t\right),u_t\left(t\right)\right)}^{\top },$

where $u\left(t\right)=u\left(\cdot ,t\right)$ and $u_t\left(t\right)=\frac{\partial }{\partial t}u\left(\cdot ,t\right)$ , and the initial state

$z_0={\left(u_0,u_1\right)}^{\top }$ (2.9)

encodes the initial displacement and velocity prescribed by $u\left(x,0\right)=u_0\left(x\right)$ and $u_t\left(x,0\right)=u_1\left(x\right)$ .

Accordingly, system (1.1) admits the following first-order representation in the dual space ${\left[D\left(\mathcal{A}\right)\right]}^{\prime }$ :

$\{\begin{array}{l}\begin{array}{ll}\dot{z}\left(t\right)=\mathcal{A}z\left(t\right)+BU\left(t\right),\hfill & t>0,\hfill \end{array}\\ z\left(0\right)=z_0.\end{array}$ (2.10)

Here, $B=\left(\begin{array}{c}0\\ c^2\delta \left(x\right)\end{array}\right)$ , where $\delta \left(x\right)$ denotes the Dirac delta distribution.

Theorem 2.1 (Well-posedness of the nonhomogeneous abstract evolution equation) Let $\mathscr{H}=H^1\left(0,1\right)\times L^2\left(0,1\right)$ be equipped with the inner product introduced in (2.1), and consider the operator $\mathcal{A}:D\left(\mathcal{A}\right)\left(\subset \mathscr{H}\right)\to \mathscr{H}$ specified by (2.2) - (2.3). Assume that $\mathcal{A}$ is the generator of an isometric $C_0$ -semigroup ${\left\{S\left(t\right)\right\}}_{t\ge 0}$ on $\mathscr{H}$ , and that $B$ is admissible as a control operator for $\mathcal{A}$ . Equivalently, $B^{*}$ is an admissible observation operator for ${\mathcal{A}}^{*}$ . Then, for every initial datum $z_0={\left(u_0,u_1\right)}^{\top }\in \mathscr{H}$ and every input $U\in L_{\text{loc}}^2\left(\left[0,\infty \right);ℝ\right)$ , the abstract system (2.10) possesses exactly one mild solution $z\in C\left(\left[0,\infty \right);\mathscr{H}\right).$ This solution is represented by the variation-of-constants formula

$z\left(t\right)=S\left(t\right)z_0+{\displaystyle {\int }_0^{t}S\left(t-s\right)BU\left(s\right)\text{d}s}.$ (2.11)

Moreover, for each $T>0$ , there exists a constant $C_T>0$ such that

${\Vert z\left(t\right)\Vert }_{\mathscr{H}}\le C_T\left({\Vert z_0\Vert }_{\mathscr{H}}+{\Vert U\Vert }_{L^2\left(0,T\right)}\right),t\in \left[0,T\right].$ (2.12)

Consequently, the solution varies continuously with respect to both the initial datum $z_0$ and the control input $U$ .

Proof. The proof proceeds in two steps.

Step 1. $\mathcal{A}$ generates an isometric $C_0$ -semigroup on $\mathscr{H}$ .

According to the Lumer-Phillips theorem, it is enough to establish that $\mathcal{A}$ is a densely defined dissipative operator and that the range of $\lambda I-\mathcal{A}$ covers $\mathscr{H}$ for some $\lambda >0$ . The density of $D\left(\mathcal{A}\right)$ in $\mathscr{H}$ follows directly from the definition of the domain. Moreover, for any $z\in D\left(\mathcal{A}\right)$ , we have

$\mathrm{Re}{\langle \mathcal{A}z,z\rangle }_{\mathscr{H}}=0,$ (2.13)

which implies that $\mathcal{A}$ is dissipative. Therefore, it remains only to prove the surjectivity of $\lambda I-\mathcal{A}$ . Let ${\left(h_1,h_2\right)}^{\top }\in \mathscr{H}$ be arbitrary and consider the equation.

$\left(\lambda I-\mathcal{A}\right)\left(\begin{array}{l}u\hfill \\ v\hfill \end{array}\right)=\left(\begin{array}{l}{h}_1\hfill \\ h_2\hfill \end{array}\right),$

which is equivalent to

$\{\begin{array}{l}\lambda u-v=h_1,\\ \lambda v-c^2{u}^{″}=h_2.\end{array}$

Substituting $v=\lambda u-h_1$ into the second equation yields the elliptic boundary value problem

${\lambda }^{2}u-c^2{u}^{″}=\lambda h_1+h_2,u^{\prime }\left(0\right)=0,u^{\prime }\left(1\right)=-qu\left(1\right).$ (2.14)

By the well-posedness theory for second-order boundary value problems, the above equation possesses a unique solution $u\in H^2\left(0,1\right)$ . Consequently, $v=\lambda u-h_1\in H^1\left(0,1\right)$ , which implies that ${\left(u,v\right)}^{\top }\in D\left(\mathcal{A}\right)$ is a solution of the resolvent equation. Therefore, $\lambda I-\mathcal{A}$ is onto. The Lumer-Phillips theorem then guarantees that $\mathcal{A}$ is the generator of an isometric $C_0$ -semigroup ${\left\{S\left(t\right)\right\}}_{t\ge 0}$ on $\mathscr{H}$ .

Step 2. The control operator $B$ is admissible for $\mathcal{A}$ .

Using the duality principle, the admissibility of $B$ is equivalent to that of the observation operator $B^{\text{*}}$ for ${\mathcal{A}}^{\text{*}}$ . Since $\mathcal{A}$ is skew-adjoint, one has ${\mathcal{A}}^{\text{*}}=-\mathcal{A}$ , and $B^{\text{*}}=\left(0,c^2\delta \left(x\right)\right)$ . Hence, $B^{\text{*}}{\left(u,u_t\right)}^{\top }=c^2{u}_t\left(0,t\right)$ . Accordingly, it remains to establish the existence of positive constants $T$ and $K_T$ such that every solution of the dual system satisfies

$\{\begin{array}{l}{v}_{tt}=c^2{v}_{xx},\\ v_x\left(0,t\right)=0,v_x\left(1,t\right)=-qv\left(1,t\right),\\ v\left(x,0\right)=v_0\left(x\right),v_t\left(x,0\right)=v_1\left(x\right)\end{array}$ (2.15)

satisfy

${\displaystyle {\int }_0^T{\left|v_t\left(0,t\right)\right|}^2\text{d}t}\le K_T{\Vert \left(v_0,v_1\right)\Vert }_{\mathscr{H}}^2.$ (2.16)

Define the energy

$E\left(t\right)=\frac{1}{2}{\Vert \left(v,v_t\right)\Vert }_{\mathscr{H}}^2=\frac{1}{2}\left({\displaystyle {\int }_0^1{c}^2{\left|v_x\right|}^2\text{d}x}+{\displaystyle {\int }_0^1{\left|v_t\right|}^2\text{d}x}+c^{2}q{\left|v\left(1\right)\right|}^2\right).$ (2.17)

Using the multiplier $\left(x-1\right)v_x$ , we derive the following identity:

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^1\left(x-1\right)v_x{v}_t\text{d}x}=\frac{1}{2}{v}_t^2\left(0,t\right)-\frac{1}{2}{\displaystyle {\int }_0^1\left(v_t^2+c^2{v}_x^2\right)\text{d}x}.$ (2.18)

Upon integrating both sides of equation (2.18) with respect to $t$ over the interval $\left[0,T\right]$ , one arrives at

${\displaystyle {\int }_0^T\left[\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^1\left(x-1\right)v_x{v}_t\text{d}x}\right]\text{d}t}=\frac{1}{2}{\displaystyle {\int }_0^T{v}_t^2\left(0,t\right)\text{d}t}-\frac{1}{2}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^1\left(v_t^2+c^2{v}_x^2\right)\text{d}x\text{d}t}}.$ (2.19)

Owing to the isometric dissipativity property, $E\left(t\right)=E\left(0\right)$ holds for every $t\ge 0$ . Observing that $\left|x-1\right|\le 1$ , it follows that

$\left|{\displaystyle {\int }_0^1\left(x-1\right)v_x{v}_t\text{d}x}\right|\le \frac{1}{2}{\displaystyle {\int }_0^1\left(v_x^2+v_t^2\right)\text{d}x}\le \frac{1}{\text{min}\left\{1,c^2\right\}}E\left(t\right)=\frac{1}{\text{min}\left\{1,c^2\right\}}E\left(0\right),$

which implies that the left-hand term in equation (2.19) satisfies

$\left|{\displaystyle {\int }_0^T\left[\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^1\left(x-1\right)v_x{v}_t\text{d}x}\right]\text{d}t}\right|\le \frac{1}{\text{min}\left\{1,c^2\right\}}E\left(0\right).$ (2.20)

Moreover, by (2.17), we obtain

${\displaystyle {\int }_0^1\left(v_t^2+c^2{v}_x^2\right)\text{d}x}=2E\left(t\right)-c^{2}q{\left|v\left(1\right)\right|}^2\le 2E\left(t\right)=2E\left(0\right),q>0.$ (2.21)

Substituting (2.20) and (2.21) into the integral equality (2.19) and rearranging terms, we obtain

$\begin{array}{c}\frac{1}{2}{\displaystyle {\int }_0^T{v}_t^2\left(0,t\right)\text{d}t}=\frac{1}{2}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^1\left(v_t^2+c^2{v}_x^2\right)\text{d}x\text{d}t}}+{\displaystyle {\int }_0^T\left[\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_0^1\left(x-1\right)v_x{v}_t\text{d}x}\right]\text{d}t}\\ \le \frac{1}{2}{\displaystyle {\int }_0^{T}2E\left(0\right)\text{d}t}+\frac{1}{\text{min}\left\{1,c^2\right\}}E\left(0\right)\\ =TE\left(0\right)+\frac{1}{\text{min}\left\{1,c^2\right\}}E\left(0\right).\end{array}$ (2.22)

Upon multiplying both sides of equation (2.22) by a factor of 2, it follows that

${\displaystyle {\int }_0^T{v}_t^2\left(0,t\right)\text{d}t}\le 2K_{T}E\left(0\right),$ (2.23)

where the constant $K_T=T+\frac{1}{\text{min}\left\{1,c^2\right\}}$ . By definition (2.17), we have

$E\left(0\right)=\frac{1}{2}{\Vert \left(u_0,u_1\right)\Vert }_{\mathscr{H}}^2.$ (2.24)

Combining this with inequality (2.23), we immediately obtain

${\displaystyle {\int }_0^T{\left|v_t\left(0,t\right)\right|}^2\text{d}t}\le K_T{\Vert \left(v_0,v_1\right)\Vert }_{\mathscr{H}}^2.$

Therefore, the admissibility of $B^{\text{*}}$ for ${\mathcal{A}}^{\text{*}}$ directly leads to the admissibility of $B$ as a control operator for $\mathcal{A}$ .

Since $\mathcal{A}$ generates an isometric $C_0$ -semigroup and $B$ is admissible, the standard theory of abstract control systems ensures that the nonhomogeneous evolution equation possesses a unique mild solution $z\in C\left(\left[0,\infty \right);\mathscr{H}\right)$ , which is represented by the variation-of-constants formula (2.11). Moreover, the solution satisfies the energy estimate (2.12), guaranteeing its continuous dependence on both the initial condition and the control input.

We expand the initial state $z_0={\left(u_0,u_1\right)}^{\top }\in \mathscr{H}$ with respect to the orthonormal eigenbasis ${\left\{{\tilde{\varphi }}_n^{\pm }\right\}}_{n\ge 1}$ as

$z_0={\displaystyle \sum }_{n=1}^{\infty }\left(a_n^{+}{\tilde{\varphi }}_n^{+}+a_n^{-}{\tilde{\varphi }}_n^{-}\right),$ (2.25)

where $a_n^{\pm }={\langle z_0,{\tilde{\varphi }}_n^{\pm }\rangle }_{\mathscr{H}}$ denote the initial modal coefficients.

We rely on a fundamental property satisfied by the isometric $C_0$ -semigroup ${\left\{S\left(t\right)\right\}}_{t\ge 0}$ associated with the operator $\mathcal{A}$ , that is,

$S\left(t\right){\tilde{\varphi }}_n^{\pm }={\text{e}}^{{\lambda }_n^{\pm }t}{\tilde{\varphi }}_n^{\pm }.$ (2.26)

Combining the variation-of-constants formula (2.11) for the system solution, we compute the observations contributed by the initial state term and the control term, respectively, and finally obtain the expression for the boundary observation $y\left(t\right)$ at $x=0$ as

$\begin{array}{c}y\left(t\right)=y_0\left(t\right)+y_U\left(t\right)\\ ={\displaystyle \sum }_{n=1}^{\infty }\frac{1}{{\alpha }_n}\left(a_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+a_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right)+{\displaystyle {\int }_0^{t}G\left(t-s,0\right)U\left(s\right)\text{d}s}\\ \triangleq u\left(0,t;0,u_0,u_1\right)+u\left(0,t;U,0,0\right),\end{array}$ (2.27)

where the Green’s function is expressed as

$G\left(t-s,0\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{\text{i}{\omega }_n{c}^3}{{\alpha }_n^2}\left({\text{e}}^{{\lambda }_n^{-}\left(t-s\right)}-{\text{e}}^{{\lambda }_n^{+}\left(t-s\right)}\right).$ (2.28)

Assume that $T_2>T_1>0$ and that the boundary control vanishes identically, i.e., $U\left(t\right)=0$ for every $t\in \left[0,T_2\right]$ . Under this assumption, the boundary observation takes the form

$y\left(t\right)=u\left(0,t;0,u_0,u_1\right)={\displaystyle \sum }_{n=1}^{\infty }\left(C_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+C_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right),t\in \left[T_1,T_2\right],$ (2.29)

in which the coefficients are given by $C_n^{\pm }=\frac{a_n^{\pm }}{{\alpha }_n}$ for each $n\in {\mathbb{N}}^{\text{*}}$ .

Because the initial data $u_0$ and $u_1$ are not known a priori, one cannot determine in advance which coefficients $C_n^{\pm }$ are nonzero. We define an unknown set $\mathbb{K}\subset {\mathbb{N}}^{\text{*}}$ such that $C_k^{\pm }\ne 0$ for $k\in \mathbb{K}$ and $C_k^{\pm }=0$ for $k\notin \mathbb{K}$ .

Theorem 2.2 Assume that $0<T_1<T_2<\infty$ , $u_0\in H^1\left(0,1\right)$ , and $u_1\in L^2\left(0,1\right)$ . Then the eigenvalues and corresponding coefficients ${\left\{\left({\lambda }_k^{+},C_k^{+}\right),\left({\lambda }_k^{-},C_k^{-}\right)\right\}}_{k\in \mathbb{K}}$ can be uniquely determined from the observation $y\left(t\right)$ on the interval $\left[T_1,T_2\right]$ .

Proof. The zero-input boundary response admits the nonharmonic Fourier series representation

$y\left(t\right)={\displaystyle \sum }_{n=1}^{\infty }\left(C_n^{+}{\text{e}}^{\text{i}{\omega }_{n}ct}+C_n^{-}{\text{e}}^{-\text{i}{\omega }_{n}ct}\right),t\ge 0,$ (2.30)

where the frequencies $\pm {\omega }_{n}c$ are purely imaginary and derived from the characteristic equation $\omega \tan \omega =q$ . By the properties of the Sturm-Liouville problem, the sequence ${\left\{{\omega }_n\right\}}_{n\ge 1}$ is uniformly separated: there exists a constant $\gamma >0$ depending only on $q$ such that

$\underset{n\ne m}{\text{inf}}\left|{\omega }_n-{\omega }_m\right|\ge \gamma >0.$ (2.31)

Consequently, the frequency set $\Lambda ={\left\{\pm {\omega }_{n}c\right\}}_{n\ge 1}$ is also uniformly separated with minimal spacing $c\gamma >0$ .

The regularity assumptions $u_0\in H^1\left(0,1\right)$ and $u_1\in L^2\left(0,1\right)$ imply that the coefficients satisfy

${\displaystyle \sum }_{n=1}^{\infty }{n}^2\left({\left|C_n^{+}\right|}^2+{\left|C_n^{-}\right|}^2\right)<\infty ,$ (2.32)

which guarantees absolute and uniform convergence of the series on any bounded time interval.

We now invoke the Ingham inequality for uniformly separated frequencies. If the observation interval length satisfies

$T_2-T_1>\frac{2\pi }{c\gamma },$ (2.33)

then there exists a constant $C_{\text{Ing}}>0$ such that

${\displaystyle {\int }_{T_1}^{T_2}{\left|{\displaystyle \sum }_{n=1}^{\infty }\left(C_n^{+}{\text{e}}^{\text{i}{\omega }_{n}ct}+C_n^{-}{\text{e}}^{-\text{i}{\omega }_{n}ct}\right)\right|}^2\text{d}t}\ge C_{\text{Ing}}{\displaystyle \sum }_{n=1}^{\infty }\left({\left|C_n^{+}\right|}^2+{\left|C_n^{-}\right|}^2\right).$ (2.34)

Suppose there exist two sets of spectral parameters $\left\{{\lambda }_k^{\pm },C_k^{\pm }\right\}$ and $\left\{{\tilde{\lambda }}_k^{\pm },{\tilde{C}}_k^{\pm }\right\}$ that generate identical boundary observations on $\left[T_1,T_2\right]$ . Subtracting the two representations yields a series that vanishes identically on $\left[T_1,T_2\right]$ . Applying the Ingham inequality to the difference series forces all coefficients to be zero, which implies that the two sets of frequencies (poles) must coincide exactly, and the corresponding modal coefficients are equal.

Therefore, the eigenvalues and corresponding coefficients

${\left\{\left({\lambda }_k^{+},C_k^{+}\right),\left({\lambda }_k^{-},C_k^{-}\right)\right\}}_{k\in \mathbb{K}}$ are uniquely determined by the boundary observation $y\left(t\right)$ on $\left[T_1,T_2\right]$ , provided the observation interval is sufficiently long.

Theorem 2.3 Assume that $0<T_1<T_2<T_3<\infty$ , $u_0\in H^1\left(0,1\right)$ , and $u_1\in L^2\left(0,1\right)$ . Let the control input $U\in L^2\left(0,T_3\right)$ satisfy

$\{\begin{array}{ll}U\left(t\right)=0,\hfill & t\in \left[0,T_2\right),\hfill \\ U\left(t\right)\ne 0,\hfill & \text{for}\text{almost}\text{all}t\in \left[T_2,T_3\right].\hfill \end{array}$ (2.35)

Denote the associated boundary measurement by

$y\left(t\right)=u\left(0,t;U,u_0,u_1\right),t\in \left[T_1,T_3\right].$

Then, the wave speed $c$ and the initial states $u_0\left(x\right)$ and $u_1\left(x\right)$ can be uniquely determined from the boundary observation $y\left(t\right)$ on the interval $\left[T_1,T_3\right]$ .

Proof. From equation (2.27), we define the difference observation

$\tilde{y}\left(t\right)\triangleq {\displaystyle \sum }_{k\in \mathbb{K}}\left[C_k^{+}{\text{e}}^{{\lambda }_k^{+}\left(t+T_2\right)}+C_k^{-}{\text{e}}^{{\lambda }_k^{-}\left(t+T_2\right)}\right]-y\left(t+T_2\right),t\in \left(0,T_3-T_2\right].$ (2.36)

Let $\tilde{U}\left(t\right)=U\left(t+T_2\right)$ for $t\in \left(0,T_3-T_2\right]$ , and denote the Green’s function as

$G\left(t,0\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{\text{i}{\omega }_n{c}^3}{{\alpha }_n^2}\left({\text{e}}^{{\lambda }_n^{-}t}-{\text{e}}^{{\lambda }_n^{+}t}\right)\triangleq {\displaystyle \sum }_{n=1}^{\infty }{G}_n\left({\text{e}}^{{\lambda }_n^{-}t}-{\text{e}}^{{\lambda }_n^{+}t}\right),t\in \left(0,T_3-T_2\right].$ (2.37)

Substituting these into the observation formula yields

$\tilde{y}\left(t\right)=-{\displaystyle {\int }_0^{t}G\left(t-s,0\right)\tilde{U}\left(s\right)\text{d}s},t\in \left(0,T_3-T_2\right].$ (2.38)

We first justify the continuity of the Green function $G\left(t,0\right)$ and then establish its unique identifiability from finite-time convolution data. Substituting the spectral expansion and normalization constants into the expression for $G\left(t,0\right)$ , we obtain the trigonometric series

$G\left(t,0\right)={\displaystyle \sum }_{n=1}^{\infty }\frac{2c}{{\omega }_n+\frac{1}{2}\text{sin}\left(2{\omega }_n\right)}\text{sin}\left({\omega }_{n}ct\right).$ (2.39)

From the Sturm-Liouville eigenvalue condition ${\omega }_n\tan {\omega }_n=q$ , we have ${\omega }_n~\left(n-\frac{1}{2}\right)\pi$ as $n\to \infty$ , which yields coefficient decay

$A_n:=\frac{2c}{{\omega }_n+\frac{1}{2}\text{sin}\left(2{\omega }_n\right)}=O\left(n^{-1}\right).$ (2.40)

The sequence $\left\{A_n\right\}$ is positive and monotonically decreasing toward zero. By the Dirichlet test for trigonometric series, the series converges uniformly on every compact time interval $\left[0,T\right]$ . Consequently,

$G\left(\cdot ,0\right)\in C\left(\left[0,T\right]\right)\subset L^1\left(0,T\right),\text{for}\text{any}T>0,$

and term-by-term evaluation gives $G\left(0,0\right)=0$ .

Recall the shifted input $\tilde{U}\left(t\right)=U\left(t+T_2\right)$ for $t\in \left(0,T_3-T_2\right]$ . From the assumption $U\in L^2\left(0,T_3\right)$ , we obtain $\tilde{U}\in L^2\left(0,T_3-T_2\right)\subset L^1\left(0,T_3-T_2\right)$ . Both $G\left(\cdot ,0\right)$ and $\tilde{U}$ satisfy the function-space hypotheses of the Titchmarsh convolution theorem ([24] Chap. 1): for functions $f,g\in L^1\left(0,T\right)$ , if the convolution

$\left(f\ast g\right)\left(t\right)={\displaystyle {\int }_0^{t}f\left(t-s\right)g\left(s\right)\text{d}s}$ (2.41)

vanishes for almost every $t\in \left(0,T\right)$ , then there exists $\alpha \in \left[0,T\right]$ such that $f=0$ almost everywhere on $\left(0,\alpha \right)$ and $g=0$ almost everywhere on $\left(0,T-\alpha \right)$ .

We now deduce the deconvolution uniqueness needed for our argument. Suppose two candidates $G_1,G_2\in L^1\left(0,T_3-T_2\right)$ satisfy $G_1\ast \tilde{U}=G_2\ast \tilde{U}$ almost everywhere on $\left(0,T_3-T_2\right]$ . Define $f:=G_1-G_2$ , so that $f\ast \tilde{U}=0$ a.e. Applying the Titchmarsh result, there exists $\alpha \in \left[0,T_3-T_2\right]$ such that $f=0$ a.e. on $\left(0,\alpha \right)$ and $\tilde{U}=0$ a.e. on $\left(0,T_3-T_2-\alpha \right)$ . By hypothesis, $\tilde{U}\left(t\right)\ne 0$ for almost every $t\in \left(0,T_3-T_2\right]$ , which forces $T_3-T_2-\alpha =0$ , i.e., $\alpha =T_3-T_2$ . Therefore, $f=0$ a.e. on $\left(0,T_3-T_2\right]$ , meaning $G_1=G_2$ almost everywhere. Since $G_1,G_2$ are continuous functions, almost-everywhere equality implies pointwise equality for all $t\in \left(0,T_3-T_2\right]$ . Hence, $G\left(t,0\right)$ is uniquely determined from the convolution relation.

Moreover, according to Theorem 2.2, the spectral parameters ${\left\{\left({\lambda }_k^{+},C_k^{+}\right),\left({\lambda }_k^{-},C_k^{-}\right)\right\}}_{k\in \mathbb{K}}$ are determined by the boundary measurement $\left\{y\left(t\right)|t\in \left[T_1,T_2\right]\right\}$ . Hence, the extended observation data $\left\{\tilde{y}\left(t\right)|t\in \left(0,T_3-T_2\right]\right\}$ is uniquely determined by the original measurement $\left\{y\left(t\right)|t\in \left[T_1,T_3\right]\right\}$ . Since $G_n\ne 0$ for all $n$ , applying Theorem 2.2 yields that the eigenvalue sequence ${\left\{{\lambda }_n^{+}\right\}}_{n\in {\mathbb{N}}^{*}}$ is identifiable from $\left\{G\left(t,0\right)|t\in \left(0,T_3-T_2\right]\right\}$ . Combining this result with the relation ${\lambda }_n^{\pm }=\pm i{\omega }_{n}c$ , the wave speed $c$ can consequently be identified uniquely.

We now proceed to investigate whether the initial data $u_0\left(x\right)$ and $u_1\left(x\right)$ can be uniquely determined. Substituting the identified $c$ back into equation (2.27), we rearrange it as

$y\left(t\right)+{\displaystyle {\int }_0^{t}G\left(t-s,0\right)U\left(s\right)\text{d}s}={\displaystyle \sum }_{n=1}^{\infty }\frac{1}{{\alpha }_n}\left(a_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+a_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right),t\in \left[T_1,T_3\right].$ (2.42)

The identification of $c$ from the measured output $\left\{y\left(t\right)|t\in \left[T_1,T_3\right]\right\}$ enables us to evaluate the corresponding spectral relation completely. In this relation, the unknown terms appear in the form of a complex exponential expansion. The uniqueness property of such spectral representations guarantees that the expansion coefficients $\left\{\frac{1}{{\alpha }_n}\left(a_n^{+},a_n^{-}\right)\right\}$ can be uniquely obtained from the available observation data. Therefore, the modal coefficients associated with the initial conditions are determined.

$u_0\left(x\right)={\displaystyle \sum }_{n=1}^{\infty }\langle u_0,{\varphi }_n\rangle {\varphi }_n\left(x\right),u_1\left(x\right)={\displaystyle \sum }_{n=1}^{\infty }\langle u_1,{\varphi }_n\rangle {\varphi }_n\left(x\right),$ (2.43)

Here, the Fourier coefficients $\langle u_0,{\varphi }_n\rangle$ and $\langle u_1,{\varphi }_n\rangle$ can be calculated from the identified spectral quantities $\left\{{\lambda }_n^{\pm }\right\}$ and $\left\{a_n^{\pm }\right\}$ . Hence, the initial profiles $u_0\left(x\right)$ and $u_1\left(x\right)$ are uniquely determined from the boundary measurement $\left\{y\left(t\right)|t\in \left[T_1,T_3\right]\right\}$ .

3. Numerical Computation Method

From the analysis in the preceding sections, it is evident that the core challenge of the identification problem lies in retrieving the spectral-coefficient pairs ${\left\{\left(C_n^{+},{\lambda }_n^{+}\right),\left(C_n^{-},{\lambda }_n^{-}\right)\right\}}_{n\in {\mathbb{N}}^{\text{*}}}$ from the boundary measurements $y\left(t\right)$ on the interval $\left[T_1,T_2\right]$ , based on the complex exponential series representation given in (2.27). A primary obstacle lies in the fact that an infinite number of coefficients $C_n^{\pm }$ in (2.27) may be nonzero. To address this, we adopt the matrix pencil technique in this section to extract a subset of the spectral-coefficient data $\left\{\left(C_n^{+},{\lambda }_n^{+}\right),\left(C_n^{-},{\lambda }_n^{-}\right)\right\}$ from the series truncated to its first $M$ terms in (2.27), with the remaining higher-order components treated as measurement noise.

3.1. Bounded-Rank Truncation

Let $M\in {\mathbb{N}}^{\text{*}}$ be a positive integer. We split the complex exponential series in (2.29) into two components:

$y\left(t\right)={\displaystyle \sum }_{n=1}^M\left(C_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+C_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right)+e\left(M+1,t\right),t\in \left[T_1,T_2\right],$ (3.1)

where the remainder term is defined as

$e\left(M+1,t\right)={\displaystyle \sum }_{n=M+1}^{\infty }\left(C_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+C_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right),t\in \left[T_1,T_2\right].$ (3.2)

The following theorem establishes a uniform estimate for the remainder $e\left(M+1,t\right)$ .

Theorem 3.1 Suppose that $c$ is bounded below by a positive constant $c_0$ , and let $\left(u_0,u_1\right)\in H^1\left(0,1\right)\times L^2\left(0,1\right)$ . At any time $t\in \left[T_1,T_2\right]$ , the error associated with the $\left(M+1\right)$ -st truncation level obeys

$\left|e\left(M+1,t\right)\right|\le C{M}^{-1/2},$ (3.3)

where $C$ is a positive constant whose value remains unchanged as $M$ varies.

Proof. As $n\to \infty$ , the characteristic frequencies satisfy ${\omega }_n~n\pi$ . From the definition of ${\alpha }_n$ , we have

${\alpha }_n=c\sqrt{{\omega }_n\left({\omega }_n+\frac{1}{2}\text{sin}\left(2{\omega }_n\right)\right)}~cn\pi .$ (3.4)

Applying Parseval’s theorem to the complete orthonormal system $\left\{{\tilde{\varphi }}_n^{\pm }\right\}$ of $\mathscr{H}$ , we obtain

${\Vert z_0\Vert }_{\mathscr{H}}^2={\displaystyle \sum }_{n=1}^{\infty }\left({\left|a_n^{+}\right|}^2+{\left|a_n^{-}\right|}^2\right)<\infty ,$ (3.5)

where $z_0={\left(u_0,u_1\right)}^{\top }$ . Thus, for any $M\in {\mathbb{N}}^{\text{*}}$ ,

${\displaystyle \sum }_{n=M+1}^{\infty }{\left|a_n^{\pm }\right|}^2\le {\Vert z_0\Vert }_{\mathscr{H}}^2<\infty .$ (3.6)

For any $t\in \left[T_1,T_3\right]$ , note that $\left|{\text{e}}^{{\lambda }_n^{\pm }t}\right|=1$ . By the triangle inequality,

$\begin{array}{c}\left|e\left(M+1,t\right)\right|=\left|{\displaystyle \sum }_{n=M+1}^{\infty }\left(C_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+C_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right)\right|\\ \le {\displaystyle \sum }_{n=M+1}^{\infty }\left(\left|C_n^{+}{\text{e}}^{{\lambda }_n^{+}t}\right|+\left|C_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right|\right)\\ ={\displaystyle \sum }_{n=M+1}^{\infty }\left(\left|C_n^{+}\right|+\left|C_n^{-}\right|\right)\\ ={\displaystyle \sum }_{n=M+1}^{\infty }\left(\frac{\left|a_n^{+}\right|}{{\alpha }_n}+\frac{\left|a_n^{-}\right|}{{\alpha }_n}\right).\end{array}$ (3.7)

The tail of the coefficient series can be controlled by the Cauchy-Schwarz estimate as follows

${\displaystyle \sum }_{n=M+1}^{\infty }\frac{\left|a_n^{\pm }\right|}{{\alpha }_n}\le \frac{1}{c_0\pi }{\left({\displaystyle \sum }_{n=M+1}^{\infty }\frac{1}{n^2}\right)}^{\frac{1}{2}}{\left({\displaystyle \sum }_{n=M+1}^{\infty }{\left|a_n^{\pm }\right|}^2\right)}^{\frac{1}{2}}\le \frac{{\Vert z_0\Vert }_{\mathscr{H}}{M}^{-1/2}}{c_0}.$ (3.8)

Let $C=\frac{2{\Vert z_0\Vert }_{\mathscr{H}}}{c_0}$ . Then

$\left|e\left(M+1,t\right)\right|\le C{M}^{-1/2},$

where $C$ is a positive constant whose value remains unchanged as $M$ varies.

3.2. Identification Algorithm

Let $0<T_1<T_2<T_3$ be three arbitrary positive constants, and let the control function $U\left(t\right)$ be chosen as the following on-off control:

$U\left(t\right)=\{\begin{array}{ll}0,\hfill & t\in \left(0,T_2\right),\hfill \\ 1,\hfill & t\in \left[T_2,T_3\right].\hfill \end{array}$ (3.9)

The corresponding boundary observation data is denoted by $\left\{y\left(t\right)=u\left(0,t;U,u_0,u_1\right)|t\in \left[T_1,T_3\right]\right\}$ .

It can be seen from the proofs of Theorem 2.2 and Theorem 2.3 on identifiability analysis that the core of the identification algorithm lies in estimating the conjugate modal pairs ${\left\{\left({\lambda }_k^{+},C_k^{+}\right),\left({\lambda }_k^{-},C_k^{-}\right)\right\}}_{k\in \mathbb{K}}$ from the finite-time observation $\left\{y\left(t\right)={\displaystyle {\sum }_{k\in \mathbb{K}}\left(C_k^{+}{\text{e}}^{{\lambda }_k^{+}t}+C_k^{-}{\text{e}}^{{\lambda }_k^{-}t}\right)}|t\in \left[T_1,T_2\right]\right\}$ . As indicated by (2.27), the boundary output $y\left(t\right)$ of system (1.1) can be decomposed into two components $u\left(0,t;0,u_0,u_1\right)$ and $u\left(0,t;U,0,0\right)$ . Once ${\left\{\left({\lambda }_k^{+},C_k^{+}\right),\left({\lambda }_k^{-},C_k^{-}\right)\right\}}_{k\in \mathbb{K}}$ are identified, the uncontrolled free response component, which is determined by the unknown initial states, $u\left(0,t;0,u_0,u_1\right)$ , can be eliminated from the boundary observation $y\left(t\right)$ . This allows the identification problem for the wave speed $c$ to be equivalently transformed into a problem with zero initial states. After estimating the wave speed $c$ , the remaining task consists of reconstructing the initial states.

Based on the above analysis, the identification algorithm for the wave speed and initial states can be summarized into the following four steps.

3.2.1. Step 1: Matrix-Pencil Recovery of Selected Eigenvalues of the System Operator $A$ from Uncontrolled Boundary Data

The boundary input is deactivated over the interval $\left[0,T_2\right]$ , that is, $U\left(t\right)=0$ . Consequently, for $t\in \left[T_1,T_2\right]$ , the measured boundary observation admits the modal representation

$\begin{array}{c}u\left(0,t;0,u_0,u_1\right)={\displaystyle \sum }_{n=1}^{\infty }\left(C_n^{+}{\text{e}}^{{\lambda }_n^{+}t}+C_n^{-}{\text{e}}^{{\lambda }_n^{-}t}\right)\\ ={\displaystyle \sum }_{k=0}^{M-1}\left(C_{n_k}^{+}{\text{e}}^{{\lambda }_{n_k}^{+}t}+C_{n_k}^{-}{\text{e}}^{{\lambda }_{n_k}^{-}t}\right)+r_M\left(t\right),t\in \left[T_1,T_2\right].\end{array}$ (3.10)

Let $K_c:=\#\mathbb{K}$ , where $\mathbb{K}$ is the index set defined previously and $K_c$ may be infinite. We enumerate the elements of $\mathbb{K}$ as

$\mathbb{K}=\left\{n_0,n_1,\cdots \right\},$

so that ${\left\{C_{n_k}^{\pm }\right\}}_{k=0}^{K_c-1}$ comprises all the nonzero spectral coefficient pairs in ${\left\{C_n^{\pm }\right\}}_{n\in {\mathbb{N}}^{\text{*}}}$ . The contribution of the observable modes beyond the first $M$ pairs is represented by

$r_M\left(t\right):={\displaystyle \sum }_{k=M}^{K_{\text{c}}-1}\left(C_{n_k}^{+}{\text{e}}^{{\lambda }_{n_k}^{+}t}+C_{n_k}^{-}{\text{e}}^{{\lambda }_{n_k}^{-}t}\right),$ (3.11)

with the upper limit interpreted as infinity when $K_{\text{c}}=\infty$ .

The truncation estimate established in Theorem 3.1 implies that

$\underset{t\in \left[T_1,T_2\right]}{\text{sup}}\left|r_M\left(t\right)\right|\le C{M}^{-1/2},$ (3.12)

for a positive constant $C$ that does not depend on $M$ . Hence, $r_M\left(t\right)$ converges uniformly to zero on $\left[T_1,T_2\right]$ as $M\to \infty$ . The uncontrolled boundary observation can therefore be approximated by the finite exponential expansion

$u\left(0,t;0,u_0,u_1\right)\approx {\displaystyle \sum }_{k=0}^{M-1}\left(C_{n_k}^{+}{\text{e}}^{{\lambda }_{n_k}^{+}t}+C_{n_k}^{-}{\text{e}}^{{\lambda }_{n_k}^{-}t}\right),t\in \left[T_1,T_2\right].$ (3.13)

To identify the exponents in (3.13), the boundary observation is sampled on a uniform temporal grid

$t_i=T_1+i\Delta t_1,i=0,1,\cdots ,N_1-1,\Delta t_1=\frac{T_2-T_1}{N_1-1}.$

Thus, $N_1$ samples are available. Setting

$y_i:=u\left(0,t_i;0,u_0,u_1\right),$ (3.14)

together with

$R_k^{\pm }:=C_{n_k}^{\pm }{\text{e}}^{{\lambda }_{n_k}^{\pm }{T}_1},z_k^{\pm }:={\text{e}}^{{\lambda }_{n_k}^{\pm }\Delta t_1},$ (3.15)

gives the discrete observation model

$y_i={\displaystyle \sum }_{k=0}^{M-1}\left[R_k^{+}{\left(z_k^{+}\right)}^i+R_k^{-}{\left(z_k^{-}\right)}^i\right]+{\epsilon }_i,i=0,1,\cdots ,N_1-1,$ (3.16)

where ${\epsilon }_i:=r_M\left(t_i\right)$ is the error generated by the neglected spectral tail. In particular, Theorem 3.1 yields

$\underset{0\le i\le N_1-1}{\text{max}}\left|{\epsilon }_i\right|=\mathcal{O}\left(M^{-1/2}\right).$ (3.17)

Therefore, the sampled boundary observation has the form of a finite sum of discrete exponentials perturbed by an additive error, which permits the direct application of a matrix-pencil construction.

Since the initial state $z_0={\left(u_0,u_1\right)}^{\text{T}}$ is real-valued and the normalized eigenvectors satisfy ${\tilde{\varphi }}_n^{-}=\overline{{\tilde{\varphi }}_n^{+}}$ , the initial modal coefficients obey $a_n^{-}=\overline{a_n^{+}}$ . Moreover, because ${\alpha }_n$ is positive and real and $C_n^{\pm }=a_n^{\pm }/{\alpha }_n$ , it follows that $C_n^{-}=\overline{C_n^{+}}$ . Consequently,

$C_n^{+}=0\iff C_n^{-}=0.$

Thus, each observable mode contributes a pair of nonzero exponential terms. Provided that the corresponding discrete poles are mutually distinct, the truncated representation in (3.16) contains $P=2M$ discrete poles.

Choose the pencil parameter $L$ such that $P\le L\le N_1-P$ . In the numerical implementation, $L$ may be set to $N_1/3$ , provided that this value lies within the admissible range. The shifted Hankel matrices constructed from the measured samples are

$Y_0=\left[\begin{array}{cccc}{y}_0& y_1& \cdots & y_{L-1}\\ y_1& y_2& \cdots & y_L\\ ⋮& ⋮& \ddots & ⋮\\ y_{N_1-L-1}& y_{N_1-L}& \cdots & y_{N_1-2}\end{array}\right],Y_1=\left[\begin{array}{cccc}{y}_1& y_2& \cdots & y_L\\ y_2& y_3& \cdots & y_{L+1}\\ ⋮& ⋮& \ddots & ⋮\\ y_{N_1-L}& y_{N_1-L+1}& \cdots & y_{N_1-1}\end{array}\right].$ (3.18)

In the absence of perturbations, $Y_0$ has rank $P$ , and the discrete poles are determined by the shift relation between $Y_0$ and $Y_1$ . For measured data, the effective model order is inferred from the singular-value decomposition

$Y_0=U\Sigma V^{\text{*}},\Sigma =\text{diag}\left({\sigma }_1,{\sigma }_2,\cdots \right),{\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0.$ (3.19)

Given a prescribed relative threshold ${\epsilon }_{\text{svd}}>0$ , the number of identifiable discrete poles is estimated by

$\widehat{P}=\#\left\{j:\frac{{\sigma }_j}{{\sigma }_1}\ge {\epsilon }_{\text{svd}}\right\}.$ (3.20)

After imposing the paired spectral structure, let $\widehat{M}$ denote the number of retained complete pole pairs. If all the identified poles admit valid conjugate partners, then

$\widehat{P}=2\widehat{M}.$ (3.21)

Let $U_{\widehat{P}}$ , ${\Sigma }_{\widehat{P}}$ , and $V_{\widehat{P}}$ denote the singular-vector and singular-value matrices associated with the $\widehat{P}$ dominant singular values. The reduced matrix pencil is then represented by

$Z_{\text{E}}={\Sigma }_{\widehat{P}}^{-1}{U}_{\widehat{P}}^{\text{*}}{Y}_1{V}_{\widehat{P}}.$ (3.22)

Its eigenvalues ${\left\{{\widehat{z}}_j\right\}}_{j=1}^{\widehat{P}}$ provide estimates of the discrete poles $\left\{z_k^{+},z_k^{-}\right\}$ . The corresponding continuous-time exponents are determined through

${\widehat{\lambda }}_j=\frac{\log {\widehat{z}}_j}{\Delta t_1},j=1,2,\cdots ,\widehat{P},$ (3.23)

where the logarithmic branch is selected consistently with the spectral location of the system operator.

The mapping from discrete-time poles to continuous-time exponents via the complex logarithm is multi-valued. Consistent selection of the logarithmic branch alone does not guarantee unique identification of the continuous-time exponents from the sampled poles. To prevent aliasing of the retained modal frequencies over the prescribed wave-speed search interval $\left[c_{\text{min}},c_{\text{max}}\right]$ , the sampling interval must satisfy the following anti-aliasing condition.

Let ${\omega }_{\widehat{M}}$ denote the angular wave number of the highest retained mode among the identified conjugate pole pairs. Since the discrete poles are given by

${\widehat{z}}_k^{\pm }=\exp \left(\pm \text{i}{\omega }_{k}c\Delta t_1\right),$ (3.24)

the sampling interval is selected such that

${\omega }_{\widehat{M}}{c}_{\max }\Delta t_1<\pi .$ (3.25)

Equivalently,

$\Delta t_1<\frac{\pi }{{\omega }_{\widehat{M}}{c}_{\max }}.$ (3.26)

Under this condition, all retained modal frequencies remain inside the principal branch of the complex logarithm, and the continuous-time exponents are uniquely determined by the sampled poles.

After arranging the identified exponents into their + and − branches, one obtains

${\left\{{\widehat{\lambda }}_{n_k}^{+},{\widehat{\lambda }}_{n_k}^{-}\right\}}_{k=0}^{\widehat{M}-1},$

which constitutes the finite spectral information required in the subsequent identification steps. Once these exponents have been determined, the associated modal coefficients can be computed from (3.16) by a linear least-squares fit. The underlying matrix-pencil principle may be found in [32].

3.2.2. Step 2: Linear Least-Squares Estimation of the Spectral Coefficients

The eigenvalue estimates obtained in Step 1 are held fixed. Since the initial state is real-valued, the spectral coefficients associated with the two conjugate branches are not independent. In particular, for each retained mode,

${\widehat{\lambda }}_{n_k}^{-}=\overline{{\widehat{\lambda }}_{n_k}^{+}},C_{n_k}^{-}=\overline{C_{n_k}^{+}},k=0,1,\cdots ,\widehat{M}-1.$

Therefore, only the coefficients of the positive-frequency branch are treated as independent unknowns. Let

$C={\left[\begin{array}{llll}{C}_{n_0}^{+}\hfill & C_{n_1}^{+}\hfill & \cdots \hfill & C_{n_{\widehat{M}-1}}^{+}\hfill \end{array}\right]}^{\text{T}}\in {ℂ}^{\widehat{M}},$ (3.27)

and define the observation vector by

$y={\left[\begin{array}{llll}{y}_0\hfill & y_1\hfill & \cdots \hfill & y_{N_1-1}\hfill \end{array}\right]}^{\text{T}}\in {ℝ}^{N_1}.$ (3.28)

The zero-input boundary observation can then be written as

$y\left(t_i\right)=2\mathrm{Re}\left({\displaystyle \sum }_{k=0}^{\widehat{M}-1}{C}_{n_k}^{+}{\text{e}}^{{\widehat{\lambda }}_{n_k}^{+}{t}_i}\right)+{\epsilon }_i,i=0,1,\cdots ,N_1-1.$ (3.29)

Define the regression matrix $\Phi \in {ℂ}^{N_1\times \widehat{M}}$ by

${\Phi }_{i+1,k+1}={\text{e}}^{{\widehat{\lambda }}_{n_k}^{+}{t}_i},i=0,1,\cdots ,N_1-1,k=0,1,\cdots ,\widehat{M}-1.$ (3.30)

Provided that $\Phi$ has full column rank, the independent spectral coefficient estimates are uniquely determined from the constrained least-squares problem

$\widehat{C}=\text{arg}\underset{C\in {ℂ}^{\widehat{M}}}{\text{min}}{\Vert y-2\mathrm{Re}\left(\Phi C\right)\Vert }_2^2.$ (3.31)

Consequently, the reconstructed zero-input boundary response is

$y_0^{*}\left(t\right)=2\mathrm{Re}\left({\displaystyle \sum }_{k=0}^{\widehat{M}-1}{\widehat{C}}_{n_k}^{+}{\text{e}}^{{\widehat{\lambda }}_{n_k}^{+}t}\right),t\in \left[T_1,T_3\right],$ (3.32)

which is real-valued by construction.

3.2.3. Step 3: Identification of the Wave Speed via Controlled Residual Fitting

A uniform sampling grid on $\left[T_2,T_3\right]$ is specified by

$t_i=T_2+i\text{Δ}{t}_2,i=0,1,\cdots ,N_2-1,$

where

$\text{Δ}{t}_2=\frac{T_3-T_2}{N_2-1}.$

Thus, $t_0=T_2$ and $t_{N_2-1}=T_3$ . The boundary control is chosen as

$U\left(t\right)=1,t\in \left[T_2,T_3\right].$

From (2.31) and (2.32), we obtain

$\begin{array}{cc}\begin{array}{c}y\left(t_i\right)=y_0\left(t_i\right)+y_1\left(t_i;c\right)\\ \approx y_0^{\text{*}}\left(t_i\right)+y_1\left(t_i;c\right),\end{array}& i=0,1,\cdots ,N_2-1.\end{array}$ (3.33)

We define a new temporal variable

${\tau }_i=t_i-T_2.$

Then

$\begin{array}{c}{y}_1\left(t_i;c\right)={\displaystyle {\int }_{T_2}^{t_i}G\left(t_i-s,0\right)\text{d}s}\\ ={\displaystyle {\int }_0^{{\tau }_i}{\displaystyle \sum }_{n=1}^{\infty }\frac{\text{i}{\omega }_n{c}^3}{{\alpha }_n^2}\left({\text{e}}^{{\lambda }_n^{-}\left({\tau }_i-s\right)}-{\text{e}}^{{\lambda }_n^{+}\left({\tau }_i-s\right)}\right)\text{d}s}\\ ={\displaystyle \sum }_{n=1}^{\infty }\frac{c^2}{{\alpha }_n^2}\left(2-{\text{e}}^{{\lambda }_n^{-}{\tau }_i}-{\text{e}}^{{\lambda }_n^{+}{\tau }_i}\right)\\ =2{\displaystyle \sum }_{n=1}^{\infty }{\beta }_n\left[1-\text{cos}\left({\omega }_{n}c{\tau }_i\right)\right].\end{array}$ (3.34)

Here, ${\left\{{\omega }_n\right\}}_{n\ge 1}$ are the spectral parameters determined by the Sturm-Liouville operator, and

${\beta }_n=\frac{1}{{\omega }_n\left({\omega }_n+\frac{1}{2}\text{sin}\left(2{\omega }_n\right)\right)}.$

On the controlled time interval $\left[T_2,T_3\right]$ , we define the residual signal by

$r\left(t_i\right):=y\left(t_i\right)-y_0^{\text{*}}\left(t_i\right),i=0,1,\cdots ,N_2-1.$ (3.35)

It follows from (3.33) that

$r\left(t_i\right)\approx y_1\left(t_i;c\right).$

For numerical evaluation, the infinite series in (3.34) is truncated after $Q$ terms:

$y_{1,Q}\left(t_i;c\right)=2{\displaystyle \sum }_{n=1}^Q{\beta }_n\left[1-\text{cos}\left({\omega }_{n}c{\tau }_i\right)\right].$ (3.36)

The wave-speed estimate is selected from the set of minimizers of the bounded nonlinear least-squares criterion

$c^{*}=\underset{c\in \left[c_{\min },c_{\max }\right]}{\arg \min }{J}_Q\left(c\right),$ (3.37)

where

$J_Q\left(c\right)={\displaystyle \sum }_{i=1}^{N_2-1}{\left|r\left(t_i\right)-y_{1,Q}\left(t_i;c\right)\right|}^2.$ (3.38)

A two-stage strategy is adopted. First, the objective function is evaluated over a coarse grid on $\left[c_{\text{min}},c_{\text{max}}\right]$ to identify a suitable initial candidate. A bounded scalar minimization method is then applied locally to obtain a refined estimate of $c^{\text{*}}$ .

3.2.4. Step 4: Reconstruction of the Initial States via Regularized Least Squares

Using the wave-speed estimate $c^{\text{*}}$ determined in Step 3, the initial state $\left(u_0,u_1\right)$ is identified from the zero-input boundary observations collected over $\left[T_1,T_2\right]$ .

Since $U\left(t\right)=0$ on $\left[T_1,T_2\right]$ , the boundary observation on this interval is generated only by the free evolution of the unknown initial states. By the truncated modal representation of the zero-input response, and replacing the unknown wave speed $c$ by its estimate $c^{\text{*}}$ , we approximate the observation data by

$y\left(t_i\right)\approx {\displaystyle \sum }_{n=1}^{P_i}\left[A_n\cos \left({\omega }_n{c}^{*}{t}_i\right)+B_n\sin \left({\omega }_n{c}^{*}{t}_i\right)\right],i=0,1,\cdots ,N_1-1.$ (3.39)

Here $P_i$ denotes the number of physical modes retained for the identification of the initial states, and $A_n,B_n$ are the real modal coefficients to be determined.

Let

$Y=\left[\begin{array}{c}y\left(t_0\right)\\ y\left(t_1\right)\\ ⋮\\ y\left(t_{N_1-1}\right)\end{array}\right]\in {ℝ}^{N_1},$ (3.40)

and define the unknown coefficient vector

$\Theta ={\left[\begin{array}{lllllll}{A}_1\hfill & B_1\hfill & A_2\hfill & B_2\hfill & \cdots \hfill & A_{P_i}\hfill & B_{P_i}\hfill \end{array}\right]}^{\top }\in {ℝ}^{2P_i}.$ (3.41)

We construct the matrix

$\Psi \in {ℝ}^{N_1\times 2P_i}$

by

${\Psi }_{i+1,2n-1}=\cos \left({\omega }_n{c}^{*}{t}_i\right),{\Psi }_{i+1,2n}=\sin \left({\omega }_n{c}^{*}{t}_i\right),$

for

$i=0,1,\cdots ,N_1-1,n=1,2,\cdots ,P_i.$

Then the identification of the coefficients $A_n$ and $B_n$ is reduced to the overdetermined linear system

$\Psi \Theta \approx Y.$ (3.42)

Equivalently, $\Theta$ is obtained by solving the linear least-squares problem

$\underset{\Theta \in {ℝ}^{2P_i}}{\min }{\Vert \Psi \Theta -Y\Vert }_2^2.$ (3.43)

Since the underlying initial-state identification problem is ill-posed, the discretized least-squares system may be severely ill-conditioned and sensitive to measurement noise and numerical errors. Therefore, we solve the above system by truncated singular value decomposition. Let

$\Psi =U\Sigma V^{\top }={\displaystyle \sum }_{k=1}^r{\sigma }_k{u}_k{v}_k^{\top },r:=\text{rank}\left(\Psi \right),$ (3.44)

where ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r>0$ are the nonzero singular values of Ψ. For a given truncation index $K$ , the TSVD regularized solution is

${\Theta }_K={\displaystyle \sum }_{k=1}^K\frac{u_k^{\top }Y}{{\sigma }_k}{v}_k.$ (3.45)

The truncation index is selected by the generalized cross-validation criterion

$K^{*}=\text{arg}\underset{1\le K\le r}{\text{min}}\text{G}CV\left(K\right),$ (3.46)

where

$\text{GCV}\left(K\right)=\frac{{\Vert \Psi {\Theta }_K-Y\Vert }_2^2}{{\left(N_1-K\right)}^2}.$ (3.47)

The final regularized coefficient vector is then defined as

${\Theta }^{*}:={\Theta }_{K^{*}}.$ (3.48)

Writing

${\Theta }^{*}={\left[\begin{array}{lllllll}{A}_1^{*}\hfill & B_1^{*}\hfill & A_2^{*}\hfill & B_2^{*}\hfill & \cdots \hfill & A_{P_i}^{*}\hfill & B_{P_i}^{*}\hfill \end{array}\right]}^{\top },$ (3.49)

we obtain the identified initial displacement and velocity as

$u_0^{*}\left(x\right):={\displaystyle \sum }_{n=1}^{P_i}{A}_n^{*}\cos \left({\omega }_{n}x\right)$ (3.50)

and

$u_1^{*}\left(x\right):={\displaystyle \sum }_{n=1}^{P_i}{\omega }_n{c}^{*}{B}_n^{*}\cos \left({\omega }_{n}x\right).$ (3.51)

Indeed, in the representation

$A_n\cos \left({\omega }_n{c}^{*}t\right)+B_n\sin \left({\omega }_n{c}^{*}t\right),$

the coefficient $A_n$ corresponds to the $n$ -th modal coefficient of the initial displacement, while ${\omega }_n{c}^{\text{*}}{B}_n$ corresponds to the $n$ -th modal coefficient of the initial velocity. Hence, the above formulas provide the finite-dimensional identifications of $u_0\left(x\right)$ and $u_1\left(x\right)$ .

4. Error Analysis

This section investigates the deterministic perturbation arising in the matrix-pencil identification of the spectral quantities in Step 1. The zero-input boundary observation is represented by an infinite modal expansion, whereas the matrix-pencil procedure operates on a finite exponential sequence. Consequently, truncating the modal expansion introduces a residual term into the sampled data and perturbs the associated shifted Hankel matrices. By combining the uniform truncation estimate established in Theorem 3.1 with perturbation results for Moore-Penrose inverses and matrix eigenvalues, we derive explicit bounds for the identified discrete poles and their corresponding continuous-time characteristic exponents.

Unlike the classical finite-dimensional exponential model with stochastic measurement noise, the perturbation considered here is generated by the neglected spectral tail of an infinite-dimensional boundary response. The analysis, therefore, quantifies how the finite-modal approximation propagates through the matrix-pencil identification.

Recall that the observation is sampled uniformly on $\left[T_1,T_2\right]$ at

$t_i=T_1+i\Delta t_1,i=0,1,\cdots ,N_1-1,\Delta t_1=\frac{T_2-T_1}{N_1-1}.$

The discrete observation model obtained in Step 1 is

$y_i={\displaystyle \sum }_{k=0}^{M-1}\left[R_k^{+}{\left(z_k^{+}\right)}^i+R_k^{-}{\left(z_k^{-}\right)}^i\right]+{\epsilon }_i,i=0,1,\cdots ,N_1-1,$ (4.1)

where

$R_k^{\pm }=C_{n_k}^{\pm }{\text{e}}^{{\lambda }_{n_k}^{\pm }{T}_1},z_k^{\pm }={\text{e}}^{{\lambda }_{n_k}^{\pm }\Delta t_1},{\epsilon }_i=r_M\left(t_i\right).$

Define the ideal finite exponential sequence by

$x_i:={\displaystyle \sum }_{k=0}^{M-1}\left[R_k^{+}{\left(z_k^{+}\right)}^i+R_k^{-}{\left(z_k^{-}\right)}^i\right].$ (4.2)

Then

$y_i=x_i+{\epsilon }_i.$

Because the initial state is real-valued, the two spectral branches occur in conjugate pairs. Provided that the sampled poles are mutually distinct, the retained $M$ physical modes generate

$P=2M$ (4.3)

nonzero discrete poles. The pencil parameter is assumed to satisfy

$P\le L\le N_1-P.$ (4.4)

Theorem 3.1 gives the uniform estimate

$\underset{0\le i\le N_1-1}{\text{max}}\left|{\epsilon }_i\right|\le \frac{C}{\sqrt{M}},$ (4.5)

where $C>0$ is the same constant as that appearing in Theorem 3.1 and is independent of $M$ .

Let $X_0$ and $X_1$ denote the shifted Hankel matrices formed from the ideal sequence ${\left\{x_i\right\}}_{i=0}^{N_1-1}$ , and let $Y_0$ and $Y_1$ be the corresponding matrices constructed from ${\left\{y_i\right\}}_{i=0}^{N_1-1}$ . Specifically,

$X_0=\left[\begin{array}{cccc}{x}_0& x_1& \cdots & x_{L-1}\\ x_1& x_2& \cdots & x_L\\ ⋮& ⋮& \ddots & ⋮\\ x_{N_1-L-1}& x_{N_1-L}& \cdots & x_{N_1-2}\end{array}\right],X_1=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_L\\ x_2& x_3& \cdots & x_{L+1}\\ ⋮& ⋮& \ddots & ⋮\\ x_{N_1-L}& x_{N_1-L+1}& \cdots & x_{N_1-1}\end{array}\right].$ (4.6)

and $Y_0,Y_1$ are defined analogously by replacing $x_i$ with $y_i$ . Hence,

$X_0,X_1,Y_0,Y_1\in {ℂ}^{\left(N_1-L\right)\times L}.$

Theorem 4.1 (Perturbation bounds for the shifted Hankel matrices) Let $L$ satisfy (4.4), and define

${\delta }_M:=C\sqrt{\frac{L\left(N_1-L\right)}{M}}.$ (4.7)

Then

${\Vert Y_0-X_0\Vert }_{\text{F}}\le {\delta }_M,{\Vert Y_1-X_1\Vert }_{\text{F}}\le {\delta }_M.$ (4.8)

Consequently,

${\Vert Y_j-X_j\Vert }_2\le {\delta }_M,j=0,1.$ (4.9)

Proof. Each entry of $Y_0-X_0$ is one of the residual samples ${\epsilon }_i$ . Since $Y_0-X_0$ contains $L\left(N_1-L\right)$ entries, (4.5) implies

$\begin{array}{c}{\Vert Y_0-X_0\Vert }_{\text{F}}^2={\displaystyle \sum }_{p=1}^{N_1-L}{\displaystyle \sum }_{q=1}^L{\left|{\left(Y_0-X_0\right)}_{p,q}\right|}^2\\ \le L\left(N_1-L\right)\underset{0\le i\le N_1-1}{\text{max}}{\left|{\epsilon }_i\right|}^2\\ \le \frac{C^{2}L\left(N_1-L\right)}{M}.\end{array}$

Taking square roots proves the first estimate in (4.8). The estimate for $Y_1-X_1$ follows from the same argument. Finally, the spectral norm is bounded above by the Frobenius norm, which gives (4.9).

Under the distinct-pole assumption and (4.4), the ideal Hankel matrix satisfies

$\text{rank}\left(X_0\right)=P.$ (4.10)

Let

${\sigma }_1\ge {\sigma }_2\ge \cdots$

be the singular values of $Y_0$ , and let $Y_{0,P}$ denote its rank-$P$ truncated SVD approximation. The analysis below assumes that the model-order selection in Step 1 identifies the correct number of discrete poles, namely,

$\widehat{P}=P.$

Define

${\tau }_P:={\Vert Y_0-Y_{0,P}\Vert }_2={\sigma }_{P+1},$ (4.11)

where ${\tau }_P=0$ if $P$ equals the number of available singular values. Introduce the relative perturbation quantity

$\rho :=\frac{{\tau }_P+{\delta }_M}{{\sigma }_P}.$ (4.12)

The ideal and perturbed matrix-pencil operators are defined by

$A_0:=X_0^{†}{X}_1,\widehat{A}:=Y_{0,P}^{†}{Y}_1,$ (4.13)

where $†$ denotes the Moore-Penrose inverse.

Suppose that

$Y_{0,P}=U_P{\Sigma }_P{V}_P^{\text{*}}.$ (4.14)

Then

$\widehat{A}=V_P{\Sigma }_P^{-1}{U}_P^{\text{*}}{Y}_1.$ (4.15)

The nonzero eigenvalues of $\widehat{A}$ coincide with those of the reduced matrix

$Z_{\text{E}}={\Sigma }_P^{-1}{U}_P^{\text{*}}{Y}_1{V}_P\in {ℂ}^{P\times P},$ (4.16)

which is precisely the matrix employed in Step 1. Indeed, this follows from the fact that the products $BC$ and $CB$ have identical nonzero eigenvalues whenever their dimensions are compatible.

Theorem 4.2 (Perturbation estimate for the recovered discrete poles) Assume that

$\rho <1$

and that the ideal matrix-pencil operator is diagonalizable:

$A_0=S\Lambda S^{-1}.$ (4.17)

Define

$\kappa \left(S\right):={\Vert S\Vert }_2{\Vert S^{-1}\Vert }_2$

and

${\eta }_A:=\frac{\phi \rho {\Vert Y_1\Vert }_2+{\delta }_M}{{\sigma }_P\left(1-\rho \right)},\phi :=\frac{1+\sqrt{5}}{2}.$ (4.18)

Then

${\Vert \widehat{A}-A_0\Vert }_2\le {\eta }_A.$ (4.19)

Furthermore, after a suitable labeling of the estimated poles,

$\left|{\widehat{z}}_j-z_j\right|\le {\eta }_z,j=1,2,\cdots ,P,$ (4.20)

where

${\eta }_z:=\kappa \left(S\right){\eta }_A.$ (4.21)

If

${\eta }_z<\frac{1}{2}\text{min}\left\{\underset{1\le j\le P}{\text{min}}\left|z_j\right|,\underset{\begin{array}{c}1\le j,k\le P\\ j\ne k\end{array}}{\text{min}}\left|z_j-z_k\right|\right\},$ (4.22)

then every identified nonzero pole can be associated uniquely with one true pole.

Proof. Using the triangle inequality, the difference between the perturbed and ideal matrix-pencil operators can be decomposed as

$\begin{array}{c}{\Vert \widehat{A}-A_0\Vert }_2={\Vert Y_{0,P}^{†}{Y}_1-X_0^{†}{X}_1\Vert }_2\\ \le {\Vert Y_{0,P}^{†}-X_0^{†}\Vert }_2{\Vert Y_1\Vert }_2+{\Vert X_0^{†}\Vert }_2{\Vert Y_1-X_1\Vert }_2.\end{array}$ (4.23)

By (4.11), Theorem 4.1, and the definition of $\rho$ ,

$\begin{array}{c}{\Vert Y_{0,P}-X_0\Vert }_2\le {\Vert Y_{0,P}-Y_0\Vert }_2+{\Vert Y_0-X_0\Vert }_2\\ \le {\tau }_P+{\delta }_M\\ =\rho {\sigma }_P.\end{array}$ (4.24)

Moreover,

${\Vert Y_{0,P}^{†}\Vert }_2=\frac{1}{{\sigma }_P}.$ (4.25)

Weyl’s singular-value perturbation inequality [33] yields

$\begin{array}{c}{\sigma }_P\left(X_0\right)\ge {\sigma }_P\left(Y_{0,P}\right)-{\Vert Y_{0,P}-X_0\Vert }_2\\ \ge \left(1-\rho \right){\sigma }_P.\end{array}$ (4.26)

Therefore,

${\Vert X_0^{†}\Vert }_2\le \frac{1}{\left(1-\rho \right){\sigma }_P}.$ (4.27)

Since

$\text{rank}\left(Y_{0,P}\right)=\text{rank}\left(X_0\right)=P,$

the same-rank perturbation estimate for Moore-Penrose inverses gives

$\begin{array}{c}{\Vert Y_{0,P}^{†}-X_0^{†}\Vert }_2\le \phi {\Vert Y_{0,P}^{†}\Vert }_2{\Vert X_0^{†}\Vert }_2{\Vert Y_{0,P}-X_0\Vert }_2\\ \le \frac{\phi \rho }{{\sigma }_P\left(1-\rho \right)}.\end{array}$ (4.28)

Substitution of (4.9), (4.27), and (4.28) into (4.23) yields

${\Vert \widehat{A}-A_0\Vert }_2\le \frac{\phi \rho {\Vert Y_1\Vert }_2+{\delta }_M}{{\sigma }_P\left(1-\rho \right)},$ (4.29)

which proves (4.19).

Applying the Bauer-Fike theorem [34] to (4.17), every eigenvalue $\widehat{z}$ of $\widehat{A}$ satisfies

$\text{dist}\left(\widehat{z},\text{spec}\left(A_0\right)\right)\le \kappa \left(S\right){\Vert \widehat{A}-A_0\Vert }_2.$ (4.30)

This proves (4.20) after a suitable labeling. Under condition (4.22), the perturbation disks centered at the distinct true poles are mutually disjoint and remain separated from the zero eigenvalue. Hence, the correspondence between the estimated and exact nonzero poles is unique.

Corollary 4.1 (Perturbation estimate for the characteristic exponents) Let

${\lambda }_j=\frac{\log z_j}{\Delta t_1},{\widehat{\lambda }}_j=\frac{\log {\widehat{z}}_j}{\Delta t_1}.$

Assume that $z_j$ and ${\widehat{z}}_j$ belong to the same analytic branch domain of the logarithm and that the line segment joining them does not intersect the origin or the selected branch cut. If

${\eta }_z<1,$

then

$\left|{\widehat{\lambda }}_j-{\lambda }_j\right|\le \frac{{\eta }_z}{\text{Δ}{t}_1\left(1-{\eta }_z\right)},j=1,2,\cdots ,P.$ (4.31)

Proof. The analytic logarithm satisfies

$\log {\widehat{z}}_j-\log z_j={\displaystyle {\int }_0^1\frac{{\widehat{z}}_j-z_j}{z_j+s\left({\widehat{z}}_j-z_j\right)}\text{d}s}.$ (4.32)

Since the characteristic exponents of the wave equation are purely imaginary,

${\lambda }_j\in \left\{{\lambda }_{n_k}^{+},{\lambda }_{n_k}^{-}\right\}=\left\{\text{i}{\omega }_{n_k}c,-\text{i}{\omega }_{n_k}c\right\},$

and therefore

$\left|z_j\right|=\left|{\text{e}}^{{\lambda }_j\text{Δ}{t}_1}\right|=1.$

For $s\in \left[0,1\right]$ ,

$\begin{array}{c}\left|z_j+s\left({\widehat{z}}_j-z_j\right)\right|\ge \left|z_j\right|-s\left|{\widehat{z}}_j-z_j\right|\\ \ge 1-{\eta }_z.\end{array}$ (4.33)

It follows from (4.32) that

$\left|\log {\widehat{z}}_j-\log z_j\right|\le \frac{{\eta }_z}{1-{\eta }_z}.$

Dividing by $\text{Δ}{t}_1$ proves (4.31).

Remark 4.1 The real-variable mean-value theorem used for positive real poles cannot be applied directly to the complex poles of the present wave equation. The integral representation (4.32) avoids this difficulty and explicitly accounts for the branch structure of the complex logarithm.

Remark 4.2 The condition $\rho <1$ ensures that the perturbation remains smaller than the weakest retained singular direction. As $\rho$ approaches one, the upper bound for ${\Vert X_0^{†}\Vert }_2$ becomes large, and the matrix-pencil recovery becomes increasingly sensitive to the truncation residual. The quantities ${\sigma }_P$ and $\kappa \left(S\right)$ further characterize the conditioning of the signal subspace and the eigenvalue problem, respectively.

Remark 4.3 For fixed $N_1$ and $L$ , the admissibility condition $2M\le L\le N_1-2M$ restricts the allowable truncation level. Accordingly, the estimate $C{M}^{-1/2}$ should be interpreted for feasible values of $M$ . Convergence of the complete discrete procedure as $M\to \infty$ would require the number of samples and the pencil dimension to increase simultaneously.

5. Numerical Simulation

Representative numerical experiments are presented to assess the proposed boundary switching-control framework for the joint identification of the wave speed and non-generic initial states in one-dimensional wave equations. The numerical study focuses on reconstruction accuracy, computational performance, and stability under measurement noise. All numerical computations are implemented on the MATLAB R2025b platform, and all calculated results are rounded to four decimal places for consistency. Considering that boundary measurement data in practical engineering are inevitably affected by sensor noise, environmental interference, and other factors, this study adopts multiplicative Gaussian white noise to simulate real measurement errors. The generation formula of noisy observation data is:

$y^{\delta }\left(t_i\right)=y\left(t_i\right)\left(1+\delta \cdot \xi \right),$ (5.1)

Here, $\delta$ specifies the relative noise amplitude, while $\xi$ is sampled from the uniform distribution on $\left[-1,1\right]$ . Three typical cases are considered in the numerical experiments to comprehensively evaluate the performance of the proposed algorithm: the noise-free case ($\delta =0$ ), the low-noise case ($\delta =0.005$ , i.e., 0.5% relative noise), and the high-noise case ($\delta =0.01$ , i.e., 1% relative noise). For each nonzero noise level, 50 independent Monte Carlo trials are performed, and the reported errors are given as the mean values with standard deviations.

The numerical investigation considers a one-dimensional wave equation subject to a controlled Neumann boundary condition at the left endpoint and a Robin boundary condition at the right endpoint. The reference wave speed is prescribed as $c^{\text{*}}=1$ , and the Robin coefficient at $x=1$ is fixed at $q=1$ . The theoretical eigenvalues of the first three odd-order modes of the system under the given boundary conditions are ${\omega }_1=0.8603$ , ${\omega }_3=6.4373$ and ${\omega }_5=12.6453$ , respectively. To verify the identification ability of the algorithm for non-generic initial states (i.e., initial states orthogonal to partial eigenfunctions of the system), the constructed initial displacement and initial velocity are strictly orthogonal to all even-order eigenfunctions, with the specific forms:

$\{\begin{array}{l}{u}_0^{*}\left(x\right)=2\cos \left({\omega }_{1}x\right)+4\cos \left({\omega }_{3}x\right)+3\cos \left({\omega }_{5}x\right),\\ u_1^{*}\left(x\right)=20\cos \left({\omega }_{1}x\right)+40\cos \left({\omega }_{3}x\right).\end{array}$ (5.2)

The synthetic boundary observation $y\left(t\right)=u\left(0,t;U,u_0^{\text{*}},u_1^{\text{*}}\right)$ is generated by solving the forward problem (1.1) with a leapfrog finite-difference scheme. The spatial interval $\left[0,1\right]$ is partitioned into 200 equal subintervals, yielding 201 grid points and a spatial mesh size of $\text{Δ}x=0.005$ . The temporal step used in the forward solver is $\text{Δ}t=0.0005$ . These discretization parameters satisfy the Courant-Friedrichs-Lewy stability requirement and therefore provide a stable numerical approximation of the forward solution.

For the inverse identification, the observation times are specified as $\left[T_1,T_2,T_3\right]=\left[0.1,4.1,8.1\right]$ , where $T_2$ denotes the control-switching instant. The resulting zero-input and controlled observation intervals, $\left[T_1,T_2\right]$ and $\left[T_2,T_3\right]$ , have the same duration. The boundary data on both intervals are sampled uniformly with $\text{Δ}{t}_1=\text{Δ}{t}_2=0.005$ . In all numerical examples, the sampling interval $\text{Δ}{t}_1$ is chosen sufficiently small to satisfy the no-aliasing condition (3.25), which guarantees that the identified continuous-time exponents are free of aliasing ambiguity. Because both endpoints of each interval are included, the corresponding numbers of samples are $N_1=N_2=801$ . For the matrix-pencil computation in Step 1, the pencil parameter is chosen as $L=N_1/3=267$ , which satisfies the admissibility requirement $P\le L\le N_1-P$ for the retained number $P$ of discrete poles.

Under the foregoing numerical settings, synthetic boundary measurements are generated by solving the forward problem with the finite-difference scheme and recording $y\left(t\right)=u\left(0,t\right)$ . The proposed identification procedure is then applied to the sampled data following Steps 1 - 4. For the noise-free data, the algorithm produces deterministic estimates of the wave speed and the initial states. For the noisy cases with $\delta =0.005$ and $\delta =0.01$ , the same procedure is repeated over 50 independent noise realizations, and the reported quantitative errors are given by the corresponding Monte Carlo means and standard deviations. The numerical results are presented below through the objective function for wave speed identification, the controlled residual fitting, the reconstructions of the initial displacement and velocity, and the associated error statistics. The numerical performance is illustrated in Figures 1-4, and the corresponding relative errors are summarized in Table 1.

Figure 1. Objective function $J\left(c\right)$ for wave speed identification.

Figure 1 displays the objective function $J\left(c\right)$ for the representative realizations under different noise levels. In all cases, the minimum is located close to the true value $c^{*}=1$ . This confirms that the controlled residual contains sufficient information to identify the wave speed.

Figure 2 compares the residual extracted from the observation data with the fitted controlled response. In the noise-free case, the two curves almost coincide. When noise is present, the residual data fluctuate around the fitted curve, while the main deterministic trend is still accurately captured. This demonstrates the validity of the controlled residual fitting used in Step 3.

Following the identification of the wave speed, the regularized least-squares method described in Step 4 is applied to the uncontrolled boundary observation to recover the initial displacement and velocity. The corresponding reconstruction results are presented in Figure 3 and Figure 4, respectively. The reconstructed curves are close to the exact initial states for all three noise levels. As expected, the reconstruction error increases as the noise level becomes larger.

Figure 2. Controlled residual fitting on $\left[T_2,T_3\right]$ .

Figure 3. Reconstruction of the initial displacement $u_0\left(x\right)$ under different noise levels.

Figure 4. Reconstruction of the initial velocity $u_1\left(x\right)$ under different noise levels.

The reconstruction accuracy is quantified by the following relative-error measures:

$E_c=\frac{\left|c^{*}-c\right|}{\left|c\right|},$

$E_{u_0}=\frac{{\Vert u_0^{*}-u_0\Vert }_{L^2\left(0,1\right)}}{{\Vert u_0\Vert }_{L^2\left(0,1\right)}},E_{u_1}=\frac{{\Vert u_1^{*}-u_1\Vert }_{L^2\left(0,1\right)}}{{\Vert u_1\Vert }_{L^2\left(0,1\right)}}.$

All numerical results are reported as either single representative realizations or Monte Carlo averages with standard deviations, depending on the noise level.

Table 1. Reconstruction errors under different noise levels.

Noise Level $\delta$	$c^{\text{*}}$	${10}^3{E}_c$	${10}^2{E}_{u_0}$	${10}^2{E}_{u_1}$
0	1.0016	1.57	2.82	3.68
0.005	1.0023 ± 0.0036	3.27 ± 2.73	6.47 ± 5.61	7.91 ± 6.63
0.01	1.001 ± 0.0074	5.96 ± 4.37	12.03 ± 8.70	14.35 ± 10.15

The quantitative results are summarized in Table 1. In the noise-free case, the remaining errors mainly come from the finite difference discretization, modal truncation, and regularization. The recovered wave speed is $c^{\text{*}}=1.0016$ , with relative error 1.5699 × 10⁻³. For noisy data, the mean relative error of the recovered wave speed increases from 3.27 × 10⁻³ at 0.5% noise to 5.9614 × 10⁻³ at 1% noise. The same increasing trend is observed for the reconstruction errors of $u_0$ and $u_1$ .

Several conclusions can be drawn from these numerical results. First, the wave speed is identified accurately in all three cases. Even for 1% multiplicative noise, the mean relative error of $c^{\text{*}}$ remains below 0.60%. Second, the reconstruction of the initial displacement is stable under moderate noise, with the mean relative error increasing from 2.8151% in the noise-free case to 6.4747% for 0.5% noise and 12.0280% for 1% noise. Third, the reconstruction of the initial velocity is slightly less accurate, with the mean relative error increasing from 3.6767% in the noise-free case to 7.9095% for 0.5% noise and 14.3520% for 1% noise. This is consistent with the fact that velocity reconstruction is more sensitive to noise and high-frequency components.

Overall, the simulations show that the proposed method can simultaneously identify the wave speed and reconstruct non-generic initial states from finite-time boundary observations. The numerical results also demonstrate the stability of the method under moderate measurement noise.

6. Conclusions

This study addressed the simultaneous recovery of the wave speed and non-generic initial states for a one-dimensional wave equation from boundary measurements collected over a finite time horizon. A boundary switching strategy was introduced to separate the observation process into a zero-input stage and a controlled stage. During the zero-input stage, the matrix-pencil method was employed to recover the relevant spectral information from the boundary response. The wave speed was subsequently determined by fitting the controlled residual through a bounded one-dimensional nonlinear least-squares problem. Once the wave speed had been estimated, the initial displacement and velocity were reconstructed from the zero-input measurements using a regularized least-squares formulation.

The proposed method separates the identification of the wave speed from the reconstruction of the initial states, which improves the stability of the inverse procedure and makes the approach suitable for non-generic initial data. Numerical experiments based on finite difference-generated observations show that the method can accurately recover the wave speed and reconstruct both initial states. The Monte Carlo results under multiplicative noise further demonstrate the robustness of the algorithm, with reconstruction errors increasing consistently as the noise level grows.

Further research will investigate extensions of the present framework to wave equations with spatially varying coefficients, broader classes of boundary conditions, and multidimensional configurations.