## Abstract

This paper investigates in-line spring–mass systems (*A*_{n}), fixed at one end and free at the other, with *n*-degrees of freedom (d.f.). The objective is to find feasible in-line systems (*B*_{n}) that are isospectral to a given system. The spring–mass systems, *A*_{n} and *B*_{n}, are represented by Jacobi matrices. An error function is developed with the help of the Jacobi matrices *A*_{n} and *B*_{n}. The problem of finding the isospectral systems is posed as an optimization problem with the aim of minimizing the error function. The approach for creating isospectral systems uses the fact that the trace of two isospectral Jacobi matrices *A*_{n} and *B*_{n} should be identical. A modification is made to the diagonal elements of the given Jacobi matrix (*A*_{n}), to create the isospectral systems. The optimization problem is solved using the firefly algorithm augmented by a local search procedure. Numerical results are obtained and resulting isospectral systems are shown for 4 d.f. and 10 d.f. systems.

## 1. Introduction

Vibrating systems, which have identical natural frequencies are known as isospectral systems. In this work, we investigate chain-like spring–mass systems as shown in figure 1. Such an in-line system may be considered as a real system made up of rigid masses and massless springs or as the discrete approximation of a continuous system. For example, a fixed-free chain-like spring–mass system may be thought of as a discrete finite degree of freedom (d.f.) approximation of a cantilever beam. Thus, isospectral spring–mass systems can yield insights into exploring many other isospectral vibratory systems. Isospectral systems are important in mechanics as they yield alternative usable designs.

The existence of isospectral systems shows that a system cannot be uniquely identified from its spectrum. This proves that for the given system, a unique solution to the inverse problem is not possible. Thus, efforts of designers to create the system from the spectra are not likely to succeed. Isospectral systems are important for the area of damage detection. In structural dynamic systems, damage causes a change in the stiffness or mass properties. The existence of isospectral systems shows that unique identification of damage is difficult from the given spectra. Isospectral systems also allow the designer to create structural dynamic systems with desirable mass and stiffness distributions, which have similar dynamic properties in terms of the spectra.

Pioneering research on *isospectral vibrating systems* was done by Gladwell (1995, 2002), Gladwell & Morassi (1995, 2010) and Gottlieb (1987, 1988). Gladwell (1995) described several ways to form spring–mass systems isospectral to a given system. A pair of isospectral fixed-free systems was obtained by the interchange . This pair was also pointed out by Ram & Elhay (1995). Gladwell (1995) also showed how the indeterminacy associated with the reduction to standard form can be used to get isospectral systems.

The construction of isospectral systems requires an understanding of inverse problems (Boley & Golub 1987; Gladwell & Zhu 1992; Nylen & Uhlig 1997; Chu 1998). Inverse problems (Gladwell 1989, 1999) concern reconstruction of a system from frequency response data and are much more complicated than the corresponding direct problem. There are several applications of inverse problems over wide areas in science and engineering including vibration theory.

The reconstruction of a spring–mass system can be performed by using the interlaced spectrum corresponding to a modified system with different end conditions (Gladwell 1989; Egana & Soto 2000). This requires the knowledge of the eigenvalues of the original system and the eigenvalues corresponding to the modified system whose last mass is fixed. It is known that a simple discrete spring–mass system can be reconstructed uniquely from the poles and zeros of the frequency response function corresponding to sinusoidal forcing at the end. However, Gladwell & Willms (1988) showed that the reconstruction is unique even if forcing is applied at an interior point (not a node of any eigenmode).

A significant step in system construction involves the reconstruction of a symmetric tridiagonal matrix, i.e. a Jacobi matrix. The reconstruction of the Jacobi matrix from spectral data needs numerical techniques or some standard tridiagonalization methods for symmetric matrices. Boor & Golub (1978) used the algorithm devised by Forsythe (1957) to reconstruct a scalar Jacobi matrix. Gladwell & Willms (1988) and Soto (1989) used the *Lanczos algorithm* and *Householder transformation*, respectively, for tridiagonalization. Egana & Soto (2000) used a modified version of the *fast orthogonal reduction algorithm* to reduce computational cost. Once the Jacobi matrix has been found, the remaining part of system construction is to extract the masses and stiffnesses from the Jacobian. The reconstruction of mass–spring systems from the Jacobian matrix can be advantageous (Nylen & Uhlig 1997). A single result for the Jacobian matrix leads to different mass-spring reconstruction results. However, all the reconstructed spring–mass system may not be unique. Hence, some scaling factors need to be introduced to make the solution unique. For example, Gladwell (1995) used a scaling factor .

