Application of the Characteristic Time Expansion Method for Estimating Nonlinear Restoring Forces

Abstract

This paper proposes a characteristic time expansion method (CTEM) for estimating nonlinear restoring forces. Because noisy data and numerical instability are the main causes of numerical developing problems in an inverse field, a polynomial to identify restoring forces is usually adopted to eliminate these problems. However, results of the way doing are undesirable for a high order of polynomial. To overcome this difficulty, the characteristic length (CL) is introduced into the power series, and a natural regularization technique is applied to ensure numerical stability and determine the existence of a solution. As compared to previous solutions presented in other researches, the proposed method is a desirable and accurate solver for the problem of restoring the force in the inverse vibration problems.

1. Introduction

Identification of nonlinear dynamical system is a kind of inverse problems and is usually encountered in engineering applications. For instance, to specify the parameters of dynamical systems is necessary in optimal processes. It is important to analyze and determine the parameters of the system using experimental testing and numerical methods. However, uses of these methods might lead to some challenging problems in the structural mechanic field because a small measurement error can cause large errors in the results of parameters’ identification.

To overcome these inverse problems, some solvers in the literature were proposed by using numerical techniques and experimental testing. For instance, various publications in [1–10] recommended uses of the damping coefficient, stiffness, and external force for solving the inverse problems. Mode shape, frequency, displacement, and velocity at different times can also be used to estimate these properties successfully [11]. Huang [12] has employed a conjugate gradient method (CGM) to solve the nonlinear inverse vibration problem for the estimation of the time-dependent stiffness coefficient. Recently, Liu [13, 14] has developed a Lie-group shooting method to study the inverse vibration problem for estimating the time-dependent damping and stiffness coefficients and simultaneously derived a closed-form solution to estimate the parameters.

A complete review on the developments of useful methods for the realm of nonlinear system identification can be found in [14]. Reference [15] also proposed the idea of a force state mapping method which is a simple procedure and allows a direct identification of the restoring force for nonlinear mechanical systems. This idea was further extended in [16–18]. Recently, Namdeo and Manohar [19] have modified the force state mapping technique with two alternative schemes of functional representation: (1) reproducing kernel particle method and kriging technique and (2) estimating the parameters of nonlinear system from measured time histories of response under known excitations. Although this numerical method exhibits the capability to reproduce polynomials of specified order and has been applied to mechanical experiments, yet how to ensure numerical stability and avoid noise disturbance are not reported.

The aim of this paper is to develop a simple, multistep regularization algorithm with easy numerical implementation and versatility. A simple polynomial interpolation can be considered as a fit for the time history of displacement response under known excitations. Theoretically, if the order of the polynomial interpolation is large enough, the approximation displacement response can be as accurate as possible. However, in practical applications, it is not a good method to easily increase the order to obtain a highly accurate approximation. In fact, doing so will lead to a highly ill-posed matrix with a high-order function (the Vandermonde matrix), which has been described in [20]. To resolve the ill-posed linear problems for the Vandermonde system, Beckermann [21] and Li [22] claimed that the optimal condition number of the Vandermonde matrix could be expected. Nevertheless, the results of Beckermann [21] and Li [22] show that in the best possible cases, the condition numbers of the Vandermonde matrices still grow exponentially with the order of the interplant polynomial. Because of this, no one is interpolating the data by the high-order polynomials in the usual bases but rather in the Chebyshev polynomials. Hence, how to alleviate the ill-poseness owing to high-order function becomes one of the main tasks in this paper. First of all, we introduce the characteristic length (CL) of computational time into the high-order polynomial expansion to relieve the ill-conditioning of the resulting coefficient matrix of the polynomial expansion and then ensure numerical stability. This concept was first proposed to deal with the Laplace equation using a physical quantity [13, 14, 23, 24]. Recently, the CL has been successfully extended to deal with the Laplace equation and sloshing wave problems [25–27]. Although the CL can enhance the numerical accuracy for solving ill-posed linear matrix, it cannot avoid the effect of measured errors for parameters identification problems. Therefore, how to overcome the instability of the mathematical procedure is quite important. In addition, a small disturbance of measured data has to be considered in the numerical algorithm because they could cause an error identification of the parameter. In order to identify an accurate and stable solution for longer time scales, some special techniques have been used, including of the singular value decomposition (SVD), the SVD with a regularization parameter determined by the L-curve method, and sensitivity analysis. Despite these efforts, the stability problem remains unresolved. To thoroughly overcome these difficulties, this paper further adopts the CL combined with the natural regularization technique [28] to track ill-posed linear problems in numerical procedures. One advantage of this regularization method is that it can determine whether a solution exists for a linear system with the noisy level.

