Structured Approaches to General Inverse Eigenvalue Problems

Abstract

An inverse eigenvalue problem is one where a set or subset of (generalized) eigenvalues is specified and the matrices that generate it are sought. Many methods for solving inverse eigenvalue problems are only applicable to matrices of a specific type. In this chapter, two recently proposed methods for structured (direct) solutions of inverse eigenvalue problems are presented. The presented methods are not restricted to matrices of a specific type and are thus applicable to matrices of all types. For the first method, the Cayley–Hamilton theorem is developed for the generalized eigenvalue vibration problem. For a given (desired) frequency spectrum, many solutions are possible. Hence, a discussion of the required information and suggestions for including structural constraints are given. An algorithm for solving the inverse eigenvalue problem using the generalized Cayley–Hamilton theorem is then demonstrated. An algorithm for solving partially described systems is also specified. The Cayley–Hamilton theorem algorithm is shown to be a good tool for solving inverse generalized eigenvalue problems. Examples of application of the method are given. A second method, referred to as the inverse eigenvalue determinant method, is also introduced. This method provides another direct approach to the reconstruction of the matrices of the generalized eigenvalue problem, given knowledge of its eigenvalues and various physical parameters. As for the first method, there are no restrictions on the type of matrices allowed for the inverse problem. Examples of application of the method are also given, including application-oriented examples.

1. Introduction

Many physical problems can be rewritten as generalized inverse eigenvalue problems. For example, in mechanical and structural system design, engineers are often faced with the task of designing systems, which have natural frequencies either that must fall outside a specific range or that must be at exactly certain frequencies. These design problems can be considered as inverse eigenvalue problems because the eigenvalues are directly related to the natural frequencies of the system.

Much focus has been applied to the study of discrete inverse eigenvalue problems. This has been made clear by a thorough review of the topic by Chu and Golub [1, 2]. Gladwell [3] takes a more direct route in which he considers specific inverse problems and matrix structures related to mechanical vibrations. Particularly, it appears that most of the literature focuses on system identification. One of the most common techniques in inverse eigenvalue problems is to use/measure the system’s spectrum and then constrain the system in some fashion in order to obtain a second spectrum [4–8]. This clearly indicates that a system usually exists and that it can be tested to obtain data required for the inverse problem of mathematically reconstructing the system. Although interesting, this approach cannot be used for novel engineering design to construct a system having a specific spectrum without another system on which to base the design.

In most cases, the research on inverse eigenvalue problems has focused on the existence, uniqueness and computability of a solution. Other studies are typical variations on those described above, including partially described problems where not all spectral information is known. These types of problems have been considered in the studies of Gladwell and Willms [9] and Ram and Elhay [10] and are of interest for problems requiring only certain frequencies to be specifically determined. Interestingly, Dias de Silva [11] and de Oliveira [12] have shown that an n × n matrix always exists when a minimum of n − 1 prescribed matrix entries and a prescribed characteristic polynomial are given as design information. The results of Dias de Silva and de Oliveira guarantee existence but not uniqueness of the matrix.

One of the main shortcomings of current inverse eigenvalue theory is the lack of a solution for general matrices having predefined forms but which do not fit within current known solution methods. In this chapter, two recently proposed methods for structured (direct) solutions of inverse eigenvalue problems are presented. The presented methods are not restricted to matrices of a specific type and are thus applicable to general matrices. For the first method, the Cayley–Hamilton theorem is developed for the generalized eigenvalue vibration problem. The Cayley–Hamilton theorem algorithm is shown to be a good design tool for solving inverse generalized eigenvalue problems. Examples of application of the method are given. A second method, referred to as the inverse eigenvalue determinant method, is also introduced. This method provides another direct approach to the reconstruction of the matrices of the generalized eigenvalue problem, given knowledge of its eigenvalues and various physical parameters. As for the first method, there are no restrictions on the type of matrices allowed for the inverse problem. Examples of application of the method are also given, including application-oriented examples.

2. Problem statements

In this chapter, we consider the following problems:

PROBLEM A: Given a specified set of eigenvalues, λ₁, …, λ_n construct an n-order system, described by an n × n matrix A, which has λ₁, …, λ_n as its eigenvalues.

Although this problem has been considered for specific forms of matrices (i.e. Jacobi, band or other matrix forms), a general solution approach does not currently exist. Problem A leads to the generalized inverse eigenvalue problem, also presented here:

PROBLEM B: Given a specified set of eigenvalues, λ₁, …, λ_n construct a system with n-degrees of freedom, described by two n × n matrices (M and K), which has λ₁, …, λ_n as its generalized eigenvalues: det(K − λM) = 0 for the given λ₁, …, λ_n.

For an engineer, Problem B relates directly to the design problem stated earlier, where a conservative vibrating system having specified natural frequencies is sought.

PROBLEM C: Given a specified set of eigenvalues, λ₁, …, λ_m construct a system with n degrees of freedom, where m < n, described by two n × n matrices (M and K), which has λ₁, …, λ_m as its generalized eigenvalues: det(K − λM) = 0 for the given λ₁, …, λ_m.

Problem C, referred to as a partially described system, is one for which only a select few eigenvalues are specified, as opposed to the entire spectrum being given.

3. Cayley–Hamilton Approach for Inverse Eigenvalue Problems

3.1. Basic Theory

In order to solve Problem A, a discussion of the Cayley–Hamilton theorem is required. The Cayley–Hamilton theorem states that if p(λ) is the characteristic polynomial of a square matrix A, obtained from p(λ) = det(λI − A), then substituting A for λ in the polynomial gives the zero matrix. Thus, by applying the theorem, matrix A satisfies its own characteristic polynomial, p(A) = 0 [13].

The Cayley–Hamilton theorem can be useful in inverse eigenvalue problems beyond the typical statement that a square matrix satisfies its own characteristic equation. Once the characteristic polynomial of a system is found from the specified spectral data, the Cayley–Hamilton theorem can be used to find an unknown matrix A, which represents the system. A set of unknown entries of the matrix A can be found from the set of equations that arise from the Cayley–Hamilton theorem. For example, suppose that a 2 × 2 matrix A is given by the form:

$A=\left[\begin{array}{cc}\hfill a_{11}\hfill & \hfill a_{12}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill \end{array}\right]$E1

In order to solve the inverse eigenvalue problem, the entries of matrix A can be populated using the Cayley–Hamilton theorem in conjunction with the set of specified eigenvalues. In order to demonstrate the process, suppose that the specified eigenvalues are given by −2 and −3. The characteristic polynomial can then constructed as

$p\left(\lambda \right)=\left(\lambda +2\right)\left(\lambda +3\right)={\lambda }^2+5\lambda +6$E2

Using the Cayley–Hamilton theorem, λ in equation (2) is replaced by A from equation (1) so that

$p\left(A\right)=A^2+5A+6I=\mathbf{0}$E3

where I is the identity matrix and 0 is the zero matrix. Expressing equation (3) term-by-term gives four equations with four unknowns, specifically