Despite considerable research, the exploration of isospectral vibrating systems remains an open challenge in the research world. Researchers were able to analytically obtain some isospectral systems. However, it may be possible to obtain other isospectral systems using a numerical approach. In other words, the problem can be framed as an optimization problem where the objective function may have multiple optima. Numerical approaches can complement linear algebraic approaches.

In this work, we pose the problem of finding isospectral mass–spring inline systems as an optimization problem where the objective function is a non-negative error function. The problem becomes similar to finding multiple solutions of a least-square problem. Previously, researchers (Chu 1995; Chen & Chu 1996) had proposed several numerical techniques on the least-squares solution of inverse eigenvalue problems. An iterative method, called *lift and projection* (LP), was proposed by Chu (1995) to construct a Hermitian matrix from its diagonal elements and eigenvalues. Chu (1995) also introduced the projected gradient method, a continuous approach to solve the inverse eigenvalue problem. Later, Chen & Chu (1996) proposed an efficient hybrid method called the LP-Newton, by first applying the LP method to attain a low order of accuracy and then switching to a comparatively faster locally convergent method (Newton method) to reach a high order of accuracy.

Recently, stochastic optimization algorithms (Fouskakis & Draper 2002; Ji & Klinowski 2006) have become very popular because of their advanced search mechanisms. These algorithms search from a population of solution points instead of from a single solution point. They need fitness function information rather than gradient-based knowledge and use probabilistic transition rules in their search for the optimal solution. In this paper, we use the firefly algorithm (FA) augmented by local search (LS) to minimize the objective function. FA (Yang 2008, 2009, 2010; Lukasik & Zak 2009) is a nature-inspired computational algorithm motivated by the sociobiology of fireflies. It has been shown by Yang (2009) that FA can capture multiple optimum points of a multi-modal function.

In this paper, we numerically find spring–mass systems that are isospectral to a given system. Organization of this paper is as follows: In §2, the mathematical model of the system is presented. In §3, optimization criterion and error function are formulated and an example of a 2 d.f. isospectral system is shown to illustrate the formulation. In §4, the firefly and LS algorithms are described. In §5, we present numerical results and resulting isospectral systems for 4 d.f. and 10 d.f. cases. Section 5 presents the conclusion.

## 2. Mathematical modelling of spring–mass system

The governing equation for an undamped spring–mass system with *n*-d.f. is given as
2.1
where the matrices *M*_{n} and *K*_{n} are the mass matrix and the stiffness matrix of order *n*, respectively. In the above equation, *λ*_{i} are the eigenvalues and are the eigenvectors ∀*i*=1,2,…,*n*. Eigenvalues are related to the frequency spectrum and eigenvectors represent the vibrating modes.

Consider a physical spring–mass system consisting of *n* masses *m*_{i}>0 connected by *n* springs of stiffness *k*_{i}>0 with fixed-free condition (*k*_{n+1}=0), as shown in figure 1. The mass matrix *M*_{n} and the stiffness matrix *K*_{n} of the system are given by
2.2
and
2.3
Here, *M*_{n} and *K*_{n} are symmetric-positive definite matrices. Moreover, *M*_{n} (mass matrix) is a diagonal matrix and *K*_{n} (stiffness matrix) is a tridiagonal matrix. The positive definite mass matrix *M*_{n} can be written as . From this factorization, the generalized eigenvalue equation (2.1) can be reduced to the standard form as follows:
2.4
where and . The eigenvalues *λ*_{i} in equation (2.4) are related to the natural frequencies . The matrix *A*_{n} is symmetric tridiagonal positive definite with the same eigenvalues of the system (*K*_{n},*M*_{n}) in equation (2.1). *A*_{n} is called the Jacobi matrix. *A*_{n} has the tridiagonal form as
2.5
where
2.6
The eigenvalues of the given matrix *A*_{n} are *λ*_{i}∀ *i*=1,2,…,*n*. Now, let us consider another Jacobi matrix *B*_{n} of order *n* having the same tridiagonal structure given as
2.7
where *c*_{i}>0 (*i*=1,2,…,*n*) and *d*_{i}>0 (*i*=1,2,…,*n*−1). Let the eigenvalues of *B*_{n} be *μ*_{i}∀ *i*=1,2,…,*n*. The Jacobi matrix *B*_{n} is called isospectral to the given Jacobi matrix *A*_{n}, if both the matrices have the same set of eigenvalues, i.e. *μ*_{i}=*λ*_{i} ∀ *i*=1,2,…,*n*. In this work, we are interested in finding the isospectral matrix (matrices) *B*_{n} for a given *A*_{n}. This problem of finding isospectral *B*_{n} matrix (matrices) is posed as an optimization problem.