Apart from the current section, Section 2 describes the mathematical formulation of the characteristic time expansion method and introduces the numerical procedure of the matrix CGM. Section 3 gives several numerical examples, including of Duffing’s oscillator, Duffing’s oscillator with negative linear stiffness, Van der Pol’s oscillator, Bouc-Wen class model, and the seat model, to compare results of our method with the analytical solutions. Finally, some concrete conclusions are drawn in Section 4.

2. Basic Formulation

A second-order ordinary differential equation (ODE) for the equation of motion is expressed as 𝑥̈+𝐻(𝑥,𝑥̇)=𝑃(𝑡), (1) where 𝑥 represents the displacement of a system response; 𝑃(𝑡) and 𝐻(𝑥,𝑥̇) are the external excitation and restoring force, respectively. In order to obtain 𝐻, a trivial rearrangement of (1) gives 𝐻(𝑥,𝑥̇)=𝑃(𝑡)−𝑥̈. (2) Here, 𝐻 can be obtained if the quantities, 𝑃(𝑡) and 𝑥̈, on the right-hand side of (2) are known. In general, it is easier to measure the displacement at some discrete sampling times than to directly measure velocities and accelerations. Therefore, if 𝑥₁(𝑡)=𝑔(𝑡) is denoted as the measured displacement, the differentiation of displacements can be expressed as follows: 𝑥₁(𝑡)(𝑡),𝑥̇=𝑔₁(𝑡)=𝑥₂(𝑡),𝑥̇₂(𝑡)=𝑥₃(𝑡). (3)(4)(5) This is however a set of index-three differential algebraic equations (DAEs) [29], which is hard to solve because the amplification of small errors and perturbations in the displacement cause severe numerically ill-posed condition.

2.1. The Characteristic Time Expansion Method

The polynomial interpolation is defined as the interpolation of a given set of data by a polynomial. In other words, given some data points, the aim is to find a polynomial which exactly goes through these points of data.

According to (3), the displacement can be expressed as a polynomial expansion: 𝑥𝑡_𝑘=^𝑛_𝑘=1𝑎_𝑘𝑡_𝑘^𝑛−1,0≤𝑡_𝑘≤𝑡_𝑓, (6) where 𝑡_𝑘 denotes each discrete time, 𝑥(𝑡_𝑘) denotes the displacement at each time, 𝑡_𝑓 denotes the final time, and 𝑎_𝑘 denotes the unknown coefficient. In many engineering applications, one wants to interpolate the data as accurate as possible. But this is limited by the interpolation of 𝑛 data with (𝑛−1)-order polynomials, where the resultant Vandermonde matrices are highly ill-conditioned as measured by the Lebesgue constant 2𝑛/[𝑒(𝑛−1)ln𝑛].

In this study, we introduce the characteristic length (CL) into (6) and express as follows 𝑥𝑡_𝑘=𝑎₀+^𝑚_𝑘=1𝑎_𝑘𝑡_𝑘𝑇₀^𝑘,0≤𝑡_𝑘≤𝑡_𝑓,𝑡_𝑓<𝑇₀, (7) where 𝑇₀ denotes the CL. Differentiation of (7) yields velocity and acceleration and they are expressed as follows: 𝑥̇𝑡_𝑘=^𝑚_𝑘=1𝑘𝑇₀𝑎_𝑘𝑡_𝑘𝑇₀^𝑘−1,𝑥̈𝑡_𝑘=^𝑚_𝑘=1𝑘(𝑘−1)𝑇²₀𝑎_𝑘𝑡_𝑘𝑇₀^𝑘−2. (8) The polynomial expansion in (7) can be used to describe the displacement of a system. Hence, (7) can be expressed as a linear equation system with 𝑛=𝑚+1: 1𝑡₁𝑇₀𝑡₁𝑇₀²⋯𝑡₁𝑇₀^𝑚1𝑡₂𝑇₀𝑡₂𝑇₀²⋯𝑡₂𝑇₀^𝑚1𝑡₃𝑇₀𝑡₃𝑇₀²⋯𝑡₃𝑇₀^𝑚1𝑡⋮⋮⋮⋱⋮_𝑚𝑇₀𝑡_𝑚𝑇₀²⋯𝑡_𝑚𝑇₀^𝑚𝑎₀𝑎₁𝑎₂⋮𝑎_𝑚=𝑥𝑡₁𝑥𝑡₂𝑥𝑡₃⋮𝑥𝑡_𝑚. (9) We denote the above equation by 𝐑𝐜=𝐛₁, (10) where 𝐜=[𝑎₀,𝑎₁,𝑎₂,…,𝑎_𝑚]^𝐓 is the vector of unknown coefficients.

