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

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 ( k 1 and k 2 ) 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, k 1 , is positive and another one, k 2 , is zero. It has been shown that the number of the deadbeat modes tends to infinity, as k 1 ! 1 þ and k 2 ¼ 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.


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 nondissipative 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 initialboundary value problem describing the beam dynamics is reduced to the first order in time evolution equation in the state Hilbert space H. The evolution of the system is completely determined by the dynamics generator 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 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 ≥ 0 and k 2 ≥ 0, such that |k 1 | þ |k 2 | . 0 (see Proposition 2). As shown in [2], even though the operator 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 C, which means that all eigenmodes are stable or neutrally stable. (We recall that to obtain an elastic mode from an eigenvalues of 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 in formulas (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 ! 1 þ .
We have also shown that there exists a sequence of values of the parameter k 1 , i.e., k n ð Þ 1 n o ∞ n¼1 , such that for each k 1 ¼ k n ð Þ 1 there exist a finite number of deadbeat modes and each corresponds to a double eigenvalue of the dynamics generator L k 1 , k 2 .
For each value k n ð Þ 1 , the operator 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 n ð Þ 1 ), 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 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 initialboundary 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 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 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.

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] 1 2 where h x; t ð Þ is the transverse deflection, E x ð Þ is the modulus of elasticity, I x ð Þ is the area moment of inertia, ϱ x ð Þ is the linear density, and A x ð Þ is the cross-sectional area of the beam.
Assuming that the beam is clamped at the left end x ¼ 0 ð Þand free at the right end x ¼ 1 ð Þ, and applying Hamilton's variational principle to the action functional defined by (1), we obtain the equation of motion where M x; t ð Þ and Q x; t ð Þ are the moment and the shear, respectively [10]:   Table 2.
Approximations of the eigenvalues for the discrete and "continuous" operators (K , À 1). 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 where T stands for transposition. The feedback control law is given by where K is the 2 Â 2 feedback matrix. We select Graphical representation of the eigenvalues of the discrete and "continuous" operators (K , À 1). with k 1 , k 2 being the control parameters. The feedback (6) can be written as Finally, we arrive to the following initial-boundary value problem: the equation of motion (2), the boundary conditions (3), and the standard initial conditions Notice that the choice of a feedback matrix K defines whether the system is dissipative or not. Indeed, let E t ð Þ be the energy of the system, defined by representation (1). Evaluating E t t ð Þ on the solutions of Eq. (2) satisfying the left-end conditions from Eqs. (3), we obtain Taking into account Eqs. (4) and (6), we represent the right-hand side of Eq. (9) as the dot product in R 2 : With the choice of K as in Eq. (7), we have Thus the system is not dissipative for all nonnegative values of k 1 and k 2 .

Operator form of the problem
In what follows, we incorporate the cross-sectional area A x ð Þ into the density, Hilbert space of two-component vector functions U T equipped with the following norm: Assuming that EI, ρ ∈ C 2 0; 1 ð Þ are positive functions, we obtain that the closure of smooth functions U , 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: where 0 ≥ x ≥ 1, t≥ 0. The dynamics generator L k 1 , k 2 is given by the following matrix differential expression: defined on the domain For any k 1 ; k 2 ð Þ∈ R 2 , the adjoint operator L * k 1 , k 2 [13] is given by i.e., L * k 1 , k 2 is defined by the same differential expression (14) on the domain described in Eq. (15), where k 1 and k 2 are replaced by Àk 2 ð Þand Àk 1 ð Þ, respectively. It follows from Eq. (16) that L 0, 0 is self-adjoint in H and thus L 0, 0 is the dynamics generator of the clamped-free beam model. For the reader's convenience, we summarize the properties of L k 1 , k 2 from [1,2] needed for the present work.
2. For each k 1 ; k 2 ð Þ∈ R 2 , |k 1 | þ |k 2 | . 0, the operator L k 1 , k 2 is a rank-two perturbation of the self-adjoint operator L 0, 0 in the sense that the operators L À1 k 1 , k 2 and L À1 0, 0 exist and are related by the rule where T k 1 , k 2 is a rank-two operator. A similar decomposition is valid for the adjoint operator, i.e., 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 k 1 , k 2 . 0 and k 1 k 2 6 ¼ EI ρ. Let The following asymptotic approximations for the eigenvalues λ n (as |n| ! ∞) of the operator L k 1 , k 2 hold: 1. If 1 , |K| , ∞, then for |n| ! ∞ one has 2. If 0 , K , 1, then for n ! ∞ one has 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 ε . 0, one can find a positive integer N, such that all eigenvalues λ n with |n| ≥ N þ 1 satisfy the estimate for the case when 1 , |K| , ∞ and for the case when 0 , |K| , 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)  ½ , and ε 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 K defined by Eq. (19). On the sub-domain with dark gray color K such that |K| . 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).