## 3. Formulation of the optimization problem

Let *p*_{n} and *q*_{n} be two vectors of size (*n*×1), containing the eigenvalues of the matrices *A*_{n} and *B*_{n}, respectively. The eigenvalues of *A*_{n} and *B*_{n} are arranged in ascending order *λ*_{1}≤*λ*_{2}≤⋯≤*λ*_{n} and *μ*_{1}≤*μ*_{2}≤⋯≤*μ*_{n}, in the vectors *p*_{n} and *q*_{n}, respectively. We need an objective function for the optimization problem. We construct a vector *e*_{n}, which is the difference between *p*_{n} and *q*_{n} vectors. Let (∥⋅∥_{2} denotes 2-norm). We will call *F*_{n} as the *error function*, which is always non-negative for symmetric positive definite Jacobi matrices *A*_{n} and *B*_{n}. The objective function in the optimization problem is the error function (*F*_{n}) and is expressed as
3.1
The system *B*_{n} is called isospectral to the given system *A*_{n} if *μ*_{i}=*λ*_{i} for all *i*. In other words, *B*_{n} is isospectral to *A*_{n} if the error function *F*_{n}=0.

Hence, in order to find the isospectral matrix (matrices) *B*_{n} for a given *A*_{n}, we minimize this error function (objective function). Thus, the optimization problem is given by
3.2
where *C* and *D* are the upper bounds on the diagonal and super-diagonal (or sub-diagonal) elements of the isospectral matrix *B*_{n}. Once we have the values of *c*_{i} and *d*_{i}, we can form the *B*_{n} matrix and obtain its eigenvalues. The values of *C* and *D* are related to the availability of real springs and masses. The number of decision variables in the above optimization problem is *n*+(*n*−1)=(2*n*−1), i.e. *n* diagonal elements and (*n*−1) super-diagonal elements.

In order to reduce the computational effort in the optimization problem, we follow a theoretical result from linear algebra. We know that Trace(*B*_{n})=Trace(*A*_{n}), when matrix *B*_{n} is isospectral to matrix *A*_{n}. We use the above result as a constraint in the optimization problem. Thus the constraint is given by
3.3
We now show a numerical example to illustrate this optimization problem.

### (a) A motivating example

Let us consider a 2 d.f. system having masses *m*_{1}=6.25 kg, *m*_{2}=4 kg and stiffness *k*_{1}=92.5 kN m^{−1}, *k*_{2}=20 kN m^{−1}.

The mass matrix *M*_{2} of the system is given by
3.4
and the stiffness matrix *K*_{2} is given by
3.5

Using the above mass and stiffness matrices, the Jacobi matrix for this system is obtained as
The eigenvalues of the above system (*A*_{2}) are *λ*_{1}=3.8678 and *λ*_{2}=19.1322.

Consider another Jacobi matrix *B*_{2}. The structure of the *B*_{2} matrix is
We assume some random values of the elements of *B*_{2} s.t.
The above Jacobi matrix has eigenvalues *μ*_{1}=6.1690 and *μ*_{2}=17.8310. The error function value becomes *F*_{2}=(*λ*_{1}−*μ*_{1})^{2}+(*λ*_{2}−*μ*_{2})^{2}=6.9886≠0. Hence, the above Jacobi is not isospectral to the given one.

For the matrix *B*_{2} to be isospectral to the given matrix *A*_{2}, the eigenvalues of *B*_{2} must be the same as the eigenvalues of *A*_{2}. The constraint on *B*_{2} matrix is

— Sum of the diagonal elements (trace) of

*B*_{2}should be equal to sum of the diagonal elements (trace) of*A*_{2}; i.e.*c*_{1}+*c*_{2}=23.