2.2. The Matrix CGM for Ill-Posed Linear System

When a matrix is ill-posed and measured data contains noisy disturbances, it is difficult to ensure the stability of the system using the conventional regularization techniques. Therefore, Liu et al. [28] proposed a natural regularization method, which proves that a solution exists when ill-posed matrix and noisy disturbances occur. This method can be described by the following matrix equation: 𝐑^𝐓𝐔^𝐓=𝐈_𝑚,thatis,(𝐔𝐑)^𝐓=𝐈_𝑚. (11) If 𝐔 is the inversion of 𝐑, numerically, 𝐔 becomes a left-inversion of 𝐑. Then we have 𝐑𝐑^𝐓𝐔^𝐓=𝐑. (12) Let 𝐑𝐗₀=𝐲₀. (13) Given 𝐗₀, say 𝐗₀=𝚰=[𝟏,…,𝟏]^𝐓, 𝐲₀ can be directly obtained because 𝐑 is a given matrix. Hence, we have 𝐲^𝐓₀𝐔^𝐓=𝐗^𝐓₀,thatis,𝐗₀=𝐔𝐲₀. (14) When (11) and (14) are combined together, they create an over-determined system to calculate 𝐔^𝐓. The over-determined system can be written as 𝐁𝐔^𝐓=𝐈_𝑚𝐗^𝐓₀, (15) where 𝐑𝐁∶=^𝐓𝐲^𝐓₀ (16) is an 𝑛×𝑚 matrix with 𝑛=𝑚+1. Multiplying (14) by 𝐁^𝐓 yields an 𝑚×𝑚 matrix equation: 𝐑𝐑^𝐓+𝐲₀𝐲^𝐓₀𝐔^𝐓=𝐑+𝐲₀𝐗^𝐓₀. (17) Besides the primal system shown in (10), we need to solve the dual system with 𝐑^𝐓𝐲=𝐛₁. (18) Applying the operators in (17) to 𝐛₁ and utilizing the above equation, that is, 𝐲=𝐑^𝐓𝐛₁, we can obtain 𝐑𝐑^𝐓+𝐲₀𝐲^𝐓₀𝐲=𝐑𝐛₁+𝐗₀⋅𝐛₁𝐲₀, (19) where 𝐲₀=𝐑𝐗₀.

Replacing the 𝐑 in (19) by 𝐑^𝐓, we have a similar equation for the primal system in (10) 𝐑^𝐓𝐑+𝐲₀𝐲^𝐓₀𝐜=𝐑^𝐓𝐛₁+𝐗₀⋅𝐛₁𝐲₀, (20) where 𝐲₀=𝐑^𝐓𝐗₀.

Finally, when 𝐜 of (20) is calculated by the CGM, the restoring force, velocity and acceleration can be obtained from (2) and (8).

3. Numerical Examples

Example 1. In this case, we consider a Duffing oscillator [29] and a second-order ODE to describe the forced vibration of a nonlinear structure by 𝑥̈+𝛾𝑥̇+𝛽𝑥+𝛼𝑥³=𝑃(𝑡), (21) where the parameters are fixed as 𝛼=1, 𝛽=−1, and 𝛾=0.3. The restoring force can be expressed as follows: 𝐻(𝑥)=𝑥³−𝑥. (22) In order to identify the restoring force 𝐻 as a function of 𝑥, a monotonic function of 𝑡 is required. In this example, 𝑥(𝑡)=𝑡²−8 is used to obtain the external force and is given by 𝑃(𝑡)=𝑡²−8³−𝑡²+0.6𝑡+10. (23) To test the stability of the numerical method, the order of the polynomial and computational time are increased. The restoring force in the initial and final time changed very rapidly. To understand the CL effect, 𝑚=201, 𝐗₀=𝐈, and 𝜀=1×10⁻¹⁶ are fixed. The maximum estimation error of 𝐻 with different CLs, shown in Figure 1, is smaller than 10⁻⁶.

