Open access peer-reviewed chapter

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

Written By

Marianna A. Shubov and Laszlo P. Kindrat

Reviewed: 20 May 2019 Published: 20 June 2019

DOI: 10.5772/intechopen.86940

From the Edited Volume

## Functional Calculus

Edited by Kamal Shah and Baver Okutmuştur

Chapter metrics overview

View Full Metrics

## 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.

### Keywords

• matrix differential operator
• eigenvalues
• Chebyshev polynomials
• numerical scheme
• boundary control

## 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 k1 and k2, 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 H. The evolution of the system is completely determined by the dynamics generator Lk1,k2, 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 Lk1,k2.

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 k1 and k2. Having in mind the practical applications of the asymptotic formulas [3, 4, 5], we discuss the case of k10 and k20, such that k1+k2>0 (see Proposition 2). As shown in [2], even though the operator Lk1,k2 is non-dissipative, for the case k1>0 and k2=0 (or k1=0 and k2>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 Lk1,k2, 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 k1 is close to 1 (while k2=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 k11+.

We have also shown that there exists a sequence of values of the parameter k1, i.e., k1nn=1, such that for each k1=k1n there exist a finite number of deadbeat modes and each corresponds to a double eigenvalue of the dynamics generator Lk1,k2. For each value k1n, the operator Lk1,k2 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 k1n), there exists a mode shape and an associate mode shape. This fact indicates that for some values of k1 and k2, there exists a significant number of associate vectors of Lk1,k2. 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 Lk1,k2, which is a non-self-adjoint matrix differential operator depending on two parameters,k1andk2, is the main object of interest. The eigenvalues and the generalized eigenvectors of Lk1,k2 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.

k1=0, k2=0.5, EI=1, ρ=0.1, N=64, εf=1020
No.NumericalAnalyticNo.NumericalAnalytic
1.7988.1+237.00i7988.1+237.00i18.28.860+13.515i29.453+14.812i
2.7020.6+222.19i7020.6+222.19i19.123.16+29.723i123.08+29.625i
3.6115.5+207.37i6115.5+207.37i20.279.13+44.431i279.14+44.437i
4.5272.8+192.56i5272.8+192.56i21.497.61+59.250i497.61+59.250i
5.4492.5+177.75i4492.5+177.75i22.778.50+74.062i778.50+74.062i
6.3774.7+162.94i3774.7+162.94i23.1121.8+88.875i1121.8+88.875i
7.3119.3+148.12i3119.3+148.12i24.1527.6+103.69i1527.6+103.69i
8.2526.3+133.31i2526.3+133.31i25.1995.7+118.50i1995.7+118.50i
9.1995.7+118.50i1995.7+118.50i26.2526.3+133.31i2526.3+133.31i
10.1527.6+103.69i1527.6+103.69i27.3119.3+148.12i3119.3+148.12i
11.1121.8+88.875i1121.8+88.875i28.3774.7+162.94i3774.7+162.94i
12.778.5+74.062i778.5+74.062i29.4492.5+177.75i4492.5+177.75i
13.497.61+59.250i497.61+59.250i30.5272.8+192.56i5272.8+192.56i
14.279.13+44.431i279.14+44.437i31.6115.5+207.37i6115.5+207.37i
15.123.16+29.723i123.08+29.625i32.7020.6+222.19i7020.6+222.19i
16.28.860+13.515i29.453+14.812i33.7988.1+237.00i7988.1+237.00i
17.2.21017+4.6007i

### Table 1.

Approximations of the eigenvalues for the discrete and “continuous” operators (K>1).

k1=1.3, k2=1.2, EI=10, ρ=0.1, N=64, εf=1020
No.NumericalAnalyticNo.NumericalAnalytic
1.25266229.07i25266229.07i18.99.46719.816i98.17714.317i
2.22206214.76i22206214.76i19.394.1728.149i394.2628.634i
3.19344200.44i19344200.44i20.887.7542.983i887.7542.951i
4.16679186.12i16679186.12i21.1578.657.267i1578.657.268i
5.14212171.81i14212171.81i22.2466.971.586i2466.971.585i
6.11942157.49i11942157.49i23.3552.585.902i3552.585.903i
7.9869.1143.17i9869.1143.17i24.4835.6100.22i4835.6100.22i
8.7993.9128.85i7993.9128.85i25.6316.0114.54i6316.0114.54i
9.6316.0114.54i6316.0114.54i26.7993.9128.85i7993.9128.85i
10.4835.6100.22i4835.6100.22i27.9869.1143.17i9869.1143.17i
11.3552.585.902i3552.585.903i28.11942157.49i11942157.49i
12.2466.971.586i2466.971.585i29.14212171.81i14212171.81i
13.1578.657.267i1578.657.268i30.16679186.12i16679186.12i
14.887.7542.983i887.7542.951i31.19344200.44i19344200.44i
15.394.1728.149i394.2628.634i32.22206214.76i22206214.76i
16.99.46719.816i98.17714.317i33.25266229.07i25266229.07i
17.1.51017+7.7256i

### Table 2.

Approximations of the eigenvalues for the discrete and “continuous” operators (K<1).

k2=0, EI=1, ρ=1, N=64, εf=1030
No.k1=1+104k1=1+107k1=1+1010
1.222.22+155.56i176.06+264.07i106.90+372.41i
2.133.38+124.45i87.723+211.27i2.43781016+254.44i
3.64.540+93.396i3.96091018+162.37i3.07171017+123.66i
4.9.8081+58.559i5.79771019+116.23i1.30121017+4.9349i
5.7.77251021+5.0488i6.36611020+44.182i8.92321018+44.421i
6.4.93651021+38.994i6.00661020+4.9383i8.91431018+44.406i
7.7.80071021+4.8257i6.01141020+4.9313i1.30121017+4.9347i
8.9.8081+58.559i6.68011020+44.651i2.99071017+123.09i
9.64.540+93.396i2.93811018+139.20i1.12861016+234.06i
10.133.38+124.45i87.723+211.27i3.22801016+315.13i
11.222.22+155.56i176.06+264.07i106.90+372.41i

### Table 3.

Eigenvalues closest to the imaginary axis as k11+.

k2=0, EI=1, ρ=1, N=64, εf=1030
#k1=1104k1=1107k1=11010
1.175.34+140.01i129.20+237.56i60.143+338.80i
2.96.394+108.84i50.206+186.90i8.6602+240.01i
3.36.896+78.778i8.0431+121.15i0.28608+123.37i
4.6.0769+42.171i0.23455+44.410i0.0074192+44.413i
5.0.11141+4.9324i0.0035253+4.9348i0.00011148+4.9348i
6.0.11141+4.9324i0.0035253+4.9348i0.00011148+4.9348i
7.6.0769+42.171i0.23455+44.410i0.0074192+44.413i
8.36.896+78.778i8.0431+121.15i0.28608+123.37i
9.96.394+108.84i50.206+186.90i8.6602+240.01i
10.175.34+140.01i129.20+237.56i60.143+338.80i

### Table 4.

Eigenvalues closest to the imaginary axis as k11.

### 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]

1201ϱxAxht2xtExIxhxx2xtdx,E1

where hxt is the transverse deflection, Ex is the modulus of elasticity, Ix is the area moment of inertia, ϱx is the linear density, and Ax 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

ϱxAxhttxt+ExIxhxxxtxx=0,0x1,t>0,E2

and the boundary conditions

h0t=hx0t=0andM1t=Q1t=0,E3

where Mxt and Qxt are the moment and the shear, respectively [10]:

Mxt=ExIxhxxxtandQxt=Mxxt.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

Ut=Q1tM1tTandYt=ht1thxt1tT,E5

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

Ut=KYt,E6

where K is the 2×2 feedback matrix. We select

K=codiagk2k1,k1,k20,E7

with k1,k2 being the control parameters. The feedback (6) can be written as

E1I1hxx(1t)=k1ht(1t)and(ExIxhxx(xt))xx=1=k2hxt1t.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 hx0=h0x,htx0=h1x.

Notice that the choice of a feedback matrix K defines whether the system is dissipative or not. Indeed, let Et be the energy of the system, defined by representation (1). Evaluating Ett on the solutions of Eq. (2) satisfying the left-end conditions from Eqs. (3), we obtain

Ett=01ϱxAxhtxthttxt+ExIxhxxxthxxtxtx=(ExIxhxx(xt))xht(xt)|x=1+E1I1hxx1thxt1t.E9

Taking into account Eqs. (4) and (6), we represent the right-hand side of Eq. (9) as the dot product in R2:

Ett=Q1tht1t+M1thxt1t=UtYt=KYtYt.E10

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

Ett=0k2k102ht1thxt1t2ht1thxt1t=k1+k2ht1thxt1t.E11

Thus the system is not dissipative for all nonnegative values of k1 and k2.

## 2. Operator form of the problem

In what follows, we incorporate the cross-sectional area Ax into the density, write ρx instead of ϱxAx, and also abbreviate EIxExIx. Let H be the Hilbert space of two-component vector functions Ux=u0xu1xT equipped with the following norm:

UH2=1201EIxu0x2+ρxu1x2x.E12

Assuming that EI,ρC201 are positive functions, we obtain that the closure of smooth functions Ux=u0xu1xT satisfying u00=u00=0 will produce the energy space H=H0201×L201. Here H0201=uH201:u0=u0=0, 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:

Utxt=iLk1,k2UxtandUx0=u0xu1xT,E13

where 0x1,t0. The dynamics generator Lk1,k2 is given by the following matrix differential expression:

Lk1,k2=i011ρx2x2EIx2x20,E14

defined on the domain

DLk1,k2={U=u0u1TH:u0H401,u1H0201;u10=u10=0;EI1u01=k1u11EIxu0xx=1=k2u11}.E15

For any k1k2R2, the adjoint operator Lk1,k2 [13] is given by

Lk1,k2=Lk2,k1,E16

i.e., Lk1,k2 is defined by the same differential expression (14) on the domain described in Eq. (15), where k1 and k2 are replaced by k2 and k1, respectively. It follows from Eq. (16) that L0,0 is self-adjoint in H and thus L0,0 is the dynamics generator of the clamped-free beam model. For the reader’s convenience, we summarize the properties of Lk1,k2 from [1, 2] needed for the present work.

Proposition 1:

1. Lk1,k2 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 k1k2R2,k1+k2>0, the operator Lk1,k2 is a rank-two perturbation of the self-adjoint operator L0,0 in the sense that the operators Lk1,k21 and L0,01 exist and are related by the rule

Lk1,k21=L0,01+Tk1,k2,E17

where Tk1,k2 is a rank-two operator. A similar decomposition is valid for the adjoint operator, i.e.,

Lk1,k21=L0,01+Tk1,k2,Tk1,k2=Tk2,k1.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 k1,k2>0 and k1k2EIρ. Let

K=k1+k2Ak1k2/A,A=EIρ,andK1.E19

The following asymptotic approximations for the eigenvalues λn (as n) of the operator Lk1,k2 hold:

1. If 1<K<, then for n one has

λn=signKnEIρπn214ln2K+1K1+iπnlnK+1K1+Oneπn.E20

2. If 0<K<1, then for n one has

λn=signKnEIρ2n+12π214ln2K+1K1+2n+12lnK+1K1+Oneπn.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 ε>0, one can find a positive integer N, such that all eigenvalues λn with nN+1 satisfy the estimate

λnsignKnEIρπn214ln2K+1K1+iπnlnK+1K1εE22

for the case when 1<K< and

λnsignKnEIρ2n+12π214ln2K+1K1+2n+12lnK+1K1εE23

for the case when 0<K<1. The following important question holds: From which indexNcan 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 01, 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 k1,k2-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).