Our interest is to find the isospectral *B*_{2} matrix. Let us consider *c*_{1}=*a*_{1}−*ϵ* and *c*_{2}=*a*_{2}+*ϵ*; i.e. 18−*ϵ* and 5+*ϵ*. This satisfies the constraint mentioned above. We need to obtain the value of *d*_{1}. We know that isospectral matrices have the same determinant which gives as follows:
3.6

From the above equation, we obtain 3.7

Now, we assume *ϵ*=2 and obtain *d*_{1}=6.1644, where *c*_{1}=16 and *c*_{2}=7. This corresponds to an isospectral Jacobi matrix given by
where, *m*_{1}′,*m*_{2}′ are the masses and *k*_{1}′,*k*_{2}′ are the stiffness of this isospectral system . We can see that the eigenvalues of are and . Hence, the error (objective) function value becomes . The values of *m*_{1}′,*m*_{2}′,*k*_{1}′ and *k*_{2}′ are obtained as follows. We know that
3.8
3.9
and
3.10
If we choose *m*_{2}′=1, then from equation (3.8), we get *k*_{2}′=7. Using this *m*_{2}′ and *k*_{2}′ in equation (3.9), we obtain *m*_{1}′=1.2895. Now, we know *m*_{1}′ and *k*_{2}′. From equation (3.10), we obtain *k*_{1}′=13.6316. The mass matrix *M*_{2}′ and stiffness matrix *K*_{2}′ for this isospectral system are
3.11
This above system is obtained using *ϵ*=2 and *m*_{2}′=1 kg. In table 1, we present ten 2 d.f. spring–mass systems isospectral to a given system. The given 2 d.f. system has the masses *m*_{1}=6.25 kg, *m*_{2}=4 kg and stiffness *k*_{1}=92.5 kN m^{−1}, *k*_{2}=20 kN m^{−1}. These are obtained for different values of *ϵ* for *m*_{2}′=1 kg. We can also use a different value of *m*_{2}′ and obtain additional isospectral systems.

For this 2 d.f. system, the isospectral system can be described in a pictorial manner as shown in figure 2. In this figure, the *x*-axis represents the value of *c*_{1}. The value of *c*_{2} is obtained as (23−*c*_{1}). The *y*-axis represents the value of . This figure shows the variation of sub-diagonal element *d*_{1} with the diagonal elements (*c*_{1},*c*_{2}) of the Jacobi matrix *B*_{2}. Each point on the curve *PQR* in figure 2 refers to a 2 d.f. spring–mass system, which is isospectral to the given system *A*_{2}. The isospectral system corresponding with *ϵ*=2 and *m*_{2}′=1 is shown by the point *Q* in figure 2. The values of for all values of *ϵ* are shown in this figure.

We also see from this figure 2 that *d*_{1} is real if the values of *c*_{1},*c*_{2} are within the boundary points *P*(3.8678) and *R*(19.1322). The boundary points *P* and *R* are the eigenvalues of the system *A*_{2}. Hence, the value of *ϵ* should be selected in such a way that the diagonal elements are within the eigenvalue range (minimum and maximum eigenvalues). The number of decision variables in this 2 d.f. optimization problem is 3; i.e. two diagonal elements *c*_{1} and *c*_{2}, and one super-diagonal element (*d*_{1}). By selecting the value of *ϵ*, the number of decision variables reduces to one (*d*_{1}). Owing to this fact, we are able to obtain all the isospectral systems using equation (3.7).

For a general *n*-d.f. optimization problem, the number of decision variables is *n*+(*n*−1)=(2*n*−1), i.e. *n* diagonal elements and *n*−1 super-diagonal elements. In the general *n*-d.f. optimization problem also, the selection of *ϵ* will reduce the number of decision variables to *n*−1 (sub-diagonal elements). Hence, we need to use some optimization methods to obtain the isospectral system *B*_{n} for a given system *A*_{n}. The objective function in the optimization problem is the error function given in equation (3.2). The error function corresponds to the square of errors in the eigenvalues of a given Jacobi matrix *A*_{n} and another Jacobi matrix *B*_{n}. It is known that the elements of a Jacobi matrix are related to the mass and stiffness values of a spring–mass system.