Figure 1 

For Example 1, showing the error of estimation with different CLs.

It can be seen that including the CL into this case is efficient to overcome an ill-posed matrix. Furthermore, by fixing 𝑇₀=1200, the exact solutions for velocity and acceleration can be determined. The numerical results are shown in Figures 2, 3, and 4. According to the numerical results, the maximum estimation errors of 𝐻 are found to be smaller than 6×10⁻¹⁰. We can find that applying the CL and matrix regularization method with the CGM can provide highly stable and accurate solutions. In order to further test the stability of the present method under different noise levels, we also consider 𝑥_𝑖=𝑥_𝑖+𝜎𝑅(𝑖) (24) as an input into the estimation equations, where 𝑅(𝑖) is a random number in [−1,1], and 𝜎 is a noise level. With different noise levels 𝜎=1%, 3%, 5%, and 10%, the computed profile of restoring forces is shown in Figure 5. Figure 5 also shows that the maximum estimated errors of 𝐻 are smaller than 10⁻¹ with noisy disturbances. We can find that the CL can effectively overcome numerical instability under the effect of the high-order function and large noise disturbances. Hence, we can see that the present method has a highly numerical accuracy and stability under the effect of the high-order function and large noise.

(a)

(b)

Figure 2 

For Example 1, (a) comparing the exact solution and numerical one for velocity and (b) displaying the estimated error of velocity.

(a)

(b)

Figure 3 

For Example 1, (a) comparing the exact solution and numerical one for acceleration and (b) displaying the estimated error of acceleration.

(a)

(b)

Figure 4 

For Example 1, (a) comparing the exact solution and numerical one for restoring force and (b) displaying the estimated error of restoring force.

(a)

(b)

Figure 5 

For Example 1, (a) comparing estimated and exact restoring forces under different noise level and (b) displaying the error of estimation with 𝜎=1%, 3%, 5%, and 10%, respectively.

Example 2. In this case, the Van der Pol oscillator is one of the nonlinear benchmark problem, and 𝐻(𝑥,𝑥̇) is given by 𝐻(𝑥,𝑥̇)=𝑥+𝑥²−1𝑥̇. (25) In this equation, 𝑥 is given by 𝑥(𝑡)=𝑡³/3−8, and then, the external force can be obtained as 𝑃(𝑡)=𝑡³3+𝑡−8𝑡³3−8𝑡²−1𝑡²−8+2𝑡. (26) In this calculation, by fixing 𝜀=1×10⁻¹⁶, 𝐗₀=𝐈, and 𝑇_𝑓=4, the numerical accuracy and stability of different parameters can be tested, including 𝑇₀=101, 𝑚=100 and 𝑇₀=2×10⁵, 𝑚=500, respectively. The maximum estimation error of 𝐻 shown in Figures 6 and 7 are smaller than10⁻⁸.

(a)

(b)

Figure 6 

For Example 2, (a) comparing exact solution and numerical one for restoring force with𝑇0=101 and 𝑚=100 under 4 seconds and (b) displaying the estimated error of restoring force.

(a)

(b)

Figure 7 

For Example 2, (a) comparing exact solution and numerical one for restoring force with 𝑇0=2×105 and 𝑚=500 under 4 seconds and (b) displaying the estimated error of restoring force.

To test the numerical stability of increasing the computational time by 10 seconds, the parameters are fixed as 𝑇₀=120, 𝑚=100 and 𝑇₀=4×10⁵, 𝑚=500, respectively. The maximum estimation errors of 𝐻, shown in Figures 8 and 9, are smaller than 10⁻⁸. From the numerical solutions in Figures 6–9, it shows that the present method can keep the same numerical accuracy with the increase of the CL when the computational time increases.

(a)

(b)

Figure 8 

For Example 2, (a) comparing exact solution and numerical one for restoring force with 𝑇0=120 and 𝑚=100 under 10 seconds and (b) displaying the estimated error of restoring force.