An eigenvalue λn of the dynamics generator Lk1,k2 is called a deadbeat mode if λn=iβn,βn>0. If the corresponding eigenfunction is Φnx, then the evolution problem (13) has a solution given in the form eiλntΦnx=eβntΦnx, 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 , with β>0. To deal with the deadbeat modes analytically, we rewrite the spectral equation Lk1,k2Φx=λΦx in the form of an equivalent problem for an operator pencil [17] as

EIφx=λ2ρφx,φ0=φ0=0,EIφ1=k1φ1,EIφ1=k2φ1.E24

If λn and φnx are an eigenvalue and eigenfunction of the pencil (24), then λn is also an eigenvalue of Lk1,k2 with the eigenfunction Φnx=1iλnφnxφnxT.

To solve problem (24), we first redefine the spectral and control parameters to eliminate ρ and EI from Eq. (24). We define λ˜,k˜1, and k˜2 by λ=EI/ρλ˜ and k˜j=EIρkj, j=1,2. Substituting these relations into Eq. (24) and eliminating the “tilde,” we obtain the following Sturm-Liouville eigenvalue problem:

φx=λ2φx,φ0=φ0=0,φ1=k1φ1,φ1=k2φ1.E25

The solution of Eq. (25) satisfying the left-end boundary conditions φ0=φ0=0 can be written in the form