For a general *n*-d.f. systems, the objective function is related to the decision variables in a complicated way. Gradient information is also hard to extract from the objective function. Hence, we choose a gradient-free stochastic optimization technique (FA) to solve this optimization problem. The constraint on trace defined in equation (3.3) helps us to reduce the computational time and effort.

## 4. Firefly algorithm

In this section, we briefly discuss the FA developed by Yang (2008, 2009, 2010) and its applications to our problem of obtaining the isospectral system *B*_{n} for a given system *A*_{n}. FA is a population-based stochastic optimization technique. FA is a recent addition to nature-inspired computational methods. It is based on the flashing behaviour of fireflies with the following assumptions (Yang 2008, 2009, 2010).

— All fireflies are unisex so that one firefly will be attracted to other fireflies regardless of its sex.

— Attractiveness of fireflies is proportional to their brightness. Thus for any two flashing fireflies, the less brighter one will move towards the brighter one. Attractiveness and brightness both decrease as the distance between the fireflies increases. If a firefly does not find anyone brighter than itself, then it moves randomly.

— The brightness or light intensity of a firefly is affected or determined by the landscape of the objective function to be optimized.

The variation of brightness (or light intensity) and the formulation of attractiveness are the two important issues in the FA. The brightness (or light intensity) is related to the objective function value. The attractiveness of a firefly is determined by its brightness. Hence, we can see that a point with lower objective function value will be attracted by a point with higher objective function value in a maximization problem. A point is a solution to the optimization problem in a *d*-dimensional solution space (search space), where *d* is the number of decision variables.

Assume there are *N* fireflies in the search space, each associated with a coordinate *x*_{i}(*t*) that represents the point (location) of firefly *i* at iteration *t*. The location is a possible solution to the optimization problem. Each location (solution) *x*_{i},*i*=1,2,…,*N* is a *d*-dimensional vector. The brightness (light intensity) *I* of a firefly at location *x*_{i} can be chosen as
4.1
where *f* is the objective function value at that point *x*_{i}. The attractiveness *β* varies with the distance *r*_{ij} between fireflies *i* and *j* and is given by
4.2
where *β*_{0} is the attractiveness when the distance *r*_{ij}=0. The light intensity decreases with the distance between two fireflies and the light is also absorbed in the media. The absorption coefficient is *γ* in the above equation. The distance between any two fireflies *i* and *j* located at points *x*_{i} and *x*_{j} (at iteration *t*) is the Cartesian distance
4.3
A firefly *i* at location *x*_{i} will be attracted to another firefly *j* owing to the fact that the firefly *j* is brighter than the firefly *i*. The movement of firefly *i* towards firefly *j* is given by
4.4
Equation (4.4) can be written as
4.5
In the above equation, *t* is the iteration number and the second term (*β*(*t*)(*x*_{j}(*t*)−*x*_{i}(*t*))) is due to attraction. The third term (*α*(rand−0.5)) is randomization with control parameter *α*. ‘rand’ is a random number, between 0 to 1, obtained from the uniform distribution. The third term is useful in exploration of search space in an efficient manner. The reason for incorporating –0.5 is to allow the movement within the range of −0.5*α* to +0.5*α*. This value is recommended in earlier FA. A complete description of the FA is given in the literature (Yang 2008, 2009, 2010).

The control parameters of the FA are *α*, *β*_{0} and *γ*. The influence of other solutions is controlled by the value of attractiveness that depends on two parameters: initial attractiveness *β*_{0} and absorption coefficient *γ*. In general, *β*_{0}∈[0,1] has been suggested (Yang 2008, 2009). When *β*_{0}=0, the value of *β*(*t*)=0 and the movement of firefly *i* in equation (4.5) is random and is called non-cooperative search. When *β*_{0}=1, the value of *β*(*t*) in equation (4.5) is a finite value depending on the distance between fireflies *i* and *j*. The movement of firefly *i* depends on the position of firefly *j* and is called cooperative search. The value of *γ* defines the nature of variation of the attractiveness with respect to the distance. Using *γ*=0 corresponds to the constant attractiveness whereas corresponds to the complete random search. Moreover, *γ* is related to the scale of the design variables. The value of *γ* can be derived from
where *γ*_{0}∈[0,1]. The randomness control parameter can be chosen from *α*∈[0,1]. The influence of control parameters on the optimization process can be understood in detail from the previous literature (Yang 2009, 2010; Lukasik & Zak 2009).