(a)

(b)

Figure 9 

For Example 2, (a) comparing exact solution and numerical one for restoring force with 𝑇0=4×105 and 𝑚=500 under 10 seconds and (b) displaying the estimated error of restoring force.

This example demonstrates the results of fixing the parameters 𝑇_𝑓=10, 𝑇₀=1.2×10⁴, and 𝑚=201 under different noise levels with 𝜎=1%, 3%, 5%, and 10%. The computed profile of 𝐻 is plotted in Figure 10. Figure 10(a) compares the restoring force with exact one, and the maximum estimation error of 𝐻 shown in Figure 10(b) is smaller than 10⁻¹. From numerical result, we can find that the maximum error occurs because the displacement is equal to zero. That is, this present method can overcome the effect of the high-order function and large noise simultaneously. Therefore, it is found that the proposed method is accurate especially when noisy disturbances are encountered.

(a)

(b)

Figure 10 

For Example 2, (a) comparing estimated and exact restoring forces and (b) displaying the error of estimation with 𝜎=1%, 3%, 5%, and 10%, respectively.

Example 3. The Bouc-Wen class model is one of the most widely used to efficiently describe smooth hysteretic behavior in engineering application. For a structural element described by the Bouc-Wen classical model, the restoring force is written as 𝐻(𝑥,𝑥̇)=𝛼𝑤²_𝑛𝑥+(1−𝛼)𝑤²_𝑛𝑧, (27) where 𝛼 is a post- and preyield stiffness ratio, 𝑤_𝑛 denotes natural frequency, and 𝑧 is an auxiliary variable that represents inelastic behavior. The evolution of 𝑧 is determined by an auxiliary ordinary differential equation, which can be written in the form of 𝑧̇=𝐴𝑥̇−𝛽𝑥̇|𝑧|^𝑛−𝛾|𝑥̇||𝑧|^𝑛−1𝑧, (28) where 𝑧̇ denotes the derivative of 𝑧 with respect to time and 𝐴 and 𝑛 are parameters that control the scale and sharpness of the hysteresis loops, respectively. Parameters, 𝛽 and 𝛾, are used to control the shape of the hysteresis loop. In order to estimate the velocity and the restoring force of the Bouc-Wen model, the group preserving scheme (GPS), which was proposed by Liu [30], is adopted. The restoring force obtained by the GPS is referred to as the exact restoring force. We consider a system that has the parameter values as 𝐴=0.5, 𝛽=−5.0, 𝛾=5.0, 𝑛=1.4, 𝛼=0.4, 𝑤_𝑛=3.0, 𝑡₀=0.0, 𝑡_𝑓=50.0, Δ𝑡=0.01, and 𝑃(𝑡)=sin(𝑡); the initial condition of (𝑥,𝑥̇,𝑧)̇ is given as (0.0, 0.0, 0.1). The exactly computed profile of𝐻is plotted in Figure 11.

Figure 11 

For Example 3, restoring force surface based on GPS.

To obtain 𝐻 using the present method, the parameters 𝑇₀=8, 𝐗₀=𝐈, 𝑚=201, and 𝜀=1×10⁻¹⁴ are fixed. The computed profile of 𝐻 at 40 to 42 seconds is plotted in Figure 12(a), and the maximum estimated error of 𝐻, shown in Figure 12(b), is smaller than 1.4×10⁻². Further, under a noise of 𝜎=5%, the computed profile of 𝐻 at 40 to 42 seconds is plotted in Figure 13(a). The maximum estimated error of 𝐻, shown in Figure 13(b), is smaller than 2.5×10⁻². We see from Figures 12 and 13 that the maximum estimated error of 𝐻 still keep in the order of 10⁻² under a noise of𝜎=5%. That is, we can use the present method to achieve a more accurate and stable solution under a large noisy level.

(a)

(b)

Figure 12 

For Example 3, (a) comparing estimated and exact restoring forces and (b) displaying the estimation error of restoring force.

(a)

(b)

Figure 13 

For Example 3, (a) comparing estimated and exact restoring forces with 𝜎=5% and (b) displaying the estimation error of restoring force.