The deadbeat modes
An eigenvalue λ n of the dynamics generator L k 1 , k 2 is called a deadbeat mode if λ n ¼ iβ n , β n . 0. If the corresponding eigenfunction is Φ n x ð Þ, then the evolution problem (13) has a solution given in the form e iλ n t Φ n x ð Þ ¼ e Àβ n t Φ n x ð Þ, 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 iβ, with β . 0. To deal with the deadbeat modes analytically, we rewrite the spectral equation L k 1 , k 2 Φ ð Þx ð Þ ¼ λΦ x ð Þ in the form of an equivalent problem for an operator pencil [17] as If λ n and φ n x ð Þ are an eigenvalue and eigenfunction of the pencil (24), then λ n is also an eigenvalue of L k 1 , k 2 with the eigenfunction Φ n x ð Þ ¼ 1 To solve problem (24), we first redefine the spectral and control parameters to eliminate ρ and EI from Eq. (24). We defineλ,k 1 , andk 2 by λ ¼ ffiffiffiffiffiffiffiffiffi ffi EI=ρ pλ and k j ¼ ffiffiffiffiffiffiffiffi EI ρ p k j , j ¼ 1, 2. Substituting these relations into Eq. (24) and eliminating the "tilde," we obtain the following Sturm-Liouville eigenvalue problem: The solution of Eq. (25) satisfying the left-end boundary conditions φ 0 ð Þ ¼ φ 0 0 ð Þ ¼ 0 can be written in the form Substituting formula (26) into the right-end boundary conditions of Eq. (25), one gets a system for the coefficients A λ ð Þ and B λ ð Þ: 2. For k 1 ¼ 1, there exist infinitely many deadbeat modes given explicitly by 3. For any k 1 . 1, there exist a finite number N k 1 ð Þ of deadbeat modes. Each mode has the form λ ¼ μ , and then N k 1 ð Þ satisfies the estimate Hence N k 1 ð Þ ! ∞ as k 1 ! 1 þ . (By X ½ we denote the greatest integer less than or equal to X).
To prove Statement (3), we rewrite Eq. (32) in the form The left-hand side of Eq. (33) is monotonically increasing, while the right-hand side is sinusoidal, with maximum 1 þ 4 k 1 À1 and minimum À1 ð Þ, and period π. So the graphs of the left-and right-hand side have intersections only on the interval 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 halfperiod. This leads to estimate (31). ■ A graphical illustration of the result of Theorem 1 is shown in Figure 5.

Structure of the deadbeat mode set
The main result on the existence and distribution of double roots of the function H x; k 1 ð Þis presented in the statement below.
To derive asymptotic distribution of the roots of Eq. (40), we check that with P n from Eq. (34), the following approximations are valid: Evaluating g s ð Þ from Eq. (40) for s ¼ 2P À1 n and using Eq. (43), we get Representation (44) implies that there exists n 0 , such that for all n ≥ n 0 , we have g 2P À1 n À Á . 0. Taking into account that g 0 ð Þ , 0, we obtain that the root s n , n ≥ n 0 , of the function g s ð Þ is located on the interval 0; 2P À1 n À Á . To find the location of this root more precisely [18], we use linear interpolation. Namely, substituting Eq. (43) into the expression for g 0 s ð Þ from Eq. (41) yields Replacing g s ð Þ by the linear function tangential to g s ð Þ at the point 2P À1 n ; g 2P À1 n À Á À Á , and finding the root of this function, we get Having this approximation for s n , we immediately get . Then λ 0 ¼ 2ix 2 n is an eigenvalue of the operator L k n ð Þ 1 , 0 0 , such that the geometric multiplicity of λ 0 is 1 and its algebraic multiplicity is 2. Therefore there exists a unique eigenvector Φ and one associate vector Ψ, such that Proof: It suffices to show that problem (24) does not have two linearly independent eigenvectors corresponding to λ 0 [18,19]. Using contradiction argument we assume that for some λ 0 the boundary-value problem (25) with k 2 ¼ 0 has two linearly independent solutions ψ and χ. Each function satisfies the problem First we observe that ψ 1 ð Þχ 1 ð Þ 6 ¼ 0. Indeed, if ψ 1 ð Þ ¼ 0, then we have Since λ 0 is purely imaginary, Eq. (51) is not valid. We define a new function: One can readily check that g satisfies the following boundary-value problem: and therefore Eq. (54) is valid if and only if λ 2 0 . 0; however, for a deadbeat mode, λ 2 0 , 0. The obtained contradiction means that for each double mode, there is one eigenfunction and one associate function. ■