### (a) Local search

The LS is a technique used in stochastic optimization algorithms to accelerate convergence to the optimal solution. We introduce the LS procedure described in the literature (Birbil & Fang 2003) in the FA. In LS procedure, information about points in the neighbourhood of the point *x*_{i} is collected. If any of the neighbourhood points are better than the original point *x*_{i}, then this better point is stored instead of the original point *x*_{i}. A better neighbourhood point means that the objective function value of this point is better than the objective function value at point *x*_{i}. Better objective function value implies more (less) in a maximization (minimization) problem. The LS procedure, described by Birbil & Fang (2003) is followed in our study. *LSITER* and *δ* are the two parameters that are used in the LS. The maximum number of iterations that are used in the LS is *LSITER*, and the multiplier for the neighbourhood search is denoted by *δ*.

We know that the solution to the optimization problem is represented as point *x*_{i} in *d*-dimensional space, where *d* is the number of decision variables. In the LS algorithm, the maximum feasible step length *δ* is computed as *Max*_{d}(*u*_{d}−*l*_{d}), where *u* and *l* are upper and lower limits on decision variables. For a given point *x*_{i}, the LS is done by considering one dimension after another. For a given dimension *k*, represents the value of the point *x*_{i} in dimension *k*. This is moved to a new position in dimension *k* and is represented as . The value of is obtained by using two random numbers *ξ*_{1}=rand(0,1), *ξ*_{2}=rand(0,1) and *δ*. If the objective function value of point is better than the objective function value of point *x*_{i}, the point *x*_{i} is replaced by the point . More detail of this LS algorithm is given in the literature (Birbil & Fang 2003). The pseudocode of the FA and the LS procedure are given below:

*Localsearch*(*δ*, *LSITER*)

**for***k*=1:*d*all*d*dimensionsgenerate 1st random number

*ξ*_{1}=rand(0,1)**while**(counter<*LSITER*)generate 2nd random number

*ξ*_{2}=rand(0,1)**if**(*ξ*_{1}>0.5)**else****end if**

**if**()**end if**

increase counter

**end while**

**end for**

*The firefly algorithm*

Initialize the population of

*N*fireflies (solutions)*x*_{i},*i*=1,2,…,*N*Evaluate the objective function value

*f*(*x*_{i}) at each solution*x*_{i}in the populationAssign light intensity

*I*(*x*_{i}) using equation (4.1)Define control parameters

*α*,*β*_{0},*γ***While**(counter<MaxIteration)**for***i*=1:*N*all*N*fireflies**for***j*=1:*N*all*N*fireflies**if**(*I*_{j}<*I*_{i})Calculate the distance

*r*_{ij}using equation (4.3)Vary the attractiveness

*β*(*t*) using equation (4.2)Move firefly

*i*towards*j*using equation (4.5)**end if**Evaluate new solutions and update the light intensity

Find the best solution

Call

*Localsearch*(*δ*,*LSITER*) to update the best solution

**end for***j*

**end for***i*increase counter

**End While**Postprocess results and visualization

## 5. Computational results and discussions

In this section, we present the computational results obtained using the FA with LS procedure. Our objective is to find the isospectral system *B*_{n} for a given system *A*_{n}, where *n* is the number of d.f. of the system. The objective function is the error function given in equation (3.2). We present the results obtained for 4 d.f. and 10 d.f. systems using the FA.

In the FA with LS, initially all the fireflies (*N*) are placed randomly in the search space in a well-dispersed manner. The parameters used in the FA are: population size *N*=20; initial attractiveness *β*_{0}=1.0; light absorption coefficient *γ*=0.04 and randomness control parameter *α*=0.2. In the LS procedure, we have used *δ*=0.02 and *LSITER*=20.

We know that the system *B*_{n} is isospectral to system *A*_{n} when the objective function value given by equation (3.2) is zero. Hence, we stop the algorithm when the objective function value is less than a given tolerance. The tolerance value used in our computation is 10^{−8} for 4 d.f. system and 10^{−6} for 10 d.f. system.

### (a) Computational results for 4 d.f. system

Consider a 4 d.f. mass–spring in-line system i.e. *n*=4. The original Jacobi matrix(*A*_{4}) is of the form:
5.1
Consider another Jacobi matrix(*B*_{4}) of the form
5.2
Our interest is to find the elements of *B*_{4} such that it is isospectral to *A*_{4}. We use the constraint on trace defined in equation (3.3) to reduce the computational time and effort. The constraint for this 4 d.f. system is
5.3
We follow the same methodology used for the 2 d.f. system to modify the diagonal elements of *B*_{4}. The diagonal elements of *B*_{4} can be modified in the following manner by considering two variables *ϵ*_{1} and *ϵ*_{2}:
5.4
5.5
and
5.6
Note that many other ways of modifying the diagonal elements are possible; for example, {*c*_{1},*c*_{2},*c*_{3},*c*_{4}}={*a*_{1}+*ϵ*_{1},*a*_{2}−*ϵ*_{1},*a*_{3}−*ϵ*_{2},*a*_{4}+*ϵ*_{2}}. We have used the above three modifications in our computational algorithm. The key idea is to keep the trace of matrix *B*_{4} equal to the trace of matrix *A*_{4} given in equation (5.3). Thus, for the 4 d.f. system, the optimization problem is expressed as
5.7
Note that *ϵ*_{1},*ϵ*_{2}≠0 ensures a non-trivial solution. The upper bounds(*E*_{1},*E*_{2}) on the variables *ϵ*_{1},*ϵ*_{2} can be determined with the help of Horn's majorization theorem. Applicability of the Horn's majorization theorem to determine the bounds is explained in the appendix.

We now consider a 4 d.f. spring–mass system with mass and stiffness values given as: *m*_{1}=1 kg, *m*_{2}=2 kg, *m*_{3}=3 kg, *m*_{4}=4 kg, and *k*_{1}=50 kN m^{−1}, *k*_{2}=60 kN m^{−1}, *k*_{3}= 70 kN m^{−1}, *k*_{4}=80 kN m^{−1}. The Jacobi matrix *A*_{4} with these values is given by
The eigenvalues of the given system *A*_{4} are 2.2394, 30.6069, 73.8387 and 138.3150. The corresponding natural frequencies are 1.4964, 5.5323, 8.5929 and 11.7607 rad s^{−1}.

For 4 d.f. spring–mass system, the optimization problem given in equation (5.7) has five decision variables: subdiagonal elements *d*_{1},*d*_{2},*d*_{3} and the values of *ϵ*_{1},*ϵ*_{2}. The bounds on the subdiagonal elements are assumed to be as: *b*_{i}−10≤*d*_{i}≤*b*_{i}+10 for *i*=1,2,3. Using Horn's majorization theorem, bounds on the values of *ϵ*_{1},*ϵ*_{2} are: 0<*ϵ*_{1}≤28.3 and 0<*ϵ*_{2}≤17.76, as we use the modifications given by equation (5.4).

The FA is applied to the optimization problem given in equation (5.7) and gives an optimal *B*_{4} for certain values of *ϵ*_{1} and *ϵ*_{2}. For different values of *ϵ*_{1} and *ϵ*_{2}, we can get many *B*_{4} matrices that are isospectral to *A*_{4}. A sample of 10 isospectral spring–mass systems obtained using the FA, for the given system *A*_{4} is presented in table 2. The mass and stiffness values of the isospectral systems are denoted as: *m*_{1}′,*m*_{2}′,*m*_{3}′,*m*_{4}′ and *k*_{1}′,*k*_{2}′,*k*_{3}′,*k*_{4}′ in table 2.

Using the FA, we can obtain the elements of *B*_{4} such that it is isospectral to the given system *A*_{4}. The mass and stiffness values of *B*_{4} are obtained in a manner similar to that explained for the 2 d.f. system in §3. We will explain this for 4 d.f. system. Consider one of the *B*_{4} matrices obtained by the FA given below:
The eigenvalues of *B*_{4} are 2.2394, 30.6069, 73.8386 and 138.3150. Hence, the error (objective) function value becomes *F*_{4}=10^{−8}≃0. The values of *m*_{1}′,*m*_{2}′,*m*_{3}′,*m*_{4}′ and *k*_{1}′,*k*_{2}′,*k*_{3}′ and *k*_{4}′ are obtained in a similar way as described for 2 d.f. system. The steps are given below:
5.8
If we choose *m*_{4}′=1, then from equations (5.8), all the values of *m*_{1}′,*m*_{2}′,*m*_{3}′ and *k*_{1}′,*k*_{2}′,*k*_{3}′ and *k*_{4}′ can be determined. This spring–mass system is shown in table 2 as isospectral system (01).

### (b) Computational results for 10 d.f. system

Now, we apply the same method for a 10-dimensional problem. We consider a 10 d.f. chain-like mass–spring system as shown in table 3. Here, we modify only two diagonal elements of the original Jacobi matrix (*A*_{10}), in one particular manner. We use one *ϵ* and apply to the 6-th and 10-th diagonal terms; i.e. the element *a*_{6} is modified to (*a*_{6}+*ϵ*), and the element *a*_{10} is modified to (*a*_{10}−*ϵ*). The number of decision variables for this 10 d.f. case is 10; subdiagonal elements *d*_{1},…,*d*_{9} and *ϵ*. Bounds on the decision variables are obtained as *b*_{i}−10≤*d*_{i}≤*b*_{i}+10 for *i*=1,2,…,9 and 0<*ϵ*≤10. The results obtained for this 10 d.f. system using the FA are shown in tables 4 and 5.

## 6. Conclusions

The problem of finding isospectral spring–mass systems is posed as an optimization problem with the objective of minimizing an error function. The error function is defined using the Jacobi matrices of *A*_{n} and *B*_{n}. We also use the fact that the trace of two isospectral Jacobi matrices *A*_{n} and *B*_{n} should be identical. We use FA augmented with LS to minimize the error function. The stochastic optimization method is applied to 4 d.f. and 10 d.f. systems. The hybrid algorithm is able to capture multiple optimum points. The isospectral spring–mass systems (scaled) are constructed from the optimum points and all resulting systems are presented. We are able to obtain several isospectral systems using this methodology. Numerical results are presented for 4 d.f. and 10 d.f. systems. The results obtained show that this methodology is useful to obtain isospectral systems.

In this study, we have used an optimization approach to obtain isospectral spring–mass systems. This optimization approach is easy to understand and implement. The FA starts with a population of solutions and converges to an optimal solution. The decision variables in this approach, correspond to the stiffness values (*k*_{i}) and the mass values (*m*_{i}) for *i*=1,2,…,*n*. Once the values of *k*_{i} and *m*_{i} are known, we can compute the eigenvalues and obtain the objective function value. In this approach, we do not need gradient information and spectral information to obtain the optimal solution. Hence, the optimization approach is used.

## Appendix A. Applicability of the Horn's majorization theorem

We illustrate the applicability of Horn's majorization theorem or Schur–Horn theorem (Horn & Johnson 1985) to estimate the bounds on the modification variables.

According to the theorem, the sum of *k* smallest elements of the diagonal vector (vector containing diagonal elements) of a Hermitian matrix should be greater than or equal to the sum of *k* smallest elements of eigenvalue vector (vector containing all the eigenvalues) of the matrix for *k*=1,2,…,*n*−1, where *n* is the order of the matrix.

We explain the applicability of the theorem with the help of an example of 4 d.f. spring–mass system. *A*_{4} is the given Jacobi matrix as shown by equation (5.1). The eigenvalues of *A*_{4} are ordered as *λ*_{1}≤*λ*_{2}≤*λ*_{3}≤*λ*_{4}. *B*_{4} is the modified Jacobi matrix as shown by equation (5.2). Without loss of generality, we assume *c*_{4}<*c*_{3}<*c*_{2}<*c*_{1} for our ease of understanding. Then Horn's majorization theorem can be expressed mathematically as
A1

We consider the modified Jacobi matrix *B*_{4} in a form so that diag(*B*_{4})={*c*_{1},*c*_{2},*c*_{3},*c*_{4}}={*a*_{1}+*ϵ*_{1},*a*_{2}−*ϵ*_{1},*a*_{3}+*ϵ*_{2},*a*_{4}−*ϵ*_{2}}. Using these inequalities (A1), we obtain the upper bounds (*E*_{1},*E*_{2}) on the variables *ϵ*_{1},*ϵ*_{2}
2
and
A3
We also obtain the condition (*a*_{3}+*a*_{4}≥*λ*_{1}+*λ*_{2}) to be satisfied.

- Received February 17, 2011.
- Accepted May 26, 2011.

- This journal is © 2011 The Royal Society