Spectral Analysis and Numerical Investigation of a Flexible Structure with Nonconservative Boundary Data

Abstract

Analytic and numerical results of the Euler-Bernoulli beam model with a two-parameter family of boundary conditions have been presented. The co-diagonal matrix depending on two control parameters (k1 and k2) relates a two-dimensional input vector (the shear and the moment at the right end) and the observation vector (the time derivatives of displacement and the slope at the right end). The following results are contained in the paper. First, high accuracy numerical approximations for the eigenvalues of the discretized differential operator (the dynamics generator of the model) have been obtained. Second, the formula for the number of the deadbeat modes has been derived for the case when one control parameter, k1, is positive and another one, k2, is zero. It has been shown that the number of the deadbeat modes tends to infinity, as k1→1+ and k2=0. Third, the existence of double deadbeat modes and the asymptotic formula for such modes have been proven. Fourth, numerical results corroborating all analytic findings have been produced by using Chebyshev polynomial approximations for the continuous problem.

1. Introduction

The present paper is concerned with the spectral analysis and numerical investigation of the eigenvalues of the Euler-Bernoulli beam model. The beam is clamped at the left end and subject to linear feedback-type conditions with a non-dissipative feedback matrix [1, 2]. Depending on the boundary parameters $k_1$ and $k_2$, the model can be either conservative, dissipative, or completely non-dissipative. We focus on the non-dissipative case, i.e., when the energy of a vibrating system is not a decreasing (or nonincreasing) function of time. In our approach, the initial-boundary value problem describing the beam dynamics is reduced to the first order in time evolution equation in the state Hilbert space $\mathcal{H}$. The evolution of the system is completely determined by the dynamics generator ${\mathcal{L}}_{k_1,k_2}$, which is an unbounded non-self-adjoint matrix differential operator (see Eqs. (2), (3), and (8)). The eigenmodes and the mode shapes of the flexible structure are defined as the eigenvalues (up to a multiple $i$) and the generalized eigenvectors of ${\mathcal{L}}_{k_1,k_2}$.

Based on the results of [1, 2], the dynamics generator has a purely discrete spectrum, whose location on the complex plane is determined by the controls $k_1$ and $k_2$. Having in mind the practical applications of the asymptotic formulas [3, 4, 5], we discuss the case of $k_1\ge 0$ and $k_2\ge 0$, such that $\mid k_1\mid +\mid k_2\mid >0$ (see Proposition 2). As shown in [2], even though the operator ${\mathcal{L}}_{k_1,k_2}$ is non-dissipative, for the case $k_1>0$ and $k_2=0$ (or $k_1=0$ and $k_2>0$), the entire set of eigenvalues is located in the closed upper half of the complex plane $\mathbb{C}$, which means that all eigenmodes are stable or neutrally stable. (We recall that to obtain an elastic mode from an eigenvalues of ${\mathcal{L}}_{k_1,k_2}$, one should multiply the eigenvalue by a factor $i$).

In the paper we address the question of accuracy of the asymptotic formulas for the eigenvalues. Namely, under what conditions the leading asymptotic terms informulas (20)and(21)can be used for practical estimation of the actual frequencies of the flexible beam? Numerical simulations show that the accuracy of the asymptotic formulas is really high; the leading asymptotic terms can be used by practitioners almost immediately, i.e., almost from the first vibrational mode. The second question is concerned with the role of the deadbeat modes. A deadbeat mode is a purely negative elastic mode that generates a solution of the evolution equation exponentially decaying in time. The deadbeat modes are important in engineering applications. As we prove in the paper, when the boundary parameter $k_1$ is close to 1 (while $k_2=0$), the number of the deadbeat modes is so large that the corresponding mode shapes become important for the description of the beam dynamics. More precisely, the number of deadbeat modes tends to infinity as $k_1\to 1^{+}$.

We have also shown that there exists a sequence of values of the parameter $k_1$, i.e., ${\left\{k_1^{\left(n\right)}\right\}}_{n=1}^{\infty }$, such that for each $k_1=k_1^{\left(n\right)}$ there exist a finite number of deadbeat modes and each corresponds to a double eigenvalue of the dynamics generator ${\mathcal{L}}_{k_1,k_2}$. For each value $k_1^{\left(n\right)}$, the operator ${\mathcal{L}}_{k_1,k_2}$ has a two-dimensional root subspace spanned by an eigenvector and an associate vector. This result means that for a double deadbeat mode (corresponding to $k_1^{\left(n\right)}$), there exists a mode shape and an associate mode shape. This fact indicates that for some values of $k_1$ and $k_2$, there exists a significant number of associate vectors of ${\mathcal{L}}_{k_1,k_2}$. Therefore, if one can prove that the set of the generalized eigenvectors (eigenvectors and associate vectors together) forms an unconditional basis for the state space, then construction of the bi-orthogonal basis [6] would be a more complicated problem than for the case when no associate vectors exist.

Finally, we mention that the feedback control of beams is a well-studied area [6], with multiple applications to the control of robotic manipulators, long and slender aircraft wings, propeller blades, large space structure [7, 8], and the dynamics of carbon nanotubes [9]. The analysis of a classical beam model with nonstandard feedback control law that originated in engineering literature [4, 10, 11, 12] may be of interest for both analysts and practitioners.

This paper is organized as follows. In Section 1 we formulate the initial-boundary value problem for the Euler-Bernoulli beam model. In Section 2, we reformulate the problem as an evolution equation in the Hilbert space of Cauchy data (the energy space). The dynamics generator ${\mathcal{L}}_{k_1,k_2}$, which is a non-self-adjoint matrix differential operator depending on two parameters,$k_1$and$k_2$, is the main object of interest. The eigenvalues and the generalized eigenvectors of ${\mathcal{L}}_{k_1,k_2}$ correspond to the modes and the mode shapes of the beam. We also give numerical approximations and graphical representations of the eigenvalues of a discrete approximation of the main operator (see Tables 1 and 2 and Figures 1 and 2). In Section 3, we study the deadbeat modes and derive the estimates for the number of the deadbeat modes from below and above for different values of the boundary parameters (see Figure 5). Section 4 is concerned with the asymptotic approximation for the set of double deadbeat modes (see Tables 3 and 4 and Figures 6 and 7). In Section 5, we outline the numerical scheme used for the spectral analysis of the finite-dimensional approximation of the dynamics generator.

$k_1=0$, $k_2=0.5$, $\mathit{EI}=1$, $\mathit{\rho }=0.1$, $N=64$, ${\mathit{\epsilon }}_f={10}^{-20}$
No.	Numerical	Analytic	No.	Numerical	Analytic
1.	$-7988.1+237.00i$	$-7988.1+237.00i$	18.	$28.860+13.515i$	$29.453+14.812i$
2.	$-7020.6+222.19i$	$-7020.6+222.19i$	19.	$123.16+29.723i$	$123.08+29.625i$
3.	$-6115.5+207.37i$	$-6115.5+207.37i$	20.	$279.13+44.431i$	$279.14+44.437i$
4.	$-5272.8+192.56i$	$-5272.8+192.56i$	21.	$497.61+59.250i$	$497.61+59.250i$
5.	$-4492.5+177.75i$	$-4492.5+177.75i$	22.	$778.50+74.062i$	$778.50+74.062i$
6.	$-3774.7+162.94i$	$-3774.7+162.94i$	23.	$1121.8+88.875i$	$1121.8+88.875i$
7.	$-3119.3+148.12i$	$-3119.3+148.12i$	24.	$1527.6+103.69i$	$1527.6+103.69i$
8.	$-2526.3+133.31i$	$-2526.3+133.31i$	25.	$1995.7+118.50i$	$1995.7+118.50i$
9.	$-1995.7+118.50i$	$-1995.7+118.50i$	26.	$2526.3+133.31i$	$2526.3+133.31i$
10.	$-1527.6+103.69i$	$-1527.6+103.69i$	27.	$3119.3+148.12i$	$3119.3+148.12i$
11.	$-1121.8+88.875i$	$-1121.8+88.875i$	28.	$3774.7+162.94i$	$3774.7+162.94i$
12.	$-778.5+74.062i$	$-778.5+74.062i$	29.	$4492.5+177.75i$	$4492.5+177.75i$
13.	$-497.61+59.250i$	$-497.61+59.250i$	30.	$5272.8+192.56i$	$5272.8+192.56i$
14.	$-279.13+44.431i$	$-279.14+44.437i$	31.	$6115.5+207.37i$	$6115.5+207.37i$
15.	$-123.16+29.723i$	$-123.08+29.625i$	32.	$7020.6+222.19i$	$7020.6+222.19i$
16.	$-28.860+13.515i$	$-29.453+14.812i$	33.	$7988.1+237.00i$	$7988.1+237.00i$
17.	$-2.2\cdot {10}^{-17}+4.6007i$