$\left[\begin{array}{cc}\hfill a_{11}^2+a_{12}{a}_{21}+5a_{11}+6\hfill & \hfill a_{11}{a}_{12}+a_{12}{a}_{22}+5a_{12}\hfill \\ \hfill a_{11}{a}_{21}+a_{21}{a}_{22}+5a_{21}\hfill & \hfill a_{22}^2+a_{12}{a}_{21}+5a_{22}+6\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right]$E4

However, equation independence is unclear. A discussion on this point is found in section 4.1. Solving the set of equations in equation (4) leads to two solutions:

$\left[a_{11}=-a_{22}-5,a_{12}=-\frac{a_{22}^2+5a_{22}+6}{a_{21}}\right]$E5

Any values can be assigned to a₂₁ and a₂₂ in equation (5), and the matrix A will have the desired eigenvalues given in equation (2). Consequently, Problem A has been solved, although it is clear that many solutions exist. This solution is particularly useful in solving inverse eigenvalue problems because it gives a range of A matrix values for which a system produces the same eigenvalues. The factors limiting the solutions as given in equation (5) are then based on the physical limits and fixed parameters of the system. Physical limits and fixed parameters serve to produce a set of solutions that are mathematically as well as physically possible.

3.2. Generalized Cayley–Hamilton theorem

Problem B is related to Problem A but takes on a more general form, which is conducive to real physical systems since many physical problems can be written as generalized eigenvalue problems. Although the characteristic polynomial for the generalized eigenvalue problem is similar to the single matrix case, it now involves two matrices rather than one. Therefore, the lesser known generalized Cayley–Hamilton theorem must be used by Chang and Chen [14].

The generalized Cayley–Hamilton theorem is modified to include a two square matrices, K and M. The nomenclature is used as a reminder that in many engineering applications the matrices represent mass (M) and stiffness (K) matrices. For the generalized eigenvalue problem, the characteristic polynomial takes the form p(λ) = det(K − λM). By substituting K and M for λ into the characteristic polynomial, the following relationship must be satisfied.

$p\left(K,M\right)=c_n{\left(M^{-1}K\right)}^n+c_{n-1}{\left(M^{-1}K\right)}^{n-1}\dots +c_1\left(M^{-1}K\right)+c_{0}I=0$E6

where c_n is the coefficient of λⁿ in p(λ). Equation (6) is valid as long as M is non-singular.

If the matrices K and M commute (i.e. KM = MK), then the generalized Cayley–Hamilton theorem can be written as

$p\left(K,M\right)=c_n{K}^n+c_{n-1}{K}^{n-1}M\dots +c_{1}K{M}^{n-1}+c_0{M}^n=0$E7

where no other restrictions are placed on matrix M. Solving these equations leads to the solution of Problem B.

3.3. Numerical example

Once again, a simple numerical example is used to demonstrate concepts. Using the same desired eigenvalues as for the previous example gives the characteristic polynomial given in equation (2). This time, the modified generalized Cayley–Hamilton theorem is used to form an equation for unknown matrices K and M so that

$p\left(K,M\right)={\left(M^{-1}K\right)}^2+5\left(M^{-1}K\right)+6I=\mathbf{0}$E8

Assuming that the K and M matrices have a similar form to that of A in equation (1), then equation (8) can be expanded to give four equations (from the four matrix entries) with eight unknowns (four unknown entries in each of the two matrices). Solving these equations produces several results, one of which can be written as

$\left[\begin{array}{l}{k}_{11}=\frac{k_{12}{k}_{21}{k}_{22}+5k_{12}{k}_{21}{m}_{22}+6k_{12}{m}_{21}{m}_{22}+6k_{21}{m}_{12}{m}_{22}-6k_{22}{m}_{12}{m}_{21}}{\left(k_{22}+2m_{22}\right)\left(k_{22}+3m_{22}\right)}\\ m_{11}=\frac{k_{21}{k}_{22}{m}_{12}+k_{12}{k}_{22}{m}_{21}-k_{12}{k}_{21}{m}_{22}+5k_{22}{m}_{12}{m}_{21}+6m_{12}{m}_{21}{m}_{22}}{\left(k_{22}+2m_{22}\right)\left(k_{22}+3m_{22}\right)}\end{array}\right]$E9

Once again it becomes clear that only two solutions are dependent, and selecting any values for k₁₂, m₁₂, k₂₁, m₂₁, k₂₂, m₂₂, will result in K and M matrices that produce a solution with the desired eigenvalues.

Increasing the size of the square matrices increases the number of independent variables faster than it does the dependent variables. In other words, for an nth order system with n-specified eigenvalues, 2n² unknown variables are required to find K and M, but, as will be shown later in this chapter, the generalized Cayley–Hamilton theorem only produces n-independent equations.

3.4. Algorithm

In order to use the Cayley–Hamilton theorem as a tool to find the n dimensional generalized eigenvalue system problem based on knowledge of desired eigenvalues and/or physical parameters, 2n² − n additional pieces of design information are required in addition to the n-desired eigenvalues. An algorithm for solving Problem B is as follows:

ALGORITHM. Given the n-desired eigenvalues λ₁, …, λ_n of an nth order generalized eigenvalue system

1. Generate the characteristic polynomial:

   $p\left(\lambda \right)=\left(\lambda -{\lambda }_1\right)\dots \left(\lambda -{\lambda }_{n-1}\right)\left(\lambda -{\lambda }_n\right)=c_n{\lambda }^n+c_{n-1}{\lambda }^{n-1}+\dots +c_1\lambda +c_0$E101

2. Generate forms for the M and K matrices from physical parameters and using the additional 2n² − n pieces of design information, applying symmetry and any other techniques based on the discretized forward problem and leaving an n number of unknowns;

3. Generate the Cayley–Hamilton theorem equation:

   $p\left(K,M\right)=c_n{\left(M^{-1}K\right)}^n+c_{n-1}{\left(M^{-1}K\right)}^{n-1}\dots +c_1\left(M^{-1}K\right)+c_{0}I=0$E111

   or for commuting M and K matrices:

   $p\left(K,M\right)=c_n{K}^n+c_{n-1}{K}^{n-1}M\dots +c_{1}K{M}^{n-1}+c_{0}M=0$E121

4. Extract all n² equations from the Cayley–Hamilton equation in step 3;

5. Select n non-zero independent equations (the Cayley–Hamilton matrix diagonal entries work well);

6. Compute the n unknowns by solving the n-independent equations;

7. Insert the n computed values into their appropriate places in the M and K matrices;

8. Verify that a valid solution is obtained by calculating det(λM + K) = 0 and ensuring that the initially desired eigenvalues are obtained.

The output consists of the n matrix entries of the M and K matrices.

4. Analysis

In this section, results on the number of independent equations and required pieces of information used in the prior sections are discussed.

4.1. Information produced by the Cayley–Hamilton theorem