φλx=Aλcoshλxcosλx+Bλsinhλxsinλx.E26

Substituting formula (26) into the right-end boundary conditions of Eq. (25), one gets a system for the coefficients Aλ and Bλ:

Aλ1+ik1coshλ1ik1cosλ+Bλ1+ik1sinhλ1ik1sinλ=0,Aλ1ik2sinhλ1+ik2sinλ+Bλ1ik2coshλ+1+ik2cosλ=0.E27

Let Δλ be the determinant of the matrix of coefficients for Aλ and Bλ in Eqs. (27). System (27) has nontrivial solutions if and only if Δλ=0, i.e.,

1+k1k2+1k1k2coshλcosλ+ik1+k2sinhλsinλ=0.E28

Theorem 1: The following results hold in the case when k1>0 and k2=0. Similar results hold in the case when k1=0 and k2>0.

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

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

λn=μn2,μn=xn1+i,xn=π22n+1,n=0,1,2,.E29

3. For any k1>1, there exist a finite number Nk1 of deadbeat modes. Each mode has the form λ=μ2,μ=x1+i, where x is a root of the function

Hxk12+1+k1cos2x+1k1cosh2x,x>0.E30

Let Xk1=12πcosh11+4k11, and then Nk1 satisfies the estimate

2Xk1+1Nk12Xk1+3.E31