Table 1.

Approximations of the eigenvalues for the discrete and “continuous” operators ($\mathcal{K}>1$).

$k_1=1.3$, $k_2=1.2$, $\mathit{EI}=10$, $\mathit{\rho }=0.1$, $N=64$, ${\mathit{\epsilon }}_f={10}^{-20}$
No.	Numerical	Analytic	No.	Numerical	Analytic
1.	$-25266-229.07i$	$-25266-229.07i$	18.	$99.467-19.816i$	$98.177-14.317i$
2.	$-22206-214.76i$	$-22206-214.76i$	19.	$394.17-28.149i$	$394.26-28.634i$
3.	$-19344-200.44i$	$-19344-200.44i$	20.	$887.75-42.983i$	$887.75-42.951i$
4.	$-16679-186.12i$	$-16679-186.12i$	21.	$1578.6-57.267i$	$1578.6-57.268i$
5.	$-14212-171.81i$	$-14212-171.81i$	22.	$2466.9-71.586i$	$2466.9-71.585i$
6.	$-11942-157.49i$	$-11942-157.49i$	23.	$3552.5-85.902i$	$3552.5-85.903i$
7.	$-9869.1-143.17i$	$-9869.1-143.17i$	24.	$4835.6-100.22i$	$4835.6-100.22i$
8.	$-7993.9-128.85i$	$-7993.9-128.85i$	25.	$6316.0-114.54i$	$6316.0-114.54i$
9.	$-6316.0-114.54i$	$-6316.0-114.54i$	26.	$7993.9-128.85i$	$7993.9-128.85i$
10.	$-4835.6-100.22i$	$-4835.6-100.22i$	27.	$9869.1-143.17i$	$9869.1-143.17i$
11.	$-3552.5-85.902i$	$-3552.5-85.903i$	28.	$11942-157.49i$	$11942-157.49i$
12.	$-2466.9-71.586i$	$-2466.9-71.585i$	29.	$14212-171.81i$	$14212-171.81i$
13.	$-1578.6-57.267i$	$-1578.6-57.268i$	30.	$16679-186.12i$	$16679-186.12i$
14.	$-887.75-42.983i$	$-887.75-42.951i$	31.	$19344-200.44i$	$19344-200.44i$
15.	$-394.17-28.149i$	$-394.26-28.634i$	32.	$22206-214.76i$	$22206-214.76i$
16.	$-99.467-19.816i$	$-98.177-14.317i$	33.	$25266-229.07i$	$25266-229.07i$
17.	$1.5\cdot {10}^{-17}+7.7256i$

Table 2.

Approximations of the eigenvalues for the discrete and “continuous” operators ($\mathcal{K}<-1$).

Figure 1.

Graphical representation of the eigenvalues of the discrete and “continuous” operators ($\mathcal{K}>1$).

Figure 2.

Graphical representation of the eigenvalues of the discrete and “continuous” operators ($\mathcal{K}<-1$).

$k_2=0$, $\mathit{EI}=1$, $\mathit{\rho }=1$, $N=64$, ${\mathit{\epsilon }}_f={10}^{-30}$
No.	$k_1=1+{10}^{-4}$	$k_1=1+{10}^{-7}$	$k_1=1+{10}^{-10}$
1.	$-222.22+155.56i$	$-176.06+264.07i$	$-106.90+372.41i$
2.	$-133.38+124.45i$	$-87.723+211.27i$	$-2.4378\cdot {10}^{-16}+254.44i$
3.	$-64.540+93.396i$	$-3.9609\cdot {10}^{-18}+162.37i$	$-3.0717\cdot {10}^{-17}+123.66i$
4.	$-9.8081+58.559i$	$-5.7977\cdot {10}^{-19}+116.23i$	$-1.3012\cdot {10}^{-17}+4.9349i$
5.	$-7.7725\cdot {10}^{-21}+5.0488i$	$-6.3661\cdot {10}^{-20}+44.182i$	$-8.9232\cdot {10}^{-18}+44.421i$
6.	$4.9365\cdot {10}^{-21}+38.994i$	$-6.0066\cdot {10}^{-20}+4.9383i$	$8.9143\cdot {10}^{-18}+44.406i$
7.	$7.8007\cdot {10}^{-21}+4.8257i$	$6.0114\cdot {10}^{-20}+4.9313i$	$1.3012\cdot {10}^{-17}+4.9347i$
8.	$9.8081+58.559i$	$6.6801\cdot {10}^{-20}+44.651i$	$2.9907\cdot {10}^{-17}+123.09i$
9.	$64.540+93.396i$	$2.9381\cdot {10}^{-18}+139.20i$	$1.1286\cdot {10}^{-16}+234.06i$
10.	$133.38+124.45i$	$87.723+211.27i$	$3.2280\cdot {10}^{-16}+315.13i$
11.	$222.22+155.56i$	$176.06+264.07i$	$106.90+372.41i$

Table 3.

Eigenvalues closest to the imaginary axis as $k_1\to 1^{+}$.

$k_2=0$, $\mathit{EI}=1$, $\mathit{\rho }=1$, $N=64$, ${\mathit{\epsilon }}_f={10}^{-30}$
#	$k_1=1-{10}^{-4}$	$k_1=1-{10}^{-7}$	$k_1=1-{10}^{-10}$
1.	$-175.34+140.01i$	$-129.20+237.56i$	$-60.143+338.80i$
2.	$-96.394+108.84i$	$-50.206+186.90i$	$-8.6602+240.01i$
3.	$-36.896+78.778i$	$-8.0431+121.15i$	$-0.28608+123.37i$
4.	$-6.0769+42.171i$	$-0.23455+44.410i$	$-0.0074192+44.413i$
5.	$-0.11141+4.9324i$	$-0.0035253+4.9348i$	$-0.00011148+4.9348i$
6.	$0.11141+4.9324i$	$0.0035253+4.9348i$	$0.00011148+4.9348i$
7.	$6.0769+42.171i$	$0.23455+44.410i$	$0.0074192+44.413i$
8.	$36.896+78.778i$	$8.0431+121.15i$	$0.28608+123.37i$
9.	$96.394+108.84i$	$50.206+186.90i$	$8.6602+240.01i$
10.	$175.34+140.01i$	$129.20+237.56i$	$60.143+338.80i$

Table 4.

Eigenvalues closest to the imaginary axis as $k_1\to 1^{-}$.

1.1 The initial-boundary value problem for the Euler-Bernoulli beam model of a unit length

The Lagrangian of the system is defined by [10, 11]

$$\frac{1}{2}{\int }_0^1\left[\mathit{\varrho }\left(x\right)A\left(x\right)h_t^2\left(x,t\right)-E\left(x\right)I\left(x\right)h_{\mathit{xx}}^2\left(x,t\right)\right]\mathit{dx},$$E1

where $h\left(x,t\right)$ is the transverse deflection, $E\left(x\right)$ is the modulus of elasticity, $I\left(x\right)$ is the area moment of inertia, $\mathit{\varrho }\left(x\right)$ is the linear density, and $A\left(x\right)$ is the cross-sectional area of the beam.

Assuming that the beam is clamped at the left end $\left(x=0\right)$ and free at the right end $\left(x=1\right)$, and applying Hamilton’s variational principle to the action functional defined by (1), we obtain the equation of motion

$$\mathit{\varrho }\left(x\right)A\left(x\right)h_{\mathit{tt}}\left(x,t\right)+{\left(E\left(x\right)I\left(x\right)h_{\mathit{xx}}\left(x,t\right)\right)}_{\mathit{xx}}=0,0\le x\le 1,t>0,$$E2

and the boundary conditions

$$h\left(0,t\right)=h_x\left(0,t\right)=0\text{and}M\left(1,t\right)=Q\left(1,t\right)=0,$$E3

where $M\left(x,t\right)$ and $Q\left(x,t\right)$ are the moment and the shear, respectively [10]:

$$M\left(x,t\right)=E\left(x\right)I\left(x\right)h_{\mathit{xx}}\left(x,t\right)\text{and}Q\left(x,t\right)=M_x\left(x,t\right).$$E4

Now we replace the free right-end conditions from Eq. (3) with the following boundary feedback control law [2, 4]. Define the input and the output as