Deadbeat mode behavior as k 1 ! 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.
As follows from Table 3, the real parts of the eigenvalues decrease steadily as k 1 ! 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 ! 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 ! 1 þ . (Compare modes no.5 and no.7 for |k 1 À 1| ¼ 10 À4 , modes no.4 and no.7 for |k 1 À 1| ¼ 10 À6 , modes no.5 and no.8 for |k 1 À 1| ¼ 10 À8 , and modes no.4 and no.7 for |k 1 À 1| ¼ 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 ! 1 À . However the rate at which their real parts approach zero is significantly lower than in the case k 1 ! 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 ! 1 AE (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 ! 1 À and k 1 ! 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.

Outline of the numerical scheme
To carry out the numerical analysis of the differential operator L k 1 , k 2 , we use the Chebyshev collocation method and cardinal functions [14][15][16].
Recall that the Nth Chebyshev polynomial of the first kind is defined by The cardinal functions, ψ k ξ ð Þ, and the Chebyshev-Gauss-Lobatto (CGL) grid points ξ k f g are defined as follows: 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 ψ k ξ j ¼ δ kj (using the Kronecker delta). The family ψ k f g N k¼1 forms a basis in the space of polynomials of degree N À 1 ð Þ, i.e., if f is such polynomial, then f and f 0 can be written in the forms Eigenvalues with |Reλ| . 10 as k 1 ! 1 À .

Incorporating the boundary conditions
Discretization of the boundary conditions in the domain description (60) yields Let r N , z N , l N ∈ R N be auxiliary row-vectors and D n j designate the jth row of the nth derivative matrix D n . Using Eqs. (66) we represent Eqs. (65) as the following matrix equation: K is called the boundary operator. Let K N be the kernel of K, i.e., We have to identify all eigenvalues of the operator L, when its domain is restricted to K N . It is clear that K N is isomorphic to R k with k dimK N ¼ dimH N À rankK ¼ 2N À 6. Let B be the matrix consisting of column vectors that form an orthonormal basis in K N . It is clear that B T B is the identity matrix on R k and BB T is the identity matrix on K. The following result holds: if λ is an eigenvalue of the operator L, and the corresponding eigenvector Ψ satisfies Eq. (67), then the same λ is an eigenvalue of the matrix B T LB À Á . However, the inverse statement is not necessarily true. Indeed, we observe that BB T is the identity in K N , which is not equivalent to the identity in H N . Assume now that λ is an eigenvalue of B T LB with corresponding eigenvector v ∈ R k . If Ψ ¼ Bv, we have but BB T L 6 ¼ L, which indicates that fake eigenvalues may exist.

Filtering of spurious eigenvalues
In order to decide which eigenvalues of B T LB should be discarded, we impose the following condition. Let Λ be the spectrum of B T LB and V be the set of its eigenfunctions. We construct the set of "trusted" eigenvalues [14,15], for some ε f . 0 filtering precision, as where Á k k 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 Ψ defined as in Eq. (63):

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 ! 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 k n ð Þ 1 n o ∞ n¼1 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 L k 1 , k 2 . The formulas for the double deadbeat modes and asymptotics for the sequence k n ð Þ 1 n o as n ! ∞ have been derived.