Hence Nk1 as k11+. (By X we denote the greatest integer less than or equal to X).

Proof: Let μ=λ=x1+i,x>0. Taking into account the relations

2coshμcosμ=cosh1+iμ+cosh1iμ,2isinhμsinμ=cosh1+iμcosh1iμ,

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

2+1+k1cos2x+1k1cosh2x=0,x>0.E32

It can be readily seen that if 0<k1<1, then 2+1+k1cos2x>0 and 1k1cosh2x>0, which means that Eq. (32) has no solutions. Statement (1) is shown. Statement (2) follows immediately if one considers Eq. (32) for k1=1.

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

cosh2x=2k11+1+2k11cos2x,x>0,k1>1.E33

The left-hand side of Eq. (33) is monotonically increasing, while the right-hand side is sinusoidal, with maximum 1+4k11 and minimum 1, and period π. So the graphs of the left- and right-hand side have intersections only on the interval 012cosh11+4k11. 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 k1, the two graphs can have two intersections, one tangential intersection or no intersections on the second half-period. This leads to estimate (31).

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

## 4. Structure of the deadbeat mode set

The main result on the existence and distribution of double roots of the function Hxk1 is presented in the statement below.

Theorem 2: For a given k1>1, the multiplicity of each root of Hxk1 does not exceed 2. There exists a sequence k1nn=0,1,2, such that the function Hxk1 has a double root if and only if k1=k1n for some n. So the original spectral problem with k1=k1n,k2=0 has a double deadbeat mode λn=μn2=2ixn2. The following asymptotic formulas hold

xn=3π4+πn+Pn1+OPn2,wherePn=exp3π2+2πn,E34

and

k1n=1+4Pn1+OPn2.E35

Proof: If x is a double root of H, then Hxk1=Hxk1=0, i.e., separating the real and imaginary parts, we have

2+k1+1cos2xk11cosh2x=0,E36
k1+1sin2x+k11sinh2x=0.E37

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

Gx1+cos2xsinh2x+1+cosh2xsin2x=0,x>0.E38

Rewriting Eq. (37) in the form k1+1sin2x=k11sinh2x, and taking into account that k1>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.,

π2n+1<2x<3π2+2πn,n0,1,2,

Eq. (38) does not have any solutions. Indeed, in the above range of x, we have cos2x+sin2x=2sin2x+π/4<1 and sin2xcos2x<1. With such estimates we obtain that

Gx=sin2x+sinh2x+12e2xcos2x+sin2x+12e2xsin2xcos2x<sin2x<0,E39

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

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

3π2+2πn<2x<2πn+1,n0,1,2.

It is convenient to rewrite system given by (36) and (37) in the form 2x=3π/2+2πn+s, where n0,1,2, 0<s<π/2. If gsG3π/4+πn+s, then Eq. (38) generates the following equation for g:

gs1+sinssinh3π2+2πn+s1+cosh3π2+2πn+scoss=0.E40

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

g0=sinh3π2+2πn1cosh3π2+2πn<0,

and for s=π/2 we have gπ/2=2sinh2πn+1>0. Evaluating g we have

gs=sins+1+2sinscosh3π2+2πn+s>0.E41

Thus gs is a monotonically increasing function, such that g0<0<gπ/2, which means that g has a unique root on 0π/2.

Finally we show that the multiplicity of a multiple root cannot exceed 2. Using a contradiction argument, assume that there exists x0, such that in addition to Eqs. (36) and (37), one has Hx0k1=0, i.e., the multiplicity of x0 is at least 3. The system Hx0k1=0 and Hx0k1=0 can be written as

k1+1sin2x0+k11sinh2x0=0,k1+1cos2x0+k11cosh2x0=0.E42

Since k1>1, the second equation of (42) yields cos2x0<0. Also, since x0 is a multiple root, we must have sin2x0<0. Then 2x0 is in the third quadrant, which means that Gx00, as we have seen above. This contradicts our assumption that x0 is a root of Eqs. (36) and (37).

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

sin2Pn1=2Pn1+OPn3,cos2Pn1=12Pn2+OPn4,2sinh3π2+2πn+2Pn1=Pn+2+Pn1+OPn2,2cosh3π2+2πn+2Pn1=Pn+2+3Pn1+OPn2.E43

Evaluating gs from Eq. (40) for s=2Pn1 and using Eq. (43), we get

g2Pn1=1+sin2Pn1sinh3π2+2πn+2Pn11+cosh3π2+2πn+2Pn1cos2Pn1=2Pn1+OPn2.E44

Representation (44) implies that there exists n0, such that for all nn0, we have g2Pn1>0. Taking into account that g0<0, we obtain that the root sn, nn0, of the function gs is located on the interval 02Pn1. To find the location of this root more precisely [18], we use linear interpolation. Namely, substituting Eq. (43) into the expression for gs from Eq. (41) yields

g2Pn1=Pn2+O1.E45

Replacing gs by the linear function tangential to gs at the point 2Pn1g2Pn1, and finding the root of this function, we get

sn=2Pn1g2Pn1g2Pn1+OPn2=2Pn1+OPn2.E46

Having this approximation for sn, we immediately get

xn=3π4+πn+sn2=3π4+πn+Pn1+OPn2.E47

From the equation Hxnk1n=0, we obtain the formula for k1n as

k1n=sinh3π/2+2πn+sn+cossnsinh3π/2+2πn+sncossn.E48

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

Corollary 1: Let k1=k1n for some n+0, and let xn be the corresponding double root of the function Hxk1n. Then λ0=2ixn2 is an eigenvalue of the operator Lk1n,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

Lk1n,0Φ=λ0Φ,Lk1n,0Ψλ0Ψ=Φ.E49

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 k2=0 has two linearly independent solutions ψ and χ. Each function satisfies the problem

φx=λ02φx,φ0=φ0=0,φ1=0,φ1=iλ0k1φ1.E50

First we observe that ψ1χ10. Indeed, if ψ1=0, then we have

01ψσ2=01ψσψσ¯=λ0201ψσ2.E51

Since λ0 is purely imaginary, Eq. (51) is not valid. We define a new function:

gx=ψxψ1χ1χx.E52

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

φx=λ02φx,φ0=φ0=0,φ1=φ1=0,E53

and therefore

01gσ2dσ=01gσgσ¯σ=λ0201gσ2dσ.E54

Eq. (54) is valid if and only if λ02>0; however, for a deadbeat mode, λ02<0. The obtained contradiction means that for each double mode, there is one eigenfunction and one associate function.

### 4.1 Deadbeat mode behavior as k1→1

As k1 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 k1 approaches 1 from above or below.