$$U\left(t\right)={\left[-Q\left(1,t\right),M\left(1,t\right)\right]}^T\text{and}Y\left(t\right)={\left[h_t\left(1,t\right),h_{\mathit{xt}}\left(1,t\right)\right]}^T,$$E5

where $T$ stands for transposition. The feedback control law is given by

$$U\left(t\right)=\mathit{KY}\left(t\right),$$E6

where $K$ is the $2\times 2$ feedback matrix. We select

$$K=\text{codiag}\left(-k_2,-k_1\right),k_1,k_2\ge 0,$$E7

with $k_1,k_2$ being the control parameters. The feedback (6) can be written as

$${\left.E\left(1\right)I\left(1\right)h_{\mathit{xx}}(1,t)=-k_1{h}_t(1,t\left)\text{and}\right(E\left(x\right)I\left(x\right)h_{\mathit{xx}}(x,t\left)\right){}_x\right|}_{x=1}=k_2{h}_{\mathit{xt}}\left(1,t\right).$$E8

Finally, we arrive to the following initial-boundary value problem: the equation of motion (2), the boundary conditions (3), and the standard initial conditions $h\left(x,0\right)=h_0\left(x\right),h_t\left(x,0\right)=h_1\left(x\right).$

Notice that the choice of a feedback matrix $K$ defines whether the system is dissipative or not. Indeed, let $\mathcal{E}\left(t\right)$ be the energy of the system, defined by representation (1). Evaluating ${\mathcal{E}}_t\left(t\right)$ on the solutions of Eq. (2) satisfying the left-end conditions from Eqs. (3), we obtain

$$\begin{array}{cc}{\mathcal{E}}_t\left(t\right)& ={\int }_0^1\left[\mathit{\varrho }\left(x\right)A\left(x\right)h_t\left(x,t\right)h_{\mathit{tt}}\left(x,t\right)+E\left(x\right)I\left(x\right)h_{\mathit{xx}}\left(x,t\right)h_{\mathit{xxt}}\left(x,t\right)\right]x\\ & =-\left(E\left(x\right)I\left(x\right)h_{\mathit{xx}}\right(xt\left)\right){}_x{h}_t(xt){|}_{x=1}+E\left(1\right)I\left(1\right)h_{\mathit{xx}}\left(1,t\right)h_{\mathit{xt}}\left(1,t\right).\end{array}$$E9

Taking into account Eqs. (4) and (6), we represent the right-hand side of Eq. (9) as the dot product in ${\mathbb{R}}^2$:

$${\mathcal{E}}_t\left(t\right)=-Q\left(1,t\right)h_t\left(1,t\right)+M\left(1,t\right)h_{\mathit{xt}}\left(1,t\right)=U\left(t\right)\cdot Y\left(t\right)=\mathit{KY}\left(t\right)\cdot Y\left(t\right).$$E10

With the choice of $K$ as in Eq. (7), we have

$${\mathcal{E}}_t\left(t\right)=\left[\begin{array}{cc}0& -k_2\\ -k_1& 0\end{array}\right]\left[\begin{array}{l}2h_t\left(1,t\right)\\ h_{\mathit{xt}}\left(1,t\right)\end{array}\right]\cdot \left[\begin{array}{l}2h_t\left(1,t\right)\\ h_{\mathit{xt}}\left(1,t\right)\end{array}\right]=-\left(k_1+k_2\right)h_t\left(1,t\right)h_{\mathit{xt}}\left(1,t\right).$$E11

Thus the system is not dissipative for all nonnegative values of $k_1$ and $k_2$.

2. Operator form of the problem

In what follows, we incorporate the cross-sectional area $A\left(x\right)$ into the density, write $\rho \left(x\right)$ instead of $\mathit{\varrho }\left(x\right)A\left(x\right)$, and also abbreviate $\mathit{EI}\left(x\right)\equiv E\left(x\right)I\left(x\right)$. Let $\mathcal{H}$ be the Hilbert space of two-component vector functions $U\left(x\right)={\left[u_0\left(x\right),u_1\left(x\right)\right]}^T$ equipped with the following norm:

$${‖U‖}_{\mathcal{H}}^2=\frac{1}{2}{\int }_0^1\left[\mathit{EI}\left(x\right){\left|u_0″\left(x\right)\right|}^2+\rho \left(x\right){\left|u_1\left(x\right)\right|}^2\right]x.$$E12

Assuming that $\mathit{EI},\rho \in {\mathcal{C}}^2\left(0,1\right)$ are positive functions, we obtain that the closure of smooth functions $U\left(x\right)={\left[u_0\left(x\right),u_1\left(x\right)\right]}^T$ satisfying $u_0\left(0\right)=u_0^{\prime }\left(0\right)=0$ will produce the energy space $\mathcal{H}=H_0^2\left(0,1\right)\times L^2\left(0,1\right)$. Here $H_0^2\left(0,1\right)=\left\{u\in H^2\left(0,1\right):u\left(0\right)=u^{\prime }\left(0\right)=0\right\}$, and the equality of function spaces is understood in the sense of a Hilbert-space isomorphism.

Problem (2) with conditions (3) can be represented as the time evolution problem:

$$U_t\left(x,t\right)=i\left({\mathcal{L}}_{k_1,k_2}U\right)\left(x,t\right)\text{and}U\left(x,0\right)={\left[u_0\left(x\right),u_1\left(x\right)\right]}^T,$$E13

where $0\ge x\ge 1,t\ge 0$. The dynamics generator ${\mathcal{L}}_{k_1,k_2}$ is given by the following matrix differential expression:

$${\mathcal{L}}_{k_1,k_2}=-i\left[\begin{array}{cc}0& 1\\ -\frac{1}{\rho \left(x\right)}\frac{{\partial }^2}{\partial x^2}\left(\mathit{EI}\left(x\right)\frac{{\partial }^2}{\partial x^2}\right)& 0\end{array}\right],$$E14

defined on the domain

$$\begin{array}{l}\mathcal{D}\left({\mathcal{L}}_{k_1,k_2}\right)=\{U={\left(u_0,u_1\right)}^T\in \mathcal{H}:u_0\in H^4\left(0,1\right),u_1\in H_0^2\left(0,1\right)\mathrm{;}{u}_1\left(0\right)=u_1^{\prime }\left(0\right)=0\mathrm{;}\\ {\left.\mathit{EI}\left(1\right)u_0″\left(1\right)=-k_1{u}_1\left(1\right),{\left(\mathit{EI}\left(x\right)u_0″\left(x\right)\right)}^{\prime }\right|}_{x=1}=k_2{u}_1\prime \left(1\right)\}.\end{array}$$E15

For any $\left(k_1,k_2\right)\in {\mathbb{R}}^2$, the adjoint operator ${\mathcal{L}}_{k_1,k_2}^{\ast }$ [13] is given by

$${\mathcal{L}}_{k_1,k_2}^{\ast }={\mathcal{L}}_{-k_2,-k_1},$$E16

i.e., ${\mathcal{L}}_{k_1,k_2}^{\ast }$ is defined by the same differential expression (14) on the domain described in Eq. (15), where $k_1$ and $k_2$ are replaced by $\left(-k_2\right)$ and $\left(-k_1\right)$, respectively. It follows from Eq. (16) that ${\mathcal{L}}_{0,0}$ is self-adjoint in $\mathcal{H}$ and thus ${\mathcal{L}}_{0,0}$ is the dynamics generator of the clamped-free beam model. For the reader’s convenience, we summarize the properties of ${\mathcal{L}}_{k_1,k_2}$ from [1, 2] needed for the present work.

Proposition 1:

1. ${\mathcal{L}}_{k_1,k_2}$ is an unbounded operator with compact resolvent, whose spectrum consists of a countable set of normal eigenvalues (i.e., isolated eigenvalues, each of finite algebraic multiplicity [6, 13]).

2. For each $\left(k_1,k_2\right)\in {\mathbb{R}}^2,\mid k_1\mid +\mid k_2\mid >0$, the operator ${\mathcal{L}}_{k_1,k_2}$ is a rank-two perturbation of the self-adjoint operator ${\mathcal{L}}_{0,0}$ in the sense that the operators ${\mathcal{L}}_{k_1,k_2}^{-1}$ and ${\mathcal{L}}_{0,0}^{-1}$ exist and are related by the rule

$${\mathcal{L}}_{k_1,k_2}^{-1}={\mathcal{L}}_{0,0}^{-1}+{\mathcal{T}}_{k_1,k_2},$$E17

where ${\mathcal{T}}_{k_1,k_2}$ is a rank-two operator. A similar decomposition is valid for the adjoint operator, i.e.,