We show here that given n distinct eigenvalues for an nth order system, the Cayley–Hamilton theorem can produce at most n-independent equations, even though n² equations are produced through the entry-by-entry analysis. Let p(t) = c_ntⁿ + c_n − 1t^n − 1 + … c₁t + c₀ be the characteristic polynomial of an n × n matrix A. The Cayley–Hamilton theorem states that p(A) = c_nAⁿ + c_n − 1A^n − 1 + … c₁A + c₀ is the zero matrix, which gives n² equations. Now, suppose that the matrix A is diagonal and let the diagonal entries be λ₁, λ₂, …, λ_n. The characteristic polynomial is

$p\left(t\right)=\left(t-{\lambda }_1\right)\left(t-{\lambda }_2\right)\dots \left(t-{\lambda }_n\right)$E10

In this case, since A is diagonal, p(A) is also a diagonal matrix with exactly n equations. In fact, the ith diagonal entry is p(λ_i).

Now, consider the case that A is not diagonal. Given n distinct eigenvalues for A, then A is diagonalizable, so that A ' = P^− 1AP is diagonal for some invertible matrix P. The characteristic polynomial of A ' is the same as the characteristic polynomial p(t) of A. In fact, the diagonal matrix A ' is the trivial solution to our design problem, and the ‘structure’ of the problem is housed in the matrix P. It is known that [15]

$p\left(A\right)=Pp\left(A\text{'}\right)P^{-1}$E11

Equation (11) states that p(A) can be obtained by combining the equations contained in p(A '). At the same time, since A ' is itself diagonal, then matrix p(A ') contains exactly n equations. Thus, this states that p(A) is a combination of exactly n-independent equations and we can expect that although p(A) has n² equations, only n of them are independent.

It is also important to point out that the n-independent equations obtained from the Cayley–Hamilton theorem are each nth order polynomials. For a polynomial system with n unknowns and also n equations, then Bézout’s theorem states that the problem has nⁿ complex solutions [16].

4.2. Required information

It is clear that the spectral information (eigenvalues) alone is not enough to solve the problem. The method presented, along with all other methods, is limited by the fact that an nth order system can produce at most n-independent equations, even though n² equations are obtained through applying the generalized Cayley–Hamilton theorem, as discussed above. If the matrices are completely unknown, then there may be as many as 2n² unknown entries in the M and K matrices. The Cayley–Hamilton method can produce at most n-independent variables and the remaining equations must be specified in other ways. For creating physically realistic systems, this generally entails preconditioning or constraining the structure of the matrices.

For solving discrete conservative vibration problems, Equation (6) can be used as long as the M matrix is non-singular, and equation (7) can be used if M and K commute. It becomes very clear by considering the example in equation (9) that more information is required in order to build a suitable vibrating system based on M and K. From an engineer’s point of view, any system, which can produce similar eigenvalues, has the potential of being a suitable solution, as long as it fits within the physical criteria set at the outset of the project. It is obvious that the example in Equation (9) can produce an infinite number of possible solutions. Equation (9) represents one of the several solutions to the Equation 8. As the system’s order increases, so does the number of potential solutions (Bézout’s theorem). Thus, engineers have many solutions at their disposal for creating suitable and optimized designs.

The inverse problem is then not limited by the eigenvalues, and in fact the eigenvalues alone do not contain enough information from which to build a physical system. Therefore, other information is required to solve the problem.

4.3. Structural constraints

One of the easiest ways to limit the number of matrix entries and include structural constraints is by including zero entries or by incorporating symmetry into the matrices. These constraints are interesting because they follow directly from real systems. For example, for many simply connected spring-mass systems, only diagonal entries are present in the mass matrix, and the stiffness matrix is expected to be symmetric and therefore can be mathematically forced to be symmetric by imposing that symmetric off-diagonal terms be the same. Physically, the structure of these matrices reflects the fact that the degrees of freedom of the system are coupled through stiffness element couplings only and no mass coupling is present. Furthermore, the symmetry of the stiffness matrix is a consequence of Newton’s third law. The system is further constrained by the fact that the off-diagonal terms in the stiffness matrix are not independent, but are related to the diagonal terms, once again reducing the number of unknowns. In this case, incorporating physical constraints into the mathematical structure of the problem can reduce the number of unknowns from 2n² to n².

In this case, the engineer still has the freedom to choose from several systems, which would satisfy the requirements. Further constraints may come in the form of available components, such as stiffeners, which must fit within specified physical dimensions or the financial budget.

These observations highlight the importance of the forward problem. Continuous systems are typically discretized using various methods, which produce a model of the system in matrix form. The form of the matrix is heavily dependent on the choice of discretization method. The inverse problem is consequently affected because it seeks to create a matrix that matches the form as chosen at the outset by the choice of discretization approach. Although the Cayley–Hamilton theorem does not discriminate in its ability to solve these various matrix forms, it is possible that certain discretization methods lead to simpler forms or matrices containing fewer variables. In this sense, an inverse problem may be easier to fully define.

4.4. Partially described systems

Another aspect that affects the amount of information required is the information available at the outset of the problem. Thus far, the entire spectral set has been used as input information, as well as specification of matrix entries when necessary. However, as stated in the study of Chu [1], often only portions of the entire spectrum are available. This is termed as a partially described inverse eigenvalue problem. This is true whether the subset of available spectral data arises from the design requirements or from the experimental data obtained for system identification. Regardless of the information available, the Cayley–Hamilton theorem can be used to produce n pieces of information for an nth order system stemming from n degrees of freedom. If certain eigenvalues are missing, the Cayley–Hamilton theorem can still be used.

Typically, in order to completely solve an nth order inverse problem using the Cayley–Hamilton theorem, only n unknown values should be present in the problem, regardless of the manner of their appearance along the solution path. Problem C makes full use of this detail during its solution by preconditioning the unknown matrix to account for this. For example, for a 3DOF system, it would be expected that only three unknowns should be present. Consider for example, the following system that has three known eigenvalues λ₁ = 0.47, λ₂ = 4.66, λ₃ = 10.87 and three unknown entries in the matrix given by

$A=\left[\begin{array}{ccc}\hfill a_1+a_2\hfill & \hfill -a_2\hfill & \hfill 0\hfill \\ \hfill -a_2\hfill & \hfill a_2+a_3\hfill & \hfill -a_3\hfill \\ \hfill 0\hfill & \hfill -a_3\hfill & \hfill a_3\hfill \end{array}\right]$E12

Using the Cayley–Hamilton theorem as described in Section 3, the matrix can be determined such that a₁ = 9.43, a₂ = 1.23, a₃ = 2.06, which gives a matrix solution for the system as

$A=\left[\begin{array}{ccc}\hfill 10.66\hfill & \hfill -1.23\hfill & \hfill 0\hfill \\ \hfill -1.23\hfill & \hfill 3.28\hfill & \hfill -2.06\hfill \\ \hfill 0\hfill & \hfill -2.06\hfill & \hfill 2.06\hfill \end{array}\right]$E13

As in Problem C, it is often the case that only partial spectral information is available. Therefore, the Cayley–Hamilton theorem can be used if the amount of missing spectral information is replaced by the same amount of matrix entry information. So, if only λ₃ = 10.87 is known then two of the three a matrix entries must be known. If a₂ = 1.23 and a₃ = 2.06 are known in the example above, then solving the Cayley–Hamilton equations gives a₁ = 9.43, λ₁ = 0.47, λ₂ = 4.66. Consequently, solving Problem C is not much different from solving Problem A. The same is true if the system is made up of a matrix M and a matrix K, except in this case, the generalized Cayley–Hamilton theorem must be used.

It becomes clear that one of the largest difficulties in using inverse eigenvalue theory to solve many physical problems is setting up the appropriate matrix form of the problem. This is achieved by consideration of the modelling approach chosen to model the forward problem. In most cases, the forward solution will utilize some form of discretization, such as finite elements, global elements, finite difference, and so on. The method chosen has a large impact on the structure of the matrix, making the solution easier or more difficult depending on the situation. Also, in discretizing, the method chosen to relate material parameters has a large effect on the number of independent variables. Therefore, it is extremely important to ensure that simplification methods are properly considered.

4.5. Polynomial solvers

The characteristic equation for an nth order system is an nth order polynomial. As discussed above, the Cayley–Hamilton theorem will produce n-independent equations. For both partially described and fully described systems, a system of nth order polynomials will need to be solved for the unknown quantities. This will require a good polynomial solver or polynomial root-finding algorithm. For the problems presented in the rest of this chapter, the symbolic computer algebra software Maple by Maplesoft, Waterloo, Canada, was used to solve the problems although any modern numerical or symbolic software package could similarly be employed. For our simulations, the Maple built-in function fsolve was used, as well as the freely available package DirectSearch. The details of the DirectSearch package can be found in the study of Moiseev [17]. It was observed that DirectSearch produced many more (valid) solutions than did fsolve (recall that a unique solution to most of these problems does not exist and thus multiple solutions are expected). The existence of multiple solutions is desirable from an engineering design point of view in that it implies that there are many solutions possible and it is up to the designer to choose the most physically realisable one, within the problem constraints. For the Cayley–Hamilton (and determinant method presented in the rest of the chapter) to be effective, a good polynomial solver is required so that the multiple possible solutions can indeed be found.

5. Example with the Cayley–Hamilton inverse method

In this section, we demonstrate by way of example how the Cayley–Hamilton method can be used to solve a practical engineering design problem by writing it as an inverse vibrations problem and thus as an inverse eigenvalue problem.

5.1. Problem statement

Given a desired fundamental frequency, construct a brace–plate system as described by a mass matrix M and a stiffness matrix K. All dimensional (geometric) properties of the brace–plate system are assumed to be specified and fixed except for the thickness of the brace h_c and the design variable, which must be found. The plate and brace material properties are also assumed to be known.

It is known from vibrations theory that the natural frequencies of a system specified by stiffness matrix K and mass matrix M are the square roots of the generalized eigenvalues of the matrices (K, M). Hence, a specified fundamental frequency of the system is the same as specifying the smallest generalized eigenvalue of (K, M).

5.2. Forward model

The first step in the solution process is to generate a finite-dimensional model of the continuous problem. The model is based on an orthotropic plate structurally reinforced by a brace in the weaker plate direction. The model is shown in Figure 1.

Figure 1.

Orthotropic plate reinforced with a rectangular brace.

The forward model is discretized using the assumed shape method. The assumed shape method is an energy method, which uses global plate elements with the kinetic and strain energy plate equations in order to determine the system’s equations of motion, from which the mass and stiffness matrices are extracted [18]. For the details of the development of the large mass and stiffness matrices, the reader is referred to the study of Dumond and Baddour[19]. The system is assumed simply supported, conservative and the material properties are assumed orthotropic. The forward model is created assuming that the mechanical properties are all related to Young’s moduli in the y direction.

5.3. Inverse model

The goal is to reconstruct the brace–plate system from a desired fundamental frequency. The generalized Cayley–Hamilton theorem inverse eigenvalue method is used as explained in Section 3.2.

A cross section of the fundamental modeshape is shown in Figure 2. It is clear that the brace affects the maximum amplitude of this modeshape, thus also affecting the associated frequency. In order to adjust the fundamental frequency of the brace–plate system to a desired value, it is necessary to adjust the thickness of the brace [20].

Figure 2.

Cross section of the brace–plate system’s fundamental modeshape.

5.4. Modelling considerations

The mechanical properties vary based on E_y, the Young’s modulus in the y direction, and the brace thickness controls the brace–plate system’s fundamental frequency. Therefore, the forward model is created using the assumed shape method while leaving these two parameters as variables. Thus, the mass matrix M is a function of h_c, the height of the brace–plate system at the (assumed fixed) location of the brace, and the stiffness matrix K is a function of h_c and also E_y. Here, 2×2 trial functions are used in the assumed shape method, so it follows that fourth order square matrices are created. The trial functions used are those of the simply supported rectangular plate given by

$w\left(x,y,t\right)={\displaystyle \sum _{n_x=1}^{m_x}{\displaystyle \sum _{n_y=2}^{m_y}sin\left(\frac{n_x\pi x}{L_x}\right)sin\left(\frac{n_y\pi y}{L_y}\right)}}{q}_{n_x{n}_y}\left(t\right)$E14

where m are the modal numbers, q is the time function and w is the displacement variable normal to the plate. The displacement variable w is then used directly in creating the kinetic and strain energy equations of the simply supported rectangular plate. These equations are broken into three sections as shown in Figure 1 in order to take into account presence of the brace. This procedure is described in detail in the study of [19]. It is assumed that the E_y is known and used as input information into the stiffness matrix. This leaves h_c as the only unknown parameter, appearing in both the mass and stiffness matrices.

In order to solve these matrices from the desired fundamental frequency, we must first create the characteristic polynomial using the desired frequency,

$p\left(\lambda \right)=\left(\lambda -a\right)\cdot \left(\lambda -b_1\right)\cdot \left(\lambda -b_2\right)\cdot \left(\lambda -b_3\right)$E15

where a is the specified fundamental frequency, and b₁, b₂ and b₃ are unknown values, which will also need to be found as part of solving the problem. Since we have assumed 2×2 trial functions, so that the mass and stiffness matrices are both 4×4, the characteristic polynomial must be fourth order, as shown in equation (15). Subsequently, p(λ) is expanded so that the polynomials coefficients can be found. Once the polynomial is created, the Cayley–Hamilton equation can be written by substituting (M^-1K) for λ into equation (15):

$p\left(K,M\right)=c_4{\left(M^{-1}K\right)}^4+c_3{\left(M^{-1}K\right)}^3+c_3{\left(M^{-1}K\right)}^2+c_1\left(M^{-1}K\right)+c_{0}I=0$E16

where c_n are the coefficients of λ in p(λ) determined through equation (15). Equation (16) produces 16 equations, of which only 4 are independent. According to Bézout’s theorem, selecting and solving the equations on the main diagonal for the four unknowns (h_b, b₁, b₂ and b₃) produce 4⁴=256 possible solutions. From the set of all possible solutions, complex-valued solutions can be immediately eliminated as not being physically meaningful. Clearly, further constraints must be added to the solution in order to get a solution which fits within the desired physical limits. These physical limits are based on the maximum and minimum brace dimensions, which are required to compensate for the range of plate stiffnesses used during the analysis, as well as the range of natural frequencies, which can be obtained using these system dimensions. Thus, the following constraints are implemented into the solution:

$\begin{array}{l}0.013\le h_b\le 0.016\mathrm{m}\\ 1\times {10}^7\le b_1\le 9\times {10}^8\mathrm{rad}/\mathrm{s}\\ 1\times {10}^7\le b_2\le 9\times {10}^8\mathrm{rad}/\mathrm{s}\\ 1\times {10}^7\le b_3\le 9\times {10}^8\mathrm{rad}/\mathrm{s}\end{array}$E17

Solving the four equations obtained from equation (16) within the constraints provided by equation (17) yields, a physically realistic solution satisfies the desired fundamental frequency, as well as the system’s parameters.

5.5. Results

The material properties of the brace and plate used during the analysis are given in Table 1.

Material properties	Values
Density, μ (kg/m^3)	403.2
Young’s modulus, Ey (MPa)	850
Young’s modulus, Ex (MPa)	Ey/0.078
Shear modulus, Gxy (MPa)	Ex×0.064
Poisson’s ratio, νxy	0.372
Poisson’s ratio, νyx	νLR×Ey/Ex

Table 1.

Material properties.

The dimensions used for the model throughout the analysis of the brace–plate system are shown in Table 2.

Dimensions	Values
Length, Lx (m)	0.24
Length, Ly (m)	0.18
Length, Lb (m)	0.012
Reference, x1 (m)	Lx/2–Lb/2
Reference, x2 (m)	x1+Lb
Thickness, hp (m)	0.003
Thickness, hb (m)	0.012
Thickness, hc (m)	hp+hb

Table 2.

Dimensions of brace-plate model.

These dimensions refer to those shown in Figure 1, where ‘p’ refers to the plate’s dimensions, ‘b’ refers to the brace’s dimensions and ‘c’ refers to the dimensions of the combined system.

As a basis for comparison, a plate having a Young’s moduli of E_y=850 MPa to which a brace is attached with a combined brace–plate thickness of h_c=0.015 m is found to have a fundamental natural frequency of 687 Hz, as calculated using the forward model. The analysis is then performed using the inverse Cayley–Hamilton method. As E_y of the plate is varied, the thickness of the brace–plate section is calculated such that the fundamental frequency of the brace–plate system is kept consistent at 687 Hz. The results of the computations are presented in Table 3.

Young’s modulus; ER (MPa)	Brace thickness; hc (m)	Fundamental frequency; a (Hz)
750	0.01576	687
800	0.01536	687
813	0.01527	687
850	0.01500	687
900	0.01466	687
950	0.01435	687

Table 3.

Results of the inverse model analysis.

Clearly, adjusting the thickness of the brace also has an effect on the other natural frequencies. These are presented in Table 4.

Young’s modulus; ER (MPa)	Brace thickness; hc (m)	b1 (Hz)	b2 (Hz)	b3 (Hz)
750	0.01576	774	1360	2653
800	0.01536	782	1363	2650
813	0.01527	784	1364	2650
850	0.01500	790	1366	2648
900	0.01466	798	1370	2645
950	0.01435	806	1374	2642

Table 4.

Calculated frequencies of the inverse model analysis.

Interestingly, the constraints indicated in equation (17), although physically strict, allow for more than one solution in certain cases. An example is shown in Table 5.

Young’s modulus; ER (MPa)	Brace thickness; hc (m)	a (Hz)	b1 (Hz)	b2 (Hz)	b3 (Hz)
750	0.01359	687	570	1149	2185

Table 5.

Alternate brace thickness solution satisfying the physical constraints.

5.6. Discussion

From these results, it is evident that it is possible to design a brace–plate system starting with a desired fundamental frequency in conjunction with using the proposed Cayley–Hamilton method. Table 3 clearly shows that by adjusting the thickness of the brace by small increments (10^-5 m, machine limit), it is possible to compensate for the variation in the cross-fibre stiffness (E_y) of the plate so that the fundamental frequency of the combined system is equal to that of the benchmark value of 687 Hz. The results obtained using the Cayley–Hamilton theorem algorithm match those values obtained using the forward model exactly. However, since no account has been taken of the other frequencies during the analysis, Table 4 shows that frequencies b₁–b₃ vary considerably from those values obtained for E_y=850 MPa. Therefore, it is important to ensure that there is a good understanding of what a given model can control. Moreover, it is interesting to note that within the strict physical constraints of equation (17), there is more than one brace–plate system (solution) that satisfies the Cayley–Hamilton theorem of equation (16). From Table 5, it can be seen that an alternate solution to the system exists, different from the one presented in Table 3, for a plate having an E_y of 750 MPa. In this case, by reducing the thickness of the brace, it is still possible to achieve a system having the desired frequency of 687Hz. However, the desired frequency is no longer the fundamental frequency but rather becomes the second frequency, and the fundamental has been replaced with a fundamental frequency of 570 Hz. It is important to keep this phenomenon in mind while designing a system. This is especially true if the order in the spectrum of a certain frequency associated with a certain modeshape is absolutely critical

6. The determinant method for partially described systems

Although effective, the Cayley–Hamilton method for solving inverse eigenvalue problems described in the previous section may not be the most computationally efficient method when solving partially described systems. Partially described systems are those for which only a select few eigenvalues are specified, as opposed to the entire spectrum being given. In this section, the determinant method, which first appeared in the study of Mir Hosseini and Baddour [21], is introduced.

For clarity, the method is explained for the case where two eigenvalues are required to be a certain value. The method proceeds in the same manner where the number of specified eigenvalues is larger. Solving the characteristic equation of the system requires the use of the determinant equation

$p\left(t\right)=det\left(K-\lambda M\right)=0$E18

where λ is the eigenvalue obtained by finding the roots of the polynomial equation 18. Suppose that two of the possible n roots (in this case the eigenvalues a₁ and a₂) are required to be of a certain value. Replacing λ in the characteristic polynomial with the desired a₁ and a₂ each in turn leads to two scalar equations, specifically

$\begin{array}{l}p\left(a_1\right)=det\left(K-a_{1}M\right)=0\\ p\left(a_2\right)=det\left(K-a_{2}M\right)=0\end{array}$E19

The two scalar equations in equation (19) can be simultaneously solved for up to two unknowns. This greatly reduces the computational effort involved. Constraints can be used to create a physically realistic system. Similar to the case for the Cayley–Hamilton theorem, the equations will be polynomials and the use of a good polynomial solver is required to obtain the multiple solutions expected.

6.1. Algorithm

In the case of Problem C, m eigenvalues are given for an n dimensional system, so that n − m is the number of unknown or unspecified eigenvalues for the system. The output of the determinant method consists of the remaining n − m eigenvalues as well as any of the unknown matrix entries.

Determinant Algorithm. A certain set of desired eigenvalues, a₁, …, a_m, which is smaller than the dimension of the system, are given as input information.

1. Generate the M and K matrices from the chosen forward modelling technique leaving variable m as unknown modelling parameters.

2. Generate a determinant equation from each desired frequency:

   p(a₁) = det(K − a₁M) and p(a₂) = det(K − a₂M)…p(a_m) = det(K − a_mM)where a₁, …, a_m are obtained from the desired natural frequencies;

3. Solve the set of m equations in m unknowns;

4. Insert the unknown values found from step 3 into the matrices and calculate the system natural frequencies using the forward modelling process in order to verify that the solution is consistent with the input values.

7. Example of using the determinant method

In this section, an example of how the determinant method can be used to solve an inverse vibration problem is presented. In this section, the authors demonstrate a method to impose two desired natural frequencies on a dynamical system of a simply supported Euler–Bernoulli beam to which a single lumped mass is attached by determining the magnitude and mounting position of the mass. The full details of this example are given in the study of Mir Hosseini and Baddour [21].

7.1. Modelling

7.1.1. The method of assumed modes

The method of assumed modes was used to find a finite-dimensional model (discrete model) for the beam with attached elements. As is well known for the assumed modes approach, the greater the number of modes used, the more accurate the results obtained. Usually, the vibrational modes of a related but simpler problem are used to find approximate solutions to a more complicated problem, in this example, the modeshapes of the beam with no attached elements are used. A good introduction to the assumed modes can be found in the study of Meirovitch [22].

The assumed method, when applied to a continuous, conservative vibrational system, will result in mass and stiffness matrices that represent the system. The dimensions of these matrices are determined by the degree of discretization selected for the problem. Here lies the main advantage of the assumed mode method over the popular finite element analysis (FEA). It has been shown in the study of Cha [23] that the same level of accuracy can be reached by the assumed modes method using smaller degrees of discretization than with FEA. This implies mass and stiffness matrices that are smaller. This reduced dimensionality of the system is far more relevant to solution of the inverse problem than to the forward problem.

The assumed modes method was chosen to derive the equations of motion for the case of an Euler–Bernoulli beam to which a number of discrete elements are attached. The transverse vibrations of the beam can be written as

$w\left(x,t\right)={\displaystyle \sum _{j=1}^N{\varphi }_j\left(x\right)}{\eta }_j\left(t\right)$E20

Here, w(x, t) is the transverse displacement of the beam, ϕ_j(x) is a space-dependent eigenfunction, η_j(t) is the generalized coordinate and N is the number of assumed modes chosen for the problem. It is important to note that ϕ_j(x) varies with the choice of the beam and any ϕ_j(x) should be chosen to satisfy the required geometric boundary conditions of the selected beam.

7.1.2. Derivation of equations of motion

In order to derive the equations of motion for the one-dimensional Euler–Bernoulli beam with multiple lumped point-mass attachments as shown in Figure 3, expressions for kinetic and potential energies must first be found. The kinetic energy of the beam is given by

Figure 3.

Beam with multiple lumped mass attachments.

$T=\frac{1}{2}{\displaystyle \sum _{j=1}^N{M}_j}{{\stackrel{.}{\eta }}_j}^2\left(t\right)+\frac{1}{2}{m}_1{\stackrel{.}{w}}^2\left(x_1,t\right)+....+\frac{1}{2}{m}_s{\stackrel{.}{w}}^2\left(x_s,t\right)$E21

where M_j is generalized masses of the bare beam (no attachments), an over dot indicates derivatives with respect to time and m₁…m_s are s lumped point masses positioned at x₁…x_s, respectively.

Using the same procedure, the equation for the potential energy can be written as

$V=\frac{1}{2}{\displaystyle \sum _{j=1}^N{K}_j{{\eta }_j}^2}\left(t\right)$E22

where K_j are the generalized stiffnesses of the bare beam.

Substituting equation (20) into equation (21), the following equation for kinetic energy is obtained

$T=\frac{1}{2}{\displaystyle \sum _{j=1}^N{M}_j{{\stackrel{.}{\eta }}_j}^2\left(t\right)}+\frac{1}{2}{m}_1{\left[{\displaystyle \sum _{j=1}^N{\varphi }_j\left(x_1\right){\stackrel{.}{\eta }}_j\left(t\right)}\right]}^2+....+\frac{1}{2}{m}_s{\left[{\displaystyle \sum _{j=1}^N{\varphi }_j\left(x_s\right){\stackrel{.}{\eta }}_j\left(t\right)}\right]}^2$E23

Equation (22) for the potential energy remains unchanged from a bare beam since no elastic element was added to the beam.

Having found the expressions for kinetic and potential energies in terms of ϕ and η , these are then substituted into the Lagrange’s equations to yield the equations of motion. Lagrange’s equations are given by

$\frac{d}{dt}\left(\frac{\partial T}{\partial {\stackrel{.}{\eta }}_i}\right)-\frac{\partial T}{\partial {\eta }_i}+\frac{\partial V}{\partial {\eta }_i}=0i=1,2,..,N$E24

where N corresponds to the number of generalized coordinates and hence the number of differential equations.

Substituting equations (23) and (22) into equation (24) and converting the system of equations into a matrix representation, the matrix equation of motion is given by

$\mathbf{M}\underset{¯}{\ddot{\eta }}+\mathbf{K}\underset{¯}{\eta }=\underset{¯}{0}$E25

where M and K are the system mass and stiffness matrices, respectively. The mass matrix is given by

$\mathbf{M}={\mathbf{M}}^d+m_1\cdot {\underset{¯}{\varphi }}_1\cdot {{\underset{¯}{\varphi }}_1}^T+\dots +m_s\cdot {\underset{¯}{\varphi }}_s\cdot {{\underset{¯}{\varphi }}_s}^T$E26

In equation (26), ${\underset{¯}{\varphi }}_1$…${\underset{¯}{\varphi }}_s$ are N-dimensional column vectors of the N eigenfunctions evaluated at the attachment points x₁ …x_s, so that for example

${\underset{¯}{\varphi }}_1=\left[\begin{array}{r}\hfill {\varphi }_1\left(x_1\right)\\ \hfill .\\ \hfill .\\ \hfill .\\ \hfill {\varphi }_N\left(x_1\right)\end{array}\right]$E27

M^d is a diagonal matrix whose diagonal components are the generalized masses M_i, and m₁… m_s are the masses of the lumped attachments.

For the stiffness matrix, since elastic elements are not added to the beam, the stiffness matrix remains a diagonal matrix whose elements are the generalized stiffnesses of the beam. Hence, the stiffness matrix is given by

$\mathbf{K}={\mathbf{K}}^d$E28

7.1.3. Choice of eigenfunctions

As can be seen in equations (26) and (27), eigenfunctions constitute a major component of the mass matrix of the combined system. It is solely a function of position, x and can take many different forms depending on the beam being considered. The most basic requirement of any eigenfunction simulating the vibrational behaviour of a beam is its ability to accommodate and satisfy the beam boundary conditions. Here, the eigenfunctions for an Euler–Bernoulli beam with simply supported boundary conditions are given by Cha [23]

${\varphi }_n\left(x\right)=\sqrt{\frac{2}{\rho L}}sin\left(\frac{n\pi x}{L}\right)$E29

where ρ is the mass per unit length of the beam and L represents the length of the beam.

7.1.4. Frequencies and modeshapes

In order to solve equation (25), a system of N second-order differential equations, the vector of generalized coordinates $\underset{¯}{\eta }$ is written as

$\underset{¯}{\eta }=\overrightarrow{\eta }{e}^{i\omega t}$E30

Here, ω is the frequency of vibration of the system. Also, the inclusion of the complex number “i” is justified given the fact that the system is conservative, and it is expected that the vibrations are purely oscillatory.

Substituting equation (30) into equation (25) and taking derivatives yields

$\left(-{\omega }^{2}M+K\right)\overrightarrow{\eta }=\underset{\_}{0}$E31

In order for equation (31) to have a non-trivial solution, the following equation must hold

$det\left(-{\omega }^2\mathbf{M}+\mathbf{K}\right)=0$E32

Substituting λ = ω² in equation (32), an equivalent expression is obtained:

$det\left(\mathbf{K}-\lambda \mathbf{M}\right)=0$E33

Equation (33) can be solved for either the squared frequencies of the system λ (forward problem) or the coefficients of K and M matrices (inverse problem). In this section, the case of inverse problem is investigated for the case of an Euler–Bernoulli beam with a single mass attachment det(K).

7.1.5. Inverse problem

The problem of imposing two frequencies on a dynamical system consisting of a beam with a single attached lumped mass is considered here which leads to the following system of equations:

$\left\{\begin{array}{c}\hfill det\left(\mathbf{K}-{\lambda }_a\mathbf{M}\right)=0\hfill \\ \hfill det\left(\mathbf{K}-{\lambda }_b\mathbf{M}\right)=0\hfill \end{array}\right.$E34

where λ_a and λ_b are the desired natural frequencies squared to be imposed on the system (design variables). Due to the fact that no stiffness element is added to the beam, K remains intact and can be determined using equation (28), whereas the mass matrix of the combined system M is given by

$\mathbf{M}={\mathbf{M}}^d+m_1.{\underset{¯}{\varphi }}_1.{{\underset{¯}{\varphi }}_1}^T$E35

where m₁ is the magnitude of the lumped mass and ${\underset{¯}{\varphi }}_1$ is a function of the mass position x₁ and is defined by equation (27). m₁ and x₁ are the unknown variables of the inverse problem. Substituting equation (35) into (34), a system of two equations and two unknowns is obtained whose solution is presented in the next section for the case of a simply supported beam.

In using the determinant method, each desired natural frequency is substituted into a determinant equation of the form of equation (34). Therefore, each desired natural frequency produces a single equation. Therefore, n desired natural frequencies require n “design degrees of freedom” or parameters to be controlled, such as added masses or springs, or their location on the beam. The value of the added mass (or spring) is considered one design degree of freedom and its unknown location is another. The form of the equations produced through imposing a desired frequency on the determinant depends on the choice of eigenfunctions in the assumed modes method and also how the unknown desired parameter (e.g. position of the mass) is included in those eigenfunctions. Here, we have chosen the traditional trigonometric functions as the eigenfunctions because they are the most well-known choices for the Euler–Bernoulli beam. Other choices of eigenfunctions may lead to better computational results with this method.

7.2. Results

The following assumptions are considered in defining the inverse problem:

• The acceptable mass range is a fraction of the mass of the beam, that is, m = cρL for 0 < c < 1.

• The degree of discretization using assumed modes is N=10 for the simply supported beam.

• The acceptable position range is a fraction of the length of the beam L, that is, pL, where 0 ≤ p ≤ 1 .

Maple (Maplesoft) was used to perform the numerical calculations.

The major steps in solving, as well as coding, the inverse eigenvalue problems are outlined here. First, two natural frequencies are chosen as desired input frequencies (λ_a, λ_b). These are the frequencies we seek to impose on the beam with its mass attachment. To insure a solvable problem, we chose known values from previously solved forward problem. The generalized mass and stiffness matrices must be formed. (M^d, K^d) whose diagonal elements are given by

$M_i=1$E36

$K_i=\frac{{\left(i\pi \right)}^{4}EI}{\left(\rho L^4\right)}$E37

The eigenfunction vector, $\underset{¯}{\varphi }$, must also be built using the simply supported beam eigenfunction equation (29). The stiffness matrix is unaffected because no stiffness elements are added. The mass matrix of the combined system is affected by the presence of lumped masses, whereas the stiffness matrix is not. Hence, the matrices are given by

${\mathbf{K}}_t={\mathbf{K}}^d$E38

${\mathbf{M}}_t={\mathbf{M}}^d+c\rho L{\underset{¯}{\varphi }}_{\mathit{pL}}{{\underset{¯}{\varphi }}_{\mathit{pL}}}^T$E39

where c is the mass coefficient and vector ${\underset{¯}{\varphi }}_{\mathit{pL}}$ is defined as

${\underset{¯}{\varphi }}_{\mathit{pL}}=\left[\begin{array}{c}\hfill {\varphi }_1\left(\mathit{pL}\right)\hfill \\ \hfill .\hfill \\ \hfill .\hfill \\ \hfill {\varphi }_N\left(\mathit{pL}\right)\hfill \end{array}\right]$E40

Equations (38) and (39) along with the two desired values of λ are substituted into equation (34). This gives two equations in two unknowns, which shall be solved for the two unknowns, c and p, using fsolve as well as the DirectSearch packages in Maple. Fsolve is Maple’s built-in equation-solving package. The details of the DirectSearch package can be found in the study of Moiseev [17]. The results obtained for c and p include the anticipated results (already known from the forward problem since we chose the desired λ from a known forward problem) plus additional results for c and p. To check whether the order of the frequencies in the frequency spectrum will be conserved or if the results returned by the inverse problem achieve the desired system frequencies, the parameters c and p must be substituted into the forward code in order to obtain the entire spectrum of system frequencies. From symmetry of the boundary condition about the beam midpoint, only half of the beam is considered and the results can be extended to the other half. These results are presented in the left-hand columns of Table 7.

Input given to inverse; determinant method	Solution via inverse; determinant method	Full span of frequency; squared spectrum obtained; through the solution; of forward problem
λ2 = 1210.7430, λ3 = 7703.7093; (Obtained with m = 0.2ρL; l = 0.3L in the forward code)	m = 0.2ρL; l = 0.3L; (This result was obtained; through both Direct Search; and fsolve methods)	76.8682; 1210.7430; 7703.7093; 22705.2481; 49733.6143; 1.19385×10^5; 2.29571×10^5; 3.45310×10^5; 5.95631×10^5; 9.74091×10^5
	m = 0.3895ρL; l = 0.3587L; (This result was obtained through both; Direct Search and fsolve packages)	59.0673; 1210.7430; 7703.7093; 18122.3022; 55694.1292; 1.18846×10^5; 1.91164×10^5; 3.88304×10^5; 5.89699×10^5; 8.69163×10^5

Table 6.

Inverse problem for second and third frequencies for a simply supported beam.

Input given to inverse determinant method	Solution via inverse determinant method	Full span of frequency; squared spectrum obtained; via the solution; of forward problem
λ2 = 1411.9003, λ4 = 20420.5236; (Obtained with m = 0.2ρL; l = 0.4L in the forward code)	m = 0.1516ρL; l = 0.0972L; (This result was obtained through Direct Search and fsolve packages)	94.7851; 1411.9003; 6626.2231; 20420.5236; 51205.1150; 1.11210×10^5; 2.15698×10^5; 3.82314×10^5; 6.29638×10^5; 9.73288×10^5
	m = 0.4695ρL; l = 0.4268L; (This result was obtained through Direct Search package only)	51.2467; 1411.9003; 6082.1204; 20420.5235; 57662.4808; 1.00203×10^5; 2.33797×10^5; 3.33408×10^5; 6.17507×10^5; 9.04206×10^5
	m = 0.2ρL; l = 0.4L; (This result was obtained via Direct Search package only)	71.3527; 1411.9003; 7164.6511; 20420.5236; 60880.6819; 1.05946×10^5; 2.21752×10^5; 3.77388×10^5; 5.73630×10^5; 9.74091×10^5

Table 7.

Inverse problem for second and fourth frequencies for a simply supported beam.

The squared frequencies in the left-hand columns of these tables are chosen from the previously solved forward problem, and their subscripts indicate their order in the frequency spectrum. The middle columns of these tables contain the values of the masses as well as their positions on the beam obtained after substituting the squared frequencies of the first columns into the inverse problem code. Finally, the right-hand- columns are the full span of frequency spectrum obtained after substituting the masses as well as their positions of the middle column into the forward problem. The desired squared frequencies are in bold face in the right-hand column vectors to make it easier for the reader to compare them with the desired frequencies of the left-hand column.

From Table 7, some observations can be made. First, the desired value and order of the two input frequencies were preserved by the inverse problem. Specifically, in Table 7, the two input frequencies remain as the second and fourth system frequencies when the complete frequency spectrum was found for all three possible solutions returned by the inverse problem. Second, for each mass obtained, there must be two corresponding positions that are symmetrical with respect to the middle of the beam. This follows from the symmetry of the boundary conditions of the simply supported beam. Finally, a good equation solver is imperative for this method to work properly. For example, the use of an alternative equation solver (DirectSearch) yielded additional parameters c and p that were not returned by the built-in equation solver of Maple.

In order to verify that the method works when the chosen natural frequencies are not close together, a second example with the simply supported beam is presented. In this second example, the chosen frequencies are the second and eighth natural frequencies. As shown in Table 8, the method worked well in this case. Interestingly, a larger number of possible solutions were obtained in this case, in comparison with the cases shown in Tables 6 and 7. It is observed that the larger the separation in chosen order between λ_i and λ_j, the greater the number of possible solutions to the inverse problem.

Input given to the; inverse problem	Solution via inverse; determinant method	Full span of frequency; squared spectrum obtained; via the solution; of forward problem
λ2 = 1411.900275, λ8 = 3.773881670 × 10^5; (Obtained with m = 0.2ρL; l = 0.4L in the forward code)	m = 0.2ρL; l = 0.6L	71.3527; 1411.9003; 7164.65111; 20420.52363; 60880.6819; 1.05946×10^5; 2.21752×10^5; 3.77388×10^5; 5.73630×10^5; 9.74091×10^5
	m = 0.0807ρL l = 0.6615L	86.6441 1411.9003; 7887.5523; 22430.0485; 56495.1732; 1.26101×10^5; 2.13015×10^5; 3.77388×10^5; 6.37810×10^5; 9.05601×10^5
	m = 0.0752ρL; l = 0.8447L	94.2653; 1411.9003; 6953.8085; 22707.0792; 58428.8453; 1.25669×10^5; 2.31928×10^5; 3.77388×10^5; 5.86872×10^5; 9.10470×10^5
	m = 0.0585ρL; l = 0.7133L	90.83562; 1411.9003; 7750.6323; 24432.2348; 55584.5685; 1.20289×10^5; 2.33870×10^5; 3.77388×10^5; 5.98869×10^5; 9.64184×10^5
	m = 0.0573ρL; l = 0.7944L	93.4974; 1411.9003; 7262.2914; 24320.8984; 60836.8705; 1.20776×10^5; 2.15464×10^5; 3.77388×10^5; 6.30690×10^5; 9.72121×10^5

Table 8.

Inverse problem for second and eighth frequencies for a simply supported beam.

7.3. Conclusion

The determinant method was used to demonstrate how to impose two natural frequencies on a dynamical system consisting of a beam to which a single lumped mass is attached. In this method, the known (design) variables are the two natural frequencies and the unknown variables are the magnitude of the attached mass as well as its position along the beam. The proposed method is easy to code and can accommodate any kind of eigenfunction. It was shown that the expected values of the added mass and its position are recovered from the inverse problem, in addition to the unexpected possible solutions. Even for this simple problem, a unique solution does not exist. The non-uniqueness can be considered as a benefit from the point of view of design possibilities. An investigation of the inverse problem also reveals that the order of the frequencies in the hierarchy of the whole frequency spectrum is conserved. Ideally, the degree of discretization should be the same for both forward and inverse problem. Although not discussed here, a promising observation made by Mir Hosseini and Baddour [21] through numerical simulations was that using a lower degree of discretization for the inverse problem could still give acceptable results from the point of view of engineering design, especially when lower orders of frequencies are involved.

8. Chapter summary and conclusions

In this chapter, two direct methods of solving inverse eigenvalue problems were presented and examples of the methods given. Both methods are not restricted in the type of matrices with which they can be employed. In both cases, nth order polynomial equations are obtained when working with nth order matrices, which must be solved for the inverse problem. Hence, good polynomial solvers, capable of producing the multiple solutions expected of both methods, are required. In using these inverse eigenvalue methods, careful care must be taken in setting up the forward modelling problem because the choice of modelling approach and chosen functions have a large effect on the order and solvability of the ensuing inverse problem.