As follows from Table 3, the real parts of the eigenvalues decrease steadily as k11+, to a point where the eigenvalue becomes a deadbeat mode. An increase in the number of deadbeat modes can be seen as k11+, 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 k11+. (Compare modes no.5 and no.7 for k11=104, modes no.4 and no.7 for k11=106, modes no.5 and no.8 for k11=108, and modes no.4 and no.7 for k11=1010). 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 k11. However the rate at which their real parts approach zero is significantly lower than in the case k11+. Even at k1=11010, the eigenvalue closest to the imaginary axis has a real part of about 104, 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 k11± (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 k11 and k11+. 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 Lk1,k2, 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

TNξ=cos,1ξ1whereξ=cosθ,θ0π.E55

The cardinal functions, ψkξ, and the Chebyshev-Gauss-Lobatto (CGL) grid points ξk are defined as follows:

ψkξ=1k1ξ2TN1ξckN12ξξk,ξk=cosk1πN1,for1kN,E56

where coefficients ck are such that c1=cN=2 and ck=1 for 1<k<N. The main property of cardinal functions is ψkξj=δkj (using the Kronecker delta). The family ψkk=1N forms a basis in the space of polynomials of degree N1, i.e., if f is such polynomial, then f and f can be written in the forms

fξ=k=1Nfξkψkξandfξ=k=1Nfξkψkξ.E57

If f=fξ1fξ2fξNT and g=fξ1fξ2fξNT, then g=Df, where D is the Chebyshev derivative matrix with the elements

D11=DNN=1+2N126,Dkk=ξk21ξk2for1<k<N,Dj,k=cj1j+kckξjξkforjk.E58

### 5.1 Discretization of Lk1,k2

Rescaling the independent variable x as ξ=2x1, we rewrite the operator and its domain, representations (14) and (15), in the form

Lk1,k2=i0116ρξ2ξ2EIξ2ξ20,E59

and

DLk1,k2={u0u1TH:u0H411,u1H0211;u11=u11=0;4EI1u01=k1u114EIξu0ξξ=1=k2u11},E60

where H=H0211×L211, equipped with the norm

UH2=141116EIξu0ξ2+ρξu1ξ2.E61

We approximate the action of Lk1,k2 on the finite-dimensional subspace HNH of polynomials of degree at most N1. Using the CGL grid and the cardinal functions, we substitute for u0 and u1 their truncated expansions:

u0ξk=1NΦkψkξ,Φk=u0ξk,u1ξk=1NΘkψkξ,Θk=u1ξk.E62

Let Φ and Θ be N-dim vectors and Ψ be a 2N-dim vector defined by

Φ=Φ1Φ2ΦNT,Θ=Θ1Θ2ΘNT,Ψ=ΦΘ.E63

Let L be the finite-dimensional approximation of the differential operator Lk1,k2. The discretized operator L induced by Lk1,k2 can be given by

L=i0IN×N16EIρD40,E64

where IN×N is the N×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

ΦN=0,DΦN=0,ΘN=0,DΘN=0,4EID2Φ1+k1Θ1=0,4EID3Φ1k2DΘ1=0.E65

Let rN,zN,lNRN be auxiliary row-vectors

rN=0001,zN=0000,lN=1000E66

and Djn designate the jth row of the nth derivative matrix Dn. Using Eqs. (66) we represent Eqs. (65) as the following matrix equation:

KΨrNzNDN1zNzNrNzNDN14EID12k1lN4EID13k2D11ΦΘ=0.E67

K is called the boundary operator. Let KN be the kernel of K, i.e., KN=vR2N:Kv=0. We have to identify all eigenvalues of the operator L, when its domain is restricted to KN. It is clear that KN is isomorphic to Rk with kdimKN=dimHNrankK=2N6. Let B be the matrix consisting of column vectors that form an orthonormal basis in KN. It is clear that BTB is the identity matrix on Rk and BBT 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 BTLB. However, the inverse statement is not necessarily true. Indeed, we observe that BBT is the identity in KN, which is not equivalent to the identity in HN. Assume now that λ is an eigenvalue of BTLB with corresponding eigenvector vRk. If Ψ=Bv, we have

BBTLΨ=BBTLBv=λBv=λΨ,E68

but BBTLL, which indicates that fake eigenvalues may exist.

### 5.3 Filtering of spurious eigenvalues

In order to decide which eigenvalues of BTLB should be discarded, we impose the following condition. Let Λ be the spectrum of BTLB and V be the set of its eigenfunctions. We construct the set of “trusted” eigenvalues [14, 15], for some εf>0 filtering precision, as

Λε=λΛ:LBvλλBvλC<εffor corresponding eigenvectorvλV,E69

where 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):

ΨC=π/4N1k=1N1ξk216EIξkD2Φk2+ρξkΘk2.

## 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, k1, is such that k11+ and the other one k2=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 k1nn=1 of the values of parameter k1, such that the corresponding deadbeat mode has a multiplicity 2, which yields the existence of the associate mode shapes for the operator Lk1,k2. The formulas for the double deadbeat modes and asymptotics for the sequence k1n as n have been derived.

## Acknowledgments

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

## References

1. 1. Shubov MA, Kindrat LP. Spectral analysis of the Euler–Bernoulli beam model with fully nonconservative feedback matrix. Mathematicsl Methods in the Applied Sciences. 2018;41(12):4691-4713. DOI: 10.1002/mma.4922
2. 2. Shubov MA, Shubov VI. Stability of a flexible structure with destabilizing boundary conditions. Proceedings of the Royal Society A. 2016;472(2191):20160109. DOI: 10.1098/rspa.2016.0109
3. 3. Shubov MA. Exact controllability of nonselfadjoint Euler-Bernoulli beam model via spectral decomposition method. IMA Journal of Mathematical Control and Information. 2008;25(2):185-203. DOI: 10.1093/imamci/dnm018
4. 4. Guiver CH, Opmeer MR. Non-dissipative boundary feedback for Rayleigh and Timoshenko beams. Systems and Control Letters. 2010;59(9):578-586. DOI: 10.1016/j.sysconle.2010.07.002
5. 5. Shubov MA. Exact controllability of coupled Euler-Bernoulli and Timoshenko beam model. IMA Journal of Mathematical Control and Information. 2006;23(3):279-300. DOI: 10.1093/imamci/dnm059
6. 6. Curtain RF, Zwart HJ. An Introduction to Infinite Dimensional Linear Systems Theory. New York, NY: Springer; 1995
7. 7. Dowell EH, editor. A Modern Course on Aeroelasticity. 4th Revised ed. Boston, MA: Klumer Academic Publishers; 2004
8. 8. Balakrishnan AV, Shubov MA. Asymptotic behaviour of aeroelastic modes for aircraft wing model in subsonic air flow. Proceedings of the Royal Society A. 2004;460(2044):1057-1091. DOI: 10.1098/rspa.2003.1217
9. 9. Shubov MA, Rojas-Araneza M. Four-branch vibrational spectrum of double-walled carbon nanotube model. Proceedings of the Royal Society A. 2011;467(2125):99-126. DOI: 10.1098/rspa.2010.0100
10. 10. Benaroya H. Mechanical Vibration: Analysis, Uncertainties, and Control. Upper Saddle River, NJ: Prentice Hall; 1998
11. 11. Han SM, Benaroya H, Wei T. Dynamics of transversely vibrating beams using four engineering theories. Journal of Sound and Vibration. 1999;225(5):935-988. DOI: 10.1006/jsvi.1999.2257
12. 12. Shubov MA. Generation of Gevrey class semigroup by non-selfadjoint Euler-Bernoulli beam model. Mathematicsl Methods in the Applied Sciences. 2006;29(18):2181-2199. DOI: 10.1002/mma.768
13. 13. Gohberg IT, Krein MG. Introduction to the Theory of Nonselfadjoint Operators in Hilbert Space. Translations of Mathematical Monographs No. 18. Providence, RI: AMS; 1996
14. 14. Mason JC, Handscomb DC. Chebyshev Polynomials. Boca Raton, FL: Chapman & Hall/CRC; 2003
15. 15. Trefethen LN, Bau D III. Numerical Linear Algebra. Philadelphia, PA: SIAM; 1997
16. 16. Shubov MA, Wineberg S, Holt R. Numerical investigation of aeroelastic mode distribution for aircraft wing model in subsonic air flow. Mathematical Problems in Engineering. 2010;2010:879519. DOI: 10.1155/2010/879519
17. 17. Marcus AS. Introduction to the Spectral Theory of Polynomial Operator Pencils. Translations of Mathematical Monographs No. 71. Providence, RI: AMS; 1988
18. 18. Fedoryuk MV. Asymptotic Analysis. Berlin, Heidelberg: Springer-Verlag; 1993
19. 19. Naimark MA. Linear Differential Operators. New York, NY: Ungar Pub Co; 1967

Written By

Marianna A. Shubov and Laszlo P. Kindrat

Reviewed: 20 May 2019 Published: 20 June 2019