$${\left({\mathcal{L}}_{k_1,k_2}^{-1}\right)}^{\ast }={\mathcal{L}}_{0,0}^{-1}+{\mathcal{T}}_{k_1,k_2}^{\ast },{\mathcal{T}}_{k_1,k_2}^{\ast }={\mathcal{T}}_{-k_2,-k_1}.$$E18

From now on, we assume that the structural parameters are constant. In the case of variable parameters, the spectral asymptotics will have the same leading terms and remainder terms depending on parameter smoothness.

Proposition 2: Assume that ${\mathrm{k}}_1,{\mathrm{k}}_2>0$ and ${\mathrm{k}}_1{\mathrm{k}}_2\ne \mathrm{EI}\mathrm{\rho }$. Let

$$\mathcal{K}=\frac{k_1+k_2}{A-k_1{k}_2/A},A=\sqrt{\mathit{EI}\rho },\text{and}\mid \mathcal{K}\mid \ne 1.$$E19

The following asymptotic approximations for the eigenvalues ${\lambda }_n$ (as $\mid n\mid \to \infty$) of the operator ${\mathcal{L}}_{k_1,k_2}$ hold:

1. If $1<\mid \mathcal{K}\mid <\infty$, then for $\mid n\mid \to \infty$ one has

   $${\lambda }_n=sign\left(\mathcal{K}n\right)\sqrt{\frac{\mathit{EI}}{\rho }}\left[{\left(\mathit{\pi n}\right)}^2-\frac{1}{4}{ln}^2\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)+\mathit{i\pi n}ln\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)\right]+O\left(n{\mathrm{e}}^{-\pi \mid n\mid }\right).$$E20

2. If $0<\mathcal{K}<1$, then for $n\to \infty$ one has

$$\begin{array}{l}{\lambda }_n=sign\left(\mathcal{K}n\right)\sqrt{\frac{\mathit{EI}}{\rho }}\left[{\left(\frac{2n+1}{2}\pi \right)}^2-\frac{1}{4}{ln}^2\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)\right.+\\ \left.\mathit{i\pi }\left(\frac{2n+1}{2}\right)ln\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)\right]+O\left(n{\mathrm{e}}^{-\pi \mid n\mid }\right).\end{array}$$E21

First of all, we address the question of accuracy of the asymptotic formulas (20) and (21). By its nature, formula (20) (as well as formula (21)) means that for any small $\epsilon >0$, one can find a positive integer $N$, such that all eigenvalues ${\lambda }_n$ with $\mid n\mid \ge N+1$ satisfy the estimate

$$\left|{\lambda }_n-sign\left(\mathcal{K}n\right)\sqrt{\frac{\mathit{EI}}{\rho }}\left[{\left(\mathit{\pi n}\right)}^2-\frac{1}{4}{ln}^2\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)+\mathit{i\pi n}ln\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)\right]\right|\le \epsilon$$E22

for the case when $1<\mid \mathcal{K}\mid <\infty$ and

$$\left|{\lambda }_n-sign\left(\mathcal{K}n\right)\sqrt{\frac{\mathit{EI}}{\rho }}\left[{\left(\frac{2n+1}{2}\pi \right)}^2-\frac{1}{4}{ln}^2\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)+\mathit{i\pi }\left(\frac{2n+1}{2}\right)ln\left(\frac{\mathcal{K}+1}{\mathcal{K}-1}\right)\right]\right|\le \epsilon$$E23

for the case when $0<\mid \mathcal{K}\mid <1$. The following important question holds: From which index$N$can the eigenvalues be approximated by the leading asymptotic terms with acceptable accuracy? In other words, can one claim that the asymptotic formulas (20) and (21) are valuable to practitioners, or are they just important mathematical results of the spectral analysis?

The results of numerical simulations (see Tables 1 and 2 and Figures 1 and 2) show that the asymptotic formulas are indeed quite accurate. That is, if one places on the complex plane the numerically produced sets of the eigenvalues, then the theoretically predicted distribution of eigenvalues can be seen almost immediately. To obtain these results, we used the numerical procedure based on Chebyshev polynomial approximations [14, 15, 16], as outlined in Section 5.

In Figures 1 and 2, we represent the graphical distribution of the eigenvalues corresponding to the discretized operator (“numerical” eigenvalues) and the leading asymptotic terms from Eqs. (20) and (21) (“analytic” eigenvalues). In Tables 1 and 2, the numerical values of the corresponding graphical points on Figures 1 and 2 are listed. We have used the following notations: $N=64$ is the number of grid points on $\left[0,1\right]$, and ${\epsilon }_f$ is the filtering parameter as described in Eq. (69). It can be easily seen from the graphs and tables that the two sets of data coincide almost immediately, i.e., the leading asymptotic terms in the approximations are very close to the numerically approximated eigenvalues.

Figure 3 shows the sub-domains of the $k_1,k_2$-plane, which correspond to different intervals for the values of $\mathcal{K}$ defined by Eq. (19). On the sub-domain with dark gray color $\mathcal{K}$ such that $\mid \mathcal{K}\mid >1$, i.e., to evaluate the asymptotic approximation for the eigenvalues, one needs formula (20), while on the complementary sub-domain, one needs formula (21).

Figure 3.

Regions of $\mathcal{K}$ on the $k_1,k_2$-plane, $A=\sqrt{\mathit{EI}\rho }=1$.

3. The deadbeat modes

An eigenvalue ${\lambda }_n$ of the dynamics generator ${\mathcal{L}}_{k_1,k_2}$ is called a deadbeat mode if ${\lambda }_n=i{\beta }_n,{\beta }_n>0$. If the corresponding eigenfunction is ${\mathrm{\Phi }}_n\left(x\right)$, then the evolution problem (13) has a solution given in the form ${\mathrm{e}}^{i{\lambda }_{n}t}{\mathrm{\Phi }}_n\left(x\right)={\mathrm{e}}^{-{\beta }_{n}t}{\mathrm{\Phi }}_n\left(x\right)$, which tends to zero without any oscillation.

As shown in paper [2], for the case when one of the control parameters is zero and the other one is positive, the entire set of the eigenvalues is located in the closed upper half plane. This result is not obvious since the operator is not dissipative; in fact, it requires a fairly nontrivial proof. However, due to this fact, we assume that any deadbeat mode can be given in the form $\mathit{i\beta }$, with $\beta >0$. To deal with the deadbeat modes analytically, we rewrite the spectral equation $\left({\mathcal{L}}_{k_1,k_2}\mathrm{\Phi }\right)\left(x\right)=\lambda \mathrm{\Phi }\left(x\right)$ in the form of an equivalent problem for an operator pencil [17] as

$$\begin{array}{l}\mathit{EI}\phi \prime \prime \prime \prime \left(x\right)={\lambda }^2\rho \phi \left(x\right),\phi \left(0\right)={\phi }^{\prime }\left(0\right)=0,\\ \mathit{EI}\phi \prime \prime \left(1\right)=-\mathit{i\lambda }{k}_1\phi \left(1\right),\mathit{EI}\phi \prime \prime \prime \left(1\right)=\mathit{i\lambda }{k}_2{\phi }^{\prime }\left(1\right).\end{array}$$E24

If ${\lambda }_n$ and ${\phi }_n\left(x\right)$ are an eigenvalue and eigenfunction of the pencil (24), then ${\lambda }_n$ is also an eigenvalue of ${\mathcal{L}}_{k_1,k_2}$ with the eigenfunction ${\mathrm{\Phi }}_n\left(x\right)={\left[\frac{1}{i{\lambda }_n}{\phi }_n\left(x\right),{\phi }_n\left(x\right)\right]}^T$.

To solve problem (24), we first redefine the spectral and control parameters to eliminate $\rho$ and $\mathit{EI}$ from Eq. (24). We define $\tilde{\lambda },{\tilde{k}}_1$, and ${\tilde{k}}_2$ by $\lambda =\sqrt{\mathit{EI}/\rho }\tilde{\lambda }$ and ${\tilde{k}}_j=\sqrt{\mathit{EI}\rho }{k}_j$, $j=1,2$. Substituting these relations into Eq. (24) and eliminating the “tilde,” we obtain the following Sturm-Liouville eigenvalue problem:

$$\phi \prime \prime \prime \prime \left(x\right)={\lambda }^2\phi \left(x\right),\phi \left(0\right)={\phi }^{\prime }\left(0\right)=0,\phi \prime \prime \left(1\right)=-\mathit{i\lambda }{k}_1\phi \left(1\right),\phi \prime \prime \prime \left(1\right)=\mathit{i\lambda }{k}_2{\phi }^{\prime }\left(1\right).$$E25

The solution of Eq. (25) satisfying the left-end boundary conditions $\phi \left(0\right)={\phi }^{\prime }\left(0\right)=0$ can be written in the form

$$\phi \left(\lambda ,x\right)=\mathcal{A}\left(\lambda \right)\left[cosh\left(\sqrt{\lambda }x\right)-cos\left(\sqrt{\lambda }x\right)\right]+\mathcal{B}\left(\lambda \right)\left[sinh\left(\sqrt{\lambda }x\right)-sin\left(\sqrt{\lambda }x\right)\right].$$E26

Substituting formula (26) into the right-end boundary conditions of Eq. (25), one gets a system for the coefficients $\mathcal{A}\left(\lambda \right)$ and $\mathcal{B}\left(\lambda \right)$:

$$\begin{array}{l}\mathcal{A}\left(\lambda \right)\left[\left(1+{\mathit{ik}}_1\right)cosh\sqrt{\lambda }-\left(1-{\mathit{ik}}_1\right)cos\sqrt{\lambda }\right]+\\ \mathcal{B}\left(\lambda \right)\left[\left(1+{\mathit{ik}}_1\right)sinh\sqrt{\lambda }-\left(1-{\mathit{ik}}_1\right)sin\sqrt{\lambda }\right]=0,\\ \mathcal{A}\left(\lambda \right)\left[\left(1-{\mathit{ik}}_2\right)sinh\sqrt{\lambda }-\left(1+{\mathit{ik}}_2\right)sin\sqrt{\lambda }\right]+\\ \mathcal{B}\left(\lambda \right)\left[\left(1-{\mathit{ik}}_2\right)cosh\sqrt{\lambda }+\left(1+{\mathit{ik}}_2\right)cos\sqrt{\lambda }\right]=0.\end{array}$$E27

Let $\mathrm{\Delta }\left(\lambda \right)$ be the determinant of the matrix of coefficients for $\mathcal{A}\left(\lambda \right)$ and $\mathcal{B}\left(\lambda \right)$ in Eqs. (27). System (27) has nontrivial solutions if and only if $\mathrm{\Delta }\left(\lambda \right)=0$, i.e.,

$$\left(1+k_1{k}_2\right)+\left(1-k_1{k}_2\right)cosh\sqrt{\lambda }cos\sqrt{\lambda }+i\left(k_1+k_2\right)sinh\sqrt{\lambda }sin\sqrt{\lambda }=0.$$E28

Theorem 1: The following results hold in the case when $k_1>0$ and $k_2=0$. Similar results hold in the case when $k_1=0$ and $k_2>0$.

1. For $0<k_1<1$, the deadbeat modes do not exist.

2. For $k_1=1$, there exist infinitely many deadbeat modes given explicitly by

   $${\lambda }_n={\mu }_n^2,{\mu }_n=x_n\left(1+i\right),x_n=\frac{\pi }{2}\left(2n+1\right),n=\mathrm{0,1,2},\dots .$$E29

3. For any $k_1>1$, there exist a finite number $\mathcal{N}\left(k_1\right)$ of deadbeat modes. Each mode has the form $\lambda ={\mu }^2,\mu =x\left(1+i\right)$, where $x$ is a root of the function

$$H\left(x;k_1\right)\equiv 2+\left(1+k_1\right)cos2x+\left(1-k_1\right)cosh2x,x>0.$$E30

Let $X\left(k_1\right)=\frac{1}{2\pi }{cosh}^{-1}\left(1+\frac{4}{k_1-1}\right)$, and then $\mathcal{N}\left(k_1\right)$ satisfies the estimate

$$2\left[X\left(k_1\right)\right]+1\le \mathcal{N}\left(k_1\right)\le 2\left[X\left(k_1\right)\right]+3.$$E31

Hence $\mathcal{N}\left(k_1\right)\to \infty$ as $k_1\to 1^{+}$. (By $\left[X\right]$ we denote the greatest integer less than or equal to $X$).

Proof: Let $\mu =\sqrt{\lambda }=x\left(1+i\right),x>0$. Taking into account the relations

$$\begin{array}{l}2cosh\mu cos\mu =cosh\left(1+i\right)\mu +cosh\left(1-i\right)\mu ,\\ 2isinh\mu sin\mu =cosh\left(1+i\right)\mu -cosh\left(1-i\right)\mu ,\end{array}$$

we reduce Eq. (28) to the following form:

$$2+\left(1+k_1\right)cos2x+\left(1-k_1\right)cosh2x=0,x>0.$$E32

It can be readily seen that if $0<k_1<1$, then $2+\left(1+k_1\right)cos2x>0$ and $\left(1-k_1\right)cosh2x>0$, which means that Eq. (32) has no solutions. Statement (1) is shown. Statement (2) follows immediately if one considers Eq. (32) for $k_1=1$.

To prove Statement (3), we rewrite Eq. (32) in the form

$$cosh2x=\frac{2}{k_1-1}+\left(1+\frac{2}{k_1-1}\right)cos2x,x>0,k_1>1.$$E33

The left-hand side of Eq. (33) is monotonically increasing, while the right-hand side is sinusoidal, with maximum $\left(1+\frac{4}{k_1-1}\right)$ and minimum $\left(-1\right)$, and period $\pi$. So the graphs of the left- and right-hand side have intersections only on the interval $\left[0,\frac{1}{2}{cosh}^{-1}\left(1+\frac{4}{k_1-1}\right)\right]$. There are two intersections for each full period of the right-hand side that fits into the above interval (Figure 4). As it can be seen in Figure 4, one should add at least one more intersection for the first half-period after the full periods. Depending on the value of $k_1$, the two graphs can have two intersections, one tangential intersection or no intersections on the second half-period. This leads to estimate (31). $■$

Figure 4.

Left- and right-hand side of Eq. (33) for different values of $k_1$.

A graphical illustration of the result of Theorem 1 is shown in Figure 5.

Figure 5.

Estimates and actual count of deadbeat modes based on numerical simulations.

4. Structure of the deadbeat mode set

The main result on the existence and distribution of double roots of the function $H\left(x;k_1\right)$ is presented in the statement below.

Theorem 2: For a given $k_1>1$, the multiplicity of each root of $H\left(x;k_1\right)$ does not exceed $2$. There exists a sequence $\left\{k_1^{\left(n\right)};n=\mathrm{0,1,2},\dots \right\}$, such that the function $H\left(x;k_1\right)$ has a double root if and only if $k_1=k_1^{\left(n\right)}$ for some $n$. So the original spectral problem with $k_1=k_1^{\left(n\right)},k_2=0$ has a double deadbeat mode ${\lambda }_n={\mu }_n^2=2{\mathit{ix}}_n^2$. The following asymptotic formulas hold

$$x_n=\frac{3\pi }{4}+\mathit{\pi n}+P_n^{-1}+O\left(P_n^{-2}\right),\text{where}{P}_n=exp\left\{\frac{3\pi }{2}+2\mathit{\pi n}\right\},$$E34

and

$$k_1^{\left(n\right)}=1+4P_n^{-1}+O\left(P_n^{-2}\right).$$E35

Proof: If $x$ is a double root of $H$, then $H\left(x;k_1\right)=H^{\prime }\left(x;k_1\right)=0$, i.e., separating the real and imaginary parts, we have

$$2+\left(k_1+1\right)cos2x-\left(k_1-1\right)cosh2x=0,$$E36

$$\left(k_1+1\right)sin2x+\left(k_1-1\right)sinh2x=0.$$E37

Eliminating $k_1$ from system given by (36) and (37), we obtain that the following equation has to be satisfied:

$$G\left(x\right)\equiv \left(1+cos2x\right)sinh2x+\left(1+cosh2x\right)sin2x=0,x>0.$$E38

Rewriting Eq. (37) in the form $\left(k_1+1\right)sin2x=-\left(k_1-1\right)sinh2x$, and taking into account that $k_1>1$, we obtain that if $x$ is the solution of Eq. (38), then $sin2x<0$.

Now we show that when $cos2x<0$ and $sin2x<0$, i.e.,

$$\pi \left(2n+1\right)<2x<\frac{3\pi }{2}+2\mathit{\pi n},n\in \left\{\mathrm{0,1,2},\dots \right\},$$

Eq. (38) does not have any solutions. Indeed, in the above range of $x$, we have $cos2x+sin2x=\sqrt{2}sin\left(2x+\pi /4\right)<-1$ and $\mid sin2x-cos2x\mid <1$. With such estimates we obtain that