Example 4. As in Example 3, we consider the viscosity damping effect into the Bouc-Wen classical model and estimate the restoring force described as 𝐻(𝑥,𝑥̇)=2𝜉𝑤_𝑛𝑥̇+𝛼𝑤²_𝑛𝑥+(1−𝛼)𝑤²_𝑛𝑧, (29) where 𝜉 is the viscosity damping ratio. In this example, the parameter values are fixed as 𝐴=0.8, 𝛽=4.0, 𝛾=2.1, 𝑛=1.4, 𝛼=0.4, 𝑤_𝑛=3.0, 𝜉=0.15, 𝑡₀=0.0, 𝑡_𝑓=10, Δ𝑡=0.01, 𝐗₀=𝐈, and 𝑃(𝑡)=0.4sin(𝑡); the initial condition of (𝑥,𝑥̇,𝑧)̇ is given as (0,0,0.1). The exact computed profile of 𝐻 is plotted in Figure 14.

Figure 14 

For Example 4, restoring force surface based on the GPS.

In this case, we use the same solver with the same parameters of Example 3. The computed profile of 𝐻 at 5 to 7 seconds is plotted in Figure 15(a). The maximum estimated error of 𝐻, shown in Figure 15(b), is smaller than 2.5×10⁻³. Again, with a noise of 𝜎=5%, the computed profile of 𝐻 at 5 to 7 seconds is plotted in Figure 16(a). The maximum estimated error of 𝐻, shown in Figure 16(b), is smaller than 1×10⁻². It can be seen in Figures 15 and 16 that the maximum errors are smaller than 10⁻² under a noise of 𝜎=5%. Therefore, we conclude that for the smooth hysteretic behavior by the Bouc-Wen classical model, the accurate and stable solutions in Examples 3 and 4 are available when the proposed method is adopted.

(a)

(b)

Figure 15 

For Example 4, (a) comparing estimated and exact restoring forces and (b) displaying the estimation error of restoring force.

(a)

(b)

Figure 16 

For Example 4, (a) comparing estimated and exact restoring forces with 𝜎=5% and (b) displaying the estimation error of restoring force.

Example 5. In order to test the numerical stability of the CTEM used for the restoring force problem of discontinuous type, the vehicle seat problem is considered. When the seated human body is exposed to vertical vibration, the single-degree of freedom model can be used to describe its seat-person mathematical behavior, as shown in Figure 17, and is given by [31] 𝑀𝑥̈+𝑐₁𝑥̇+𝑐₂|𝑥̇|̇𝑘𝑥+¹1+𝑘₂|𝑥|𝑥=𝑃(𝑡). (30) The parameters of the discontinuous typed vehicle seat model are given as 𝑘₁=48000Nm⁻¹, 𝑘₂=24000Nm⁻¹, 𝑐₁=300Nsm⁻¹, 𝑐₂=1500Nsm⁻¹, 𝑀=8kg, and𝑀₁=42kg. Here the external force is given by 𝑃(𝑡)=0.04cos(𝑡), and the parameters are given by 𝑇₀=11, 𝑡_𝑓=10, 𝐗₀=𝟎.𝟎𝟎𝟏, and 𝜀=1×10⁻¹⁴, respectively. The computed profile of 𝐻 by 𝑚=51 and 201 is shown in Figure 18(a). The maximum estimated error of 𝐻, shown in Figures 18(b) and 18(c), is smaller than 5×10⁻³. Numerical results show that this present method does not exhibit the numerical oscillation (Gibb’s phenomenon) when the high-order function is used. Hence, this example demonstrates that the present method has a high level of accuracy and stability for the restoring force problem of discontinuous type.

Figure 17 

For Example 5, seat-person model of single-degree of freedom system [31].

(a)

(b)

(c)

Figure 18 

For Example 5, (a) comparing estimated and exact restoring forces, (b) displaying the estimation error of restoring force with 𝑚=51, and (c) displaying the estimation error of restoring force with 𝑚=201.

4. Conclusions

In nonlinear mechanical system analysis, the inverse vibration problem is difficult to solve under the measured data with noise. This paper has successfully combined the CTEM with a natural regularization algorithm to determine the unknown restoring force. Due to inclusion of the CL to retain high accuracy and stability, the approximation method can avoid the numerical instability caused by a high-order polynomial function. In addition, when the measured data is contaminated by a large noise, the errors can be controlled by utilizing a natural regularization technique and increasing the CL. In summary, the presented method is an effective and convenient approach to solve the inverse vibration problems.