$$G\left(x\right)=sin2x+sinh2x+\frac{1}{2}{\mathrm{e}}^{2x}\left(cos2x+sin2x\right)+\frac{1}{2}{\mathrm{e}}^{-2x}\left(sin2x-cos2x\right)<sin2x<0,$$E39

which mean that Eq. (38) cannot be satisfied.

Now we consider the case when $cos2x>0$ and $sin2x<0$, i.e.,

$$\frac{3\pi }{2}+2\mathit{\pi n}<2x<2\pi \left(n+1\right),n\in \left\{\mathrm{0,1,2},\dots \right\}.$$

It is convenient to rewrite system given by (36) and (37) in the form $2x=3\pi /2+2\mathit{\pi n}+s$, where $n\in \left\{\mathrm{0,1,2},\dots \right\}$, $0<s<\pi /2$. If $g\left(s\right)\equiv G\left(3\pi /4+\mathit{\pi n}+s\right)$, then Eq. (38) generates the following equation for $g$:

$$g\left(s\right)\equiv \left(1+sins\right)sinh\left(\frac{3\pi }{2}+2\mathit{\pi n}+s\right)-\left[1+cosh\left(\frac{3\pi }{2}+2\mathit{\pi n}+s\right)\right]coss=0.$$E40

Let us show that for each $n$, Eq. (40) has a unique solution. For $s=0$ we have

$$g\left(0\right)=sinh\left(\frac{3\pi }{2}+2\mathit{\pi n}\right)-1-cosh\left(\frac{3\pi }{2}+2\mathit{\pi n}\right)<0,$$

and for $s=\pi /2$ we have $g\left(\pi /2\right)=2sinh\left(2\pi \left(n+1\right)\right)>0$. Evaluating $g^{\prime }$ we have

$$g^{\prime }\left(s\right)=sins+\left(1+2sins\right)cosh\left(\frac{3\pi }{2}+2\mathit{\pi n}+s\right)>0.$$E41

Thus $g\left(s\right)$ is a monotonically increasing function, such that $g\left(0\right)<0<g\left(\pi /2\right)$, which means that $g$ has a unique root on $\left[0,\pi /2\right]$.

Finally we show that the multiplicity of a multiple root cannot exceed 2. Using a contradiction argument, assume that there exists $x_0$, such that in addition to Eqs. (36) and (37), one has $H\prime \prime \left(x_0;k_1\right)=0$, i.e., the multiplicity of $x_0$ is at least 3. The system $H^{\prime }\left(x_0;k_1\right)=0$ and $H\prime \prime \left(x_0;k_1\right)=0$ can be written as

$$\begin{array}{ll}\left(k_1+1\right)sin2x_0+\left(k_1-1\right)sinh2x_0& =0,\\ \left(k_1+1\right)cos2x_0+\left(k_1-1\right)cosh2x_0& =0.\end{array}$$E42

Since $k_1>1$, the second equation of (42) yields $cos2x_0<0$. Also, since $x_0$ is a multiple root, we must have $sin2x_0<0$. Then $2x_0$ is in the third quadrant, which means that $G\left(x_0\right)\ne 0$, as we have seen above. This contradicts our assumption that $x_0$ is a root of Eqs. (36) and (37).

To derive asymptotic distribution of the roots of Eq. (40), we check that with $P_n$ from Eq. (34), the following approximations are valid:

$$\begin{array}{l}sin\left(2P_n^{-1}\right)=2P_n^{-1}+O\left(P_n^{-3}\right),cos\left(2P_n^{-1}\right)=1-2P_n^{-2}+O\left(P_n^{-4}\right),\\ 2sinh\left(\frac{3\pi }{2}+2\mathit{\pi n}+2P_n^{-1}\right)=P_n+2+P_n^{-1}+O\left(P_n^{-2}\right),\\ 2cosh\left(\frac{3\pi }{2}+2\mathit{\pi n}+2P_n^{-1}\right)=P_n+2+3P_n^{-1}+O\left(P_n^{-2}\right).\end{array}$$E43

Evaluating $g\left(s\right)$ from Eq. (40) for $s=2P_n^{-1}$ and using Eq. (43), we get

$$\begin{array}{l}g\left(2P_n^{-1}\right)=\left[1+sin\left(2P_n^{-1}\right)\right]sinh\left(\frac{3\pi }{2}+2\mathit{\pi n}+2P_n^{-1}\right)-\\ \left[1+cosh\left(\frac{3\pi }{2}+2\mathit{\pi n}+2P_n^{-1}\right)\right]cos\left(2P_n^{-1}\right)=2P_n^{-1}+O\left(P_n^{-2}\right).\end{array}$$E44

Representation (44) implies that there exists $n_0$, such that for all $n\ge n_0$, we have $g\left(2P_n^{-1}\right)>0$. Taking into account that $g\left(0\right)<0$, we obtain that the root $s_n$, $n\ge n_0$, of the function $g\left(s\right)$ is located on the interval $\left(0,2P_n^{-1}\right)$. To find the location of this root more precisely [18], we use linear interpolation. Namely, substituting Eq. (43) into the expression for $g^{\prime }\left(s\right)$ from Eq. (41) yields

$$g^{\prime }\left(2P_n^{-1}\right)=\frac{P_n}{2}+O\left(1\right).$$E45

Replacing $g\left(s\right)$ by the linear function tangential to $g\left(s\right)$ at the point $\left(2P_n^{-1},g\left(2P_n^{-1}\right)\right)$, and finding the root of this function, we get

$$s_n=2P_n^{-1}-\frac{g\left(2P_n^{-1}\right)}{g^{\prime }\left(2P_n^{-1}\right)}+O\left(P_n^{-2}\right)=2P_n^{-1}+O\left(P_n^{-2}\right).$$E46

Having this approximation for $s_n$, we immediately get

$$x_n=\frac{3\pi }{4}+\mathit{\pi n}+\frac{s_n}{2}=\frac{3\pi }{4}+\mathit{\pi n}+P_n^{-1}+O\left(P_n^{-2}\right).$$E47

From the equation $H\left(x_n;k_1^{\left(n\right)}\right)=0$, we obtain the formula for $k_1^{\left(n\right)}$ as

$$k_1^{\left(n\right)}=\frac{sinh\left(3\pi /2+2\mathit{\pi n}+s_n\right)+cos{s}_n}{sinh\left(3\pi /2+2\mathit{\pi n}+s_n\right)-cos{s}_n}.$$E48

Substituting formulas (43) and (46) into formula (48), we obtain representation (35).$■$

Corollary 1: Let $k_1=k_1^{\left(n\right)}$ for some $n\in {\mathrm{\mathbb{N}}}^{+}\cup \left\{0\right\}$, and let $x_n$ be the corresponding double root of the function $H\left(x;k_1^{\left(n\right)}\right)$. Then ${\lambda }_0=2{\mathit{ix}}_n^2$ is an eigenvalue of the operator ${\mathcal{L}}_{k_1^{\left(n\right)},0^{\prime }}$, such that the geometric multiplicity of ${\lambda }_0$ is 1 and its algebraic multiplicity is 2. Therefore there exists a unique eigenvector $\mathrm{\Phi }$ and one associate vector $\mathrm{\Psi }$, such that

$${\mathcal{L}}_{k_1^{\left(n\right)},0}\mathrm{\Phi }={\lambda }_0\mathrm{\Phi },{\mathcal{L}}_{k_1^{\left(n\right)},0}\mathrm{\Psi }-{\lambda }_0\mathrm{\Psi }=\mathrm{\Phi }.$$E49

Proof: It suffices to show that problem (24) does not have two linearly independent eigenvectors corresponding to ${\lambda }_0$ [18, 19]. Using contradiction argument we assume that for some ${\lambda }_0$ the boundary-value problem (25) with $k_2=0$ has two linearly independent solutions $\psi$ and $\chi$. Each function satisfies the problem

$$\phi \prime \prime \prime \prime \left(x\right)={\lambda }_0^2\phi \left(x\right),\phi \left(0\right)={\phi }^{\prime }\left(0\right)=0,\phi \prime \prime \prime \left(1\right)=0,\phi \prime \prime \left(1\right)=i{\lambda }_0{k}_1\phi \left(1\right).$$E50

First we observe that $\psi \left(1\right)\chi \left(1\right)\ne 0$. Indeed, if $\psi \left(1\right)=0$, then we have

$${\int }_0^1{\left|\psi \prime \prime \left(\sigma \right)\right|}^2\mathrm{d\sigma }={\int }_0^1\psi \prime \prime \prime \prime \left(\sigma \right)\overline{\psi \left(\sigma \right)}\mathrm{d\sigma }={\lambda }_0^2{\int }_0^1{\left|\psi \left(\sigma \right)\right|}^2\mathrm{d\sigma }.$$E51

Since ${\lambda }_0$ is purely imaginary, Eq. (51) is not valid. We define a new function:

$$g\left(x\right)=\psi \left(x\right)-\frac{\psi \left(1\right)}{\chi \left(1\right)}\chi \left(x\right).$$E52

One can readily check that $g$ satisfies the following boundary-value problem:

$$\phi \prime \prime \prime \prime \left(x\right)={\lambda }_0^2\phi \left(x\right),\phi \left(0\right)={\phi }^{\prime }\left(0\right)=0,\phi \left(1\right)=\phi \prime \prime \left(1\right)=0,$$E53

and therefore

$${\int }_0^1{\left|{g^{\prime }}^{\prime }\left(\sigma \right)\right|}^2\mathrm{d}\sigma ={\int }_0^{1}g\prime \prime \prime \prime \left(\sigma \right)\overline{g\left(\sigma \right)}\sigma ={\lambda }_0^2{\int }_0^1{\left|g\left(\sigma \right)\right|}^2\mathrm{d}\sigma .$$E54

Eq. (54) is valid if and only if ${\lambda }_0^2>0$; however, for a deadbeat mode, ${\lambda }_0^2<0$. The obtained contradiction means that for each double mode, there is one eigenfunction and one associate function.$■$

4.1 Deadbeat mode behavior as ${\mathit{k}}_{\mathbf{1}}\to \mathbf{1}$

As $k_1$ approaches $1$, the spectral branches are moving upward and toward the imaginary axis (Figures 6 and 7). As a result of this motion, eigenvalues approach the imaginary axis at different rates depending on whether $k_1$ approaches $1$ from above or below.

Figure 6.

Eigenvalues with $\mid Re\lambda \mid >10$ as $k_1\to 1^{+}$.

Figure 7.

Eigenvalues with $\mid Re\lambda \mid >10$ as $k_1\to 1^{-}$.

As follows from Table 3, the real parts of the eigenvalues decrease steadily as $k_1\to 1^{+}$, to a point where the eigenvalue becomes a deadbeat mode. An increase in the number of deadbeat modes can be seen as $k_1\to 1^{+}$, which is in agreement with Statement (3) of Theorem 1. One can see from Table 3 that there are pairs of modes such that the distance between them tends to zero as $k_1\to 1^{+}$. (Compare modes no.5 and no.7 for $\mid k_1-1\mid ={10}^{-4}$, modes no.4 and no.7 for $\mid k_1-1\mid ={10}^{-6}$, modes no.5 and no.8 for $\mid k_1-1\mid ={10}^{-8}$, and modes no.4 and no.7 for $\mid k_1-1\mid ={10}^{-10}$). Such behavior indicates convergence of the two simple deadbeat modes to a double mode, which is consistent with Theorem 2.

Analyzing Table 4, one can see that the eigenvalues get closer to the imaginary axis as $k_1\to 1^{-}$. However the rate at which their real parts approach zero is significantly lower than in the case $k_1\to 1^{+}$. Even at $k_1=1-{10}^{-10}$, the eigenvalue closest to the imaginary axis has a real part of about ${10}^{-4}$, which means that it is not a deadbeat mode (see Statement (1) of Theorem 1).

The eigenvalues near the imaginary axis approach the same double deadbeat modes in both cases when $k_1\to 1^{\pm }$ (see Statement (2) of Theorem 1). In conclusion, one can claim that the eigenvalues are indeed approaching the imaginary axis; however, the rate of this approach is different for $k_1\to 1^{-}$ and $k_1\to 1^{+}$. In the former case, an eigenvalue’s distance from the imaginary axis decreases very slowly; in the latter case, the eigenvalues quickly “jump” on the imaginary axis and turn into deadbeat modes.

5. Outline of the numerical scheme

To carry out the numerical analysis of the differential operator ${\mathcal{L}}_{k_1,k_2}$, we use the Chebyshev collocation method and cardinal functions [14, 15, 16].

Recall that the $N$th Chebyshev polynomial of the first kind is defined by

$$T_N\left(\xi \right)=cos\mathit{N\theta },-1\le \xi \le 1\text{where}\xi =cos\theta ,\theta \in \left[0,\pi \right].$$E55

The cardinal functions, ${\mathrm{\psi }}_k\left(\xi \right)$, and the Chebyshev-Gauss-Lobatto (CGL) grid points $\left\{{\xi }_k\right\}$ are defined as follows:

$${\psi }_k\left(\xi \right)={\left(-1\right)}^k\frac{\left(1-{\xi }^2\right)T_{N-1}^{\prime }\left(\xi \right)}{c_k{\left(N-1\right)}^2\left(\xi -{\xi }_k\right)},{\xi }_k=cos\frac{\left(k-1\right)\pi }{N-1},\text{for}1\le k\le N,$$E56

where coefficients $c_k$ are such that $c_1=c_N=2$ and $c_k=1$ for $1<k<N$. The main property of cardinal functions is ${\psi }_k\left({\xi }_j\right)={\delta }_{\mathit{kj}}$ (using the Kronecker delta). The family ${\left\{{\psi }_k\right\}}_{k=1}^N$ forms a basis in the space of polynomials of degree $\left(N-1\right)$, i.e., if $f$ is such polynomial, then $f$ and $f^{\prime }$ can be written in the forms

$$f\left(\xi \right)=\sum _{k=1}^{N}f\left({\xi }_k\right){\psi }_k\left(\xi \right)\text{and}{f}^{\prime }\left(\xi \right)=\sum _{k=1}^N{f}^{\prime }\left({\xi }_k\right){\psi }_k\left(\xi \right).$$E57

If $\mathit{f}={\left[f\left({\xi }_1\right),f\left({\xi }_2\right),\dots ,f\left({\xi }_N\right)\right]}^T$ and $\mathit{g}={\left[f^{\prime }\left({\xi }_1\right),f^{\prime }\left({\xi }_2\right),\dots ,f^{\prime }\left({\xi }_N\right)\right]}^T$, then $\mathit{g}=D\mathit{f}$, where $D$ is the Chebyshev derivative matrix with the elements

$$\begin{array}{c}{D}_{11}=-D_{\mathit{NN}}=\frac{1+2{\left(N-1\right)}^2}{6},D_{\mathit{kk}}=-\frac{{\xi }_k}{2\left(1-{\xi }_k^2\right)}\text{for}1<k<N,\\ D_{j,k}=\frac{c_j{\left(-1\right)}^{j+k}}{c_k\left({\xi }_j-{\xi }_k\right)}\text{for}j\ne k.\end{array}$$E58

5.1 Discretization of ${\mathcal{L}}_{{\mathit{k}}_{\mathbf{1}},{\mathit{k}}_{\mathbf{2}}}$

Rescaling the independent variable $x$ as $\xi =2x-1$, we rewrite the operator and its domain, representations (14) and (15), in the form

$${\mathcal{L}}_{k_1,k_2}=-i\left[\begin{array}{cc}0& 1\\ -\frac{16}{\rho \left(\xi \right)}\frac{{\partial }^2}{\partial {\xi }^2}\left(\mathit{EI}\left(\xi \right)\frac{{\partial }^2}{\partial {\xi }^2}\right)& 0\end{array}\right],$$E59

and

$$\begin{array}{l}\mathcal{D}\left({\mathcal{L}}_{k_1,k_2}\right)=\{{\left(u_0,u_1\right)}^T\in \mathcal{H}:u_0\in H^4\left(-1,1\right),u_1\in H_0^2\left(-1,1\right)\mathrm{;}{u}_1\left(-1\right)=u_1^{\prime }\left(-1\right)=0\mathrm{;}\\ {\left.4\mathit{EI}\left(1\right)u_0^{″}\left(1\right)=-k_1{u}_1\left(1\right),4{\left(\mathit{EI}\left(\xi \right)u_0^{″}\left(\xi \right)\right)}^{\prime }\right|}_{\xi =1}=k_2{u}_1^{\prime }\left(1\right)\},\end{array}$$E60

where $\mathcal{H}=H_0^2\left(-1,1\right)\times L^2\left(-1,1\right)$, equipped with the norm

$${‖U‖}_{\mathcal{H}}^2=\frac{1}{4}{\int }_{-1}^1\left[16\mathit{EI}\left(\xi \right){\left|u_0^{″}\left(\xi \right)\right|}^2+\rho \left(\xi \right){\left|u_1\left(\xi \right)\right|}^2\right]\mathrm{d\xi }.$$E61

We approximate the action of ${\mathcal{L}}_{k_1,k_2}$ on the finite-dimensional subspace ${\mathcal{H}}_N\subset \mathcal{H}$ of polynomials of degree at most $\left(N-1\right)$. Using the CGL grid and the cardinal functions, we substitute for $u_0$ and $u_1$ their truncated expansions:

$$u_0\left(\xi \right)\approx \sum _{k=1}^N{\mathrm{\Phi }}_k{\psi }_k\left(\xi \right),{\mathrm{\Phi }}_k=u_0\left({\xi }_k\right),u_1\left(\xi \right)\approx \sum _{k=1}^N{\mathrm{\Theta }}_k{\psi }_k\left(\xi \right),{\mathrm{\Theta }}_k=u_1\left({\xi }_k\right).$$E62

Let $\mathbf{\Phi }$ and $\mathbf{\Theta }$ be $N$-dim vectors and $\mathbf{\Psi }$ be a $2N$-dim vector defined by

$$\mathbf{\Phi }={\left[{\mathrm{\Phi }}_1,{\mathrm{\Phi }}_2,\dots ,{\mathrm{\Phi }}_N\right]}^T,\mathbf{\Theta }={\left[{\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\dots ,{\mathrm{\Theta }}_N\right]}^T,\mathbf{\Psi }=\left[\begin{array}{l}\mathbf{\Phi }\\ \mathbf{\Theta }\end{array}\right].$$E63

Let $\mathit{L}$ be the finite-dimensional approximation of the differential operator ${\mathcal{L}}_{k_1,k_2}$. The discretized operator $\mathit{L}$ induced by ${\mathcal{L}}_{k_1,k_2}$ can be given by

$$\mathit{L}=-i\left[\begin{array}{cc}\mathbf{0}& I_{N\times N}\\ -16\frac{\mathit{EI}}{\rho }{D}^4& \mathbf{0}\end{array}\right],$$E64

where $I_{N\times N}$ is the $N\times N$ identity matrix and $D$ is the derivative matrix (58).

5.2 Incorporating the boundary conditions

Discretization of the boundary conditions in the domain description (60) yields

$$\begin{array}{c}{\mathrm{\Phi }}_N=0,{\left[D\mathbf{\Phi }\right]}_N=0,{\mathrm{\Theta }}_N=0,{\left[D\mathbf{\Theta }\right]}_N=0,\\ 4\mathit{EI}{\left[D^2\mathbf{\Phi }\right]}_1+k_1{\mathrm{\Theta }}_1=0,4\mathit{EI}{\left[D^3\mathbf{\Phi }\right]}_1-k_2{\left[D\mathbf{\Theta }\right]}_1=0.\end{array}$$E65

Let $r_N,z_N,l_N\in {\mathbb{R}}^N$ be auxiliary row-vectors

$$r_N=\left[00\cdots 01\right],z_N=\left[00\cdots 00\right],l_N=\left[10\cdots 00\right]$$E66

and $D_j^n$ designate the $j$th row of the $n$th derivative matrix $D^n$. Using Eqs. (66) we represent Eqs. (65) as the following matrix equation:

$$\mathbb{K}\mathbf{\Psi }\equiv \left[\begin{array}{cc}{r}_N& z_N\\ D_N^1& z_N\\ z_N& r_N\\ z_N& D_N^1\\ 4{\mathit{EID}}_1^2& k_1{l}_N\\ 4{\mathit{EID}}_1^3& -k_2{D}_1^1\end{array}\right]\left[\begin{array}{l}\mathbf{\Phi }\\ \mathbf{\Theta }\end{array}\right]=0.$$E67

$\mathbb{K}$ is called the boundary operator. Let ${\mathcal{K}}_N$ be the kernel of $\mathbb{K}$, i.e., ${\mathcal{K}}_N=\left\{v\in {\mathbb{R}}^{2N}:\mathbb{K}v=0\right\}$. We have to identify all eigenvalues of the operator $\mathbf{L}$, when its domain is restricted to ${\mathcal{K}}_N$. It is clear that ${\mathcal{K}}_N$ is isomorphic to ${\mathbb{R}}^k$ with $k\equiv dim{\mathcal{K}}_N=dim{\mathcal{H}}_N-rank\mathbb{K}=2N-6$. Let $B$ be the matrix consisting of column vectors that form an orthonormal basis in ${\mathcal{K}}_N$. It is clear that $B^{T}B$ is the identity matrix on ${\mathbb{R}}^k$ and ${\mathit{BB}}^T$ is the identity matrix on $\mathbb{K}$. The following result holds: if $\lambda$ is an eigenvalue of the operator $\mathit{L}$, and the corresponding eigenvector $\mathbf{\Psi }$ satisfies Eq. (67), then the same $\lambda$ is an eigenvalue of the matrix $\left(B^T\mathit{L}B\right)$. However, the inverse statement is not necessarily true. Indeed, we observe that ${\mathit{BB}}^T$ is the identity in ${\mathcal{K}}_{\mathcal{N}}$, which is not equivalent to the identity in ${\mathcal{H}}_N$. Assume now that $\lambda$ is an eigenvalue of $B^T\mathit{L}B$ with corresponding eigenvector $v\in {\mathbb{R}}^k$. If $\mathbf{\Psi }=\mathit{Bv}$, we have

$${\mathit{BB}}^T\mathit{L}\mathbf{\Psi }={\mathit{BB}}^T\mathit{L}\mathit{Bv}=\mathit{\lambda Bv}=\lambda \mathbf{\Psi },$$E68

but ${\mathit{BB}}^T\mathit{L}\ne \mathit{L},$ which indicates that fake eigenvalues may exist.

5.3 Filtering of spurious eigenvalues

In order to decide which eigenvalues of $B^T\mathit{L}B$ should be discarded, we impose the following condition. Let $\mathrm{\Lambda }$ be the spectrum of $B^T\mathit{L}B$ and $V$ be the set of its eigenfunctions. We construct the set of “trusted” eigenvalues [14, 15], for some ${\epsilon }_f>0$ filtering precision, as

$${\mathrm{\Lambda }}_{\epsilon }=\left\{\lambda \in \mathrm{\Lambda }:{‖\mathit{L}{\mathit{Bv}}_{\lambda }-\lambda {\mathit{Bv}}_{\lambda }‖}_C<{\epsilon }_f,\text{for\hspace{0.17em}corresponding\hspace{0.17em}eigenvector}{v}_{\lambda }\in V\right\},$$E69

where ${‖\cdot ‖}_C$ is a discrete approximation to the integral norm defined in Eq. (61). (The subscript $C$ is short for Chebyshev). Using the CGL quadrature, we obtain the following formula for the norm of a vector $\mathbf{\Psi }$ defined as in Eq. (63):

$${‖\mathbf{\Psi }‖}_C=\frac{\pi /4}{N-1}\sum _{k=1}^N\sqrt{1-{\xi }_k^2}\left[16\mathit{EI}\left({\xi }_k\right){\left|{\left[D^2\mathbf{\Phi }\right]}_k\right|}^2+\rho \left({\xi }_k\right){\left|{\mathrm{\Theta }}_k\right|}^2\right].$$

6. Conclusions

In this work we have considered the spectral properties of the Euler-Bernoulli beam model with special feedback-type boundary conditions. The dynamics generator of the model is a non-self-adjoint matrix differential operator acting in a Hilbert space of two-component Cauchy data. This operator has been approximated by a “discrete” operator using Chebyshev polynomial approximation. We have shown that the eigenvalues of the main operator can be approximated by the eigenvalues of its discrete counterpart with high accuracy. This means that the leading asymptotic terms in formulas (20) and (21) can be used by practitioners who need the elastic modes.

Further results deal with existence and formulas of the deadbeat modes. It has been shown that for the case when one control parameter, $k_1$, is such that $k_1\to 1^{+}$ and the other one $k_2=0$, the number of deadbeat modes approaches infinity. The formula for the rate at which the number of the deadbeat modes tends to infinity has been derived. It has also been established that there exists a sequence ${\left\{k_1^{\left(n\right)}\right\}}_{n=1}^{\infty }$ of the values of parameter $k_1$, such that the corresponding deadbeat mode has a multiplicity 2, which yields the existence of the associate mode shapes for the operator ${\mathcal{L}}_{k_1,k_2}$. The formulas for the double deadbeat modes and asymptotics for the sequence $\left\{k_1^{\left(n\right)}\right\}$ as $n\to \infty$ have been derived.

Acknowledgments

Partial support of the National Science Foundation award DMS-1810826 is highly appreciated by the first author.