Power Secant Method applied to natural frequency extraction of Timoshenko beam structures

This work deals with an improved plane frame formulation whose exact dynamic stiﬀness matrix (DSM) presents, uniquely, null determinant for the natural frequencies. In comparison with the classical DSM, the formulation herein presented has some major advantages: local mode shapes are preserved in the formulation so that, for any positive frequency, the DSM will never be ill-conditioned; in the absence of poles, it is possible to employ the secant method in order to have a more computationally eﬃcient eigenvalue extraction procedure. Applying the procedure to the more general case of Timo-shenko beams, we introduce a new technique, named “power deﬂation”, that makes the secant method suitable for the transcendental nonlinear eigenvalue problems based on the improved DSM. In order to avoid overﬂow occurrences that can hinder the se-cant method iterations, limiting frequencies are formulated, with scaling also applied to the eigenvalue problem. Comparisons with results available in the literature demonstrate the strength of the proposed method. Computational eﬃciency is compared with solutions obtained both by FEM and by the Wittrick-Williams algorithm.


INTRODUCTION
Using an improved DSM, exact mode shape and natural frequency calculation of Euler-Bernoulli skeleton structures is discussed at length by Dias and Alves [4].According to the Euler-Bernoulli theory, the beam transverse vibrations are only due to bending deflection.Hence, shear deflection is not taken into account, which restricts the work in [4] to the slender members and low order modes.
For short beams, the Euler-Bernoulli beam theory is no longer appropriate, and knowledge gained from the Timoshenko beam theory should be used [5,[7][8][9]12].This theory considers concomitantly the bending and shear deflections as well as the mass and the rotatory inertia per unit length.
Timoshenko beam theory will not be developed in detail here.For the purposes of the present work, we will use the theoretical results presented in Ref. [5].A complete basic text on the theory and some of its applications can be found in Ref. [8].Additionally, greater elucidation about the theory can be gathered from [7] and [9].
We stress that in the vibration analysis, the concept of slenderness is associated with the distance between null bending moments.Hence, even for slender beams, when the modal order increases, the shear deflection and the rotatory inertia will have an increasing influence on the natural frequencies.We may conclude, therefore, that using the Timoshenko theory is always more suitable when there are short members in the model and/or there is an interest in calculating modes of a higher order.Nevertheless, as shown, the price paid for using this theory reflects the mathematical complexity, which is greater than for the simpler Euler-Bernoulli theory.
Additionally, as the very title of Ref. [5] suggests, it can be interesting to compute the effect of axial load in flexural vibrations.When a beam suffers transverse displacement, its cross-sections suffer translations and rotations in such a manner that the axial loads perform work.If such work is added to the total amount of the elastic energy, the natural frequencies are affected.Therefore, tensile axial loads increase the natural frequencies, while, inversely, compression axial loads decrease the natural frequencies.In this last case, depending on the magnitude of the load, the "dynamic buckling" phenomenon can occur, which is characterized by the total nullification of the bending rigidity of the beam.
The fundamentals of the dynamic problem herein investigated are given in sections 2 and 3, which lead to the DSM formulation given in section 4. Section 5 investigates the numerical approach necessary for the implementation of the secant method [1].In section 6, this method is modified to accommodate special characteristics of the non-linear transcendental eigenvalue problem explored here.
In this context, two particular numerical aspects are discussed: overflow occurrences and deflation of the DSM determinant.
Because the secant method is an iterative process, its implementation is limited by the computer numerical structure such that overflow occurrences should be predicted and circumvented; without such planning, this process will certainly fail.Therefore, for the calculation of the hyperbolic function arguments as well as for the DSM determinant calculation, the numerical techniques described in section 5.3 should be considered.
For linear eigenvalue problems, whose similarity has been established by the finite elements method [1], the dynamic stiffness matrix determinant is a polynomial of equal or lesser order than the number of degrees of freedom.In this case, polynomial deflation [1,13] works perfectly as an auxiliary technique of the secant method, allowing for iterative calculation of a certain eigenvalue, in a recursive process, if the lower order eigenvalues are known.
For the present DSM case, nevertheless, its determinant exhibits a non-polynomial form so that the polynomial deflation cannot be successfully utilized.Thus, we introduce in the section 6 the technique called "power deflation": in the calculation of the nth-order eigenvalue, the determinant is deflated by a product of monomial powers; each one of these containing an eigenvalue whose order is lower than n.

MEMBER EQUILIBRIUM EQUATIONS
For an elementary volume of length dx (Fig. 1), the equilibrium of lateral forces can be written as: where F C (x, t) is the shear force due to bending, m = ρA is the uniform distribution mass and v(x, t) is the transverse displacement.
Analogously, the equilibrium of the elementary volume moments gives: where M F (x, t) is the bending moment, j = ρI is the uniformly distributed rotatory inertia, θ(x, t) is the angle of rotation of the cross-section, and P is a static axial load.
The static axial load is supposedly time-invariant, constant throughout the length and positive under tensile condition.
On the other hand, the shearing condition of the elementary volume can be written as: where G = E/[2(1+υ)] is the shear modulus and A S is the shear area that is given as a fraction of the total area of the cross-section (A).The shear area characterizes the stiffness in relation to the deflection due to shear, so that if A S = ∞ there is no such deflection.Thus, eliminating F C (x, t) and M F (x, t) from these equations, it results in the following differential equation, applicable for both transverse displacement and rotation [5]: where the dot denotes v(x, t) or θ(x, t).Noted in this equation is the existence of crossed terms that are involved in the relations between the beam parameters: EI = flexural rigidity, GA S = shear rigidity, m = mass per unit length, j = rotatory inertia and P = static axial load.Additionally, using the same elementary volume of length dx (Fig. 1), the dynamic equilibrium of axial forces, from which P is excluded because it is constant, can be written as: where F N (x, t) is the normal force and u(x, t) is the axial displacement.

MEMBER EQUILIBRIUM SOLUTION
In free vibration, the displacements are synchronous and harmonic and can be written as: In Plane Rotation (6) stressing that, for the present theory, the rotation θ(x, t) cannot be directly obtained from the derivative of v(x, t), as it could in the Euler-Bernoulli theory.Substitution of expressions (6) for displacements in the equations from the previous section, after some algebraic transformations, leads to: where Ā, B, C, D, Ē and F are integration constants to be determined by using appropriate boundary conditions, θ T and θ H are parameters to be given forward and: Analogously, in free vibration, the internal forces can be written as: Latin American Journal of Solids and Structures 7(2010) 307 -333 Bending Moments (9) with where the parameters F A , F T , F H , M T and M H are given forward.
For convenience, we define the following non-dimensional parameters: where, in Eq. ( 11) we have: Thus, we finally write: where Latin American Journal of Solids and Structures 7(2010) 307 -333 with the non-dimensional parameters α defined by:

Element local dynamic stiffness matrix
Employing the convention for the end displacement and the end forces shown in Fig. 1, as well as the origin there indicated for the local reference system, and defining for J = A, T or H : the dynamic stiffness matrix of an axially loaded Timoshenko member is given by: where the matrices φ φ φ Q and φ φ φ P are, respectively: End displacements versus integration constants relationship Latin American Journal of Solids and Structures 7(2010) 307 -333 End forces versus integration constants relationship In the expressions (28) and (29), the parameters θ T , θ H , F A , F T , F H , M T and M H depend on the frequency ω and are defined by Eqs.(18).
By the usual FEM procedures, the matrix of expression (27), after applied to each member of the model, is useful to assemble the classical global DSM [5].However, for frequencies that coincide with the natural frequencies of the members in the clamped-clamped condition, the classical DSM has poles, not roots.
This can be easily seen by observing that the matrix ϕ Q (ω) has no inverse when ω is one of the natural frequencies of the member in the clamped-clamped condition so that, as consequence, the DSM determinant grows to infinity.
For this reason, the classical DSM is not suitable in an iterative method that needs determinant values to estimate the natural frequencies.To overcome this difficulty, in the following section it is shown how the matrices ϕ P (ω) and ϕ Q (ω) can be used in order to assemble an improved version of the global DSM that has no poles [4].Then, the improved DSM can be used to extract eigenvalues by the secant method [1], as it will be presented in section 6.

Improved global dynamic stiffness matrix
By applying the standard FEM procedures, the process by which the matrices φ φ φ P (ω) and φ φ φ Q (ω) are employed for the in-plane model assembly, in the light of the Timoshenko theory, follows the same steps that are described in Ref. [4] for the Euler-Bernoulli theory.For the sake of brevity, the details will not be presented here.The terms are: M Total number of members of the model; N Total number of nodal degrees of freedom; Q o Vector of the amplitudes of the nodal displacements; Θ Vector that contains the constants of integration of all the members; Ω Matrix that unites the nodal concentrated masses and springs; ψ Q Matrix that determines the connection between the members of the model; Φ P Matrix that unites the matrices φ φ φ P (ω) of all the members; Φ Q Matrix that unites the matrices φ φ φ Q (ω) of all the members.The eigenvalue problem, which is established through the exact improved dynamic stiffness matrix Ψ(ω), results in the following form [4]: being necessary to search for ω such that: For this eigenvalue problem, which does not show poles in the natural frequencies [4], it is possible to obtain the roots of Eq. ( 31) and then return to Eq. ( 30) to obtain the corresponding vectors Θ and Q o .The first of these vectors contains the integration constants, which in the latter instance determine the amplitudes of the internal displacements, and the second vector determines the amplitudes of the nodal displacements.
In these terms, for the calculation of each natural mode shape, the solution of Eq. ( 30) is based in a method that linearizes the eigenvalue problem in the vicinity of the corresponding natural frequency, as described in Ref. [4].For the sake of brevity, this method will not be reproduced here.
In Ref. [4], a simple method of determinant tracing was presented to solve Eq. ( 31).Although robust, it is not computationally efficient because the method requires successive evaluations of the determinant, at frequency intervals not greater than the tolerance specified for the natural frequencies calculations.Thus, the objective is to implement the secant method with the use of power deflation, which is computationally much more efficient.Having in mind this implementation, in the following section some particular issues about the present eigenvalue problem are discussed.
Based on transcendental member stiffness matrices, Ref. [9] defines a new dynamic stiffness matrix, called K ∞ , which also eliminates the poles from the eigenvalue problem.However, this method requires further developments in order to include rigid off-set and end release in its formulation as well as has no direct capability to determine the local mode shapes for which Q o = 0 and Θ ≠ 0.

PARTICULAR ISSUES
As we have shown, derivation of the Timoshenko theory is much more laborious than the Euler-Bernoulli theory; consequently, new physical and mathematical questions arise, which we will discuss in the following sections.

Limiting compression axial load
By imposing that the denominator in expressions (21) and ( 22), for the parameters α T and α H , respectively, should be real and finite, we may write: Latin American Journal of Solids and Structures 7(2010) 307 -333 which limits the axial load to: Thus, the lower the shear rigidity, GA S , the lower the axial load that can be applied.In practical terms, for actual materials such a limit is sufficiently high to make feasible the calculation for any realistic structure.
On the other hand, even for infinite shear rigidity, the compressive axial load can be physically limited by the buckling phenomenon.As an illustrative example, we may consider the case of a simply supported beam with no shear deflection (A S = ∞ ⇒ s = 0).For this situation of a Rayleigh beam, the equations in section 3 are applied in the appendix A and give for the natural frequencies: where, in accordance with the Euler buckling load : Hence, if the load is compressive (P < 0) and has magnitude greater or equal to the buckling load, null or complex natural frequencies will occur; a fact that characterizes the physical impossibility of dynamic equilibrium.

Transition frequency
The parameter α T , present in the arguments of the trigonometric functions, will always result in a real number because the conditions ∂ 2 < 0 and ∆ < 0 cannot occur concomitantly.
On the other hand, the parameter α H , in the arguments of the hyperbolic functions, can result in a real number or a purely imaginary number, without thereby implying any impossibility, either mathematical or physical.
In these terms, a particular frequency value exists for which the transition of α H occurs from a real number to an imaginary number.Thus, the use of expression (25) and the condition (32) imposes that: to obtain, from Eq. ( 11), the following expression for the transition frequency: The arguments of the hyperbolic functions are purely imaginary numbers for frequencies greater than the transition frequency, a fact that indicates that such functions convert to the trigonometric functions.From a computer calculation point of view, this is a favorable condition because the trigonometric functions, by being limited, do not present overflow problems.
As noted in Eq. (36), if the shear deflection is neglected (s = 0) there is no transition and overflow can occur because the hyperbolic function arguments are always real numbers.The same fact occurs when the rotatory inertia is neglected (r = 0).On the other hand, it can be noted that the transition frequency does not depend on the member length and therefore it is not affected by mesh refinement.

Treatment of overflow conditions
For computer calculation, there is always an upper limit for the largest number that can be stored in the RAM (R max ).If during some step of an iterative calculation such a limit is surpassed, this calculation is certainly compromised due to overflow.
The equations in sections 3 and 4 have parameters that increase with the frequency ω such that this may lead to overflow error, particularly because hyperbolic functions are used.
Hence, to overcome overflow occurrences, appendix B is dedicated to the establishment of five possible limiting values (Ω lim J for J = 1, 2, 3, 4, 5) for the non-dimensional frequency Ω defined by Eq. ( 11).Accordingly, each member of the model will have its own limiting frequency given by: For a framed structure, the overall limiting frequency is the smallest value given by Eq. ( 37) applied for all the members.If this overall frequency is less than the maximum desired frequency, the only remaining alternative is to promote mesh refinement.As inferred from the example in Fig. 2, the shorter the elements the greater the overall limiting frequency.
On the other hand, the more refined the mesh the greater the possibility of overflow occurrence in the calculation of the DSM determinant because the matrix order will be greater.Fortunately, as shown afterwards, scaling operation in the DSM can be employed to avoid such an occurrence.
In order to obtain a more elucidative visual effect, the example in Fig. 2 uses an intentionally low R max value.In Fig. 2a, it can be observed how each of the non-dimensional limiting frequencies Ω lim J varies with the mesh refinement.By applying Eq. (37), Fig. 2b shows the behavior of the overall limiting frequency f lim .Finally, Fig. 2c shows how mesh refinement increases the modal order of the greatest natural frequency that can be calculated.
Table 1 explores how the increase in the order of R max , up to a compatible value with the real computer capacity, affects the modal order of the greatest natural frequency that can be calculated.In general, it is assumed that only unrealistic data should cause relevant computational limitation by overflow error.For example, from the situation in the last line in Table 1, if an unrealistic axial load P = +10 100 ton is applied, the maximal calculable mode falls to an order of 4.

Dynamic stiffness scaling
We consider the dynamic stiffness matrix, which is defined in Eq. ( 30), now affected by a real and positive scaling constant α G : so that the calculation of the DSM determinant can be written as: where M and N are positive integers, which represent, respectively, the number of elements of the model and the total number of nodal degrees of freedom.Thus, we note from Eq. (39) that the scaling constant α G can be used to change the magnitude of the determinant result without, in essence, affecting the eigenvalue problem.
With no changes in the mode shapes and natural frequency results, a scaling constant α G lower than the unit can avoid the overflow occurrence.
To predict a suitable value of α G , we apply a process of trial and error, in the following terms: 1.The overflow condition is tested for a value less than R max , through the appropriate adoption of the safety coefficient S F > 1: 2. In the calculation of the first natural frequency, while the determinant is greater than V max , the α G value is diminished and the iterative process is restarted.
3. Once a natural frequency is obtained, the maximum current determinant value (D max ), if greater than V max , is used to diminish the current scaling constant by using: For some methods other than the secant method, like bisection, scaling of the dynamic stiffness determinant could be conveniently performed during calculation by maintaining its magnitude within an allowable range and storing an additional integer exponent, so avoiding the need of predicting a suitable value for α G .

POWER SECANT METHOD
The objective here is to elaborate a method that serves as an alternative to the tracing method of Ref. [4], leading to a search for the natural frequencies that are computationally more efficient.
We begin with the discussion of the simple secant method; afterwards, we address the question of polynomial deflation, and, finally, we deal with power deflation.
In principle, the simple secant method lends itself to the extraction of only one root of a function.If the function has only one root, the method is directly applied, and it does not Latin American Journal of Solids and Structures 7(2010) 307 -333 need deflation.If, on the other hand, the function has more than one root, the method can still work under the following conditions: • If the starting points are chosen below the first root, convergence to this root will occur; • If the starting points are sufficiently close to a given root, convergence to this root will occur.
In the case of functions with several roots, we can apply the deflation technique to extract them in ascending order [1,13].Thus, by removing the first root by deflation, if the starting points are chosen from below the second root, convergence to this second root will occur.In succession: being known the (n-1 ) lower roots, if by deflation these roots are removed, when choosing starting points from below the nth root, convergence to this root will occur.

Simple Secant Method
Observing the secant method in its simplest form in Fig. 3, we can develop through recursive application of the equation that represents the intersection of the secant line with the λ axis: Beginning with k = 1, so that it is necessary to arbitrate two starting points λ 1 and λ 2 , the iterations are interrupted when a tolerance test is met, as prescribed by: where d is the number of significant digits that are desired for the approximation to the sought root.If during the iterations f k = f k+1 occurs, the method fails in an unrecoverable way.In general, the method should not be directly applied in case the function has null derivative points, as when more than one root exists.In this context, the deflation technique becomes relevant.

Polynomial deflation
Lets suppose that the function f (λ), from which we desire to extract the roots, is a polynomial of order N, with N real and distinct roots.Thus, the function is known in an explicit form from its coefficients a n : or the function can be imagined in the implicit factorized form, described by monomial products involving its roots λn : As shown before, the simple secant method can be applied for calculation of the first root if the starting points that are lower than this root.Once the first root λ1 is obtained, a new function is found by using the following deflation: Now, the first root of the deflated polynomial F 2 (λ)is in the position of the second root of the original polynomial f (λ).Thus, the secant method can be now applied on the deflated polynomial.If the starting points are chosen below this first root, the method converges to the second root.
Such a deflation procedure can be generalized for the progressive calculation of all the roots.
where F n (λ) is the deflated form of the original polynomial that makes the secant method to converge to the n-order root, being necessary that all the lower (n-1 ) roots have already been calculated.
For the purpose of the eigenvalue extraction of a linear problem, and for the conditions that have been established here, the secant method with polynomial deflation can be applied with success even for the case of repeated natural frequencies [13].

Power deflation
We have seen in the previous section that the secant method with polynomial deflation works perfectly well for the cases in which the function is a polynomial with real roots.In this context, for linear and non-transcendental eigenvalue problems, similar to the problems that Latin American Journal of Solids and Structures 7(2010) 307 -333 arise from FEM models, the secant method can be employed without reservations.Being K the stiffness matrix and M the mass matrix of a FEM model with N degrees of freedom: is a polynomial of order N with a maximum of N real and non negative roots [1].The fact that the polynomial in Eq. ( 46) is not known in an explicit form, like in Eq. ( 42), does not hinder the application of the secant method because it is still possible to make numerical evaluations of f (λ).Especially for case of Eq. ( 46), such numerical evaluations are more efficiently conducted through the Gauss factorization [1]: where L(λ) is a lower triangular matrix and D(λ) is a diagonal matrix.Once performing factorization, the determinant is simply calculated from the products of the diagonal terms of the matrix D(λ): On the other hand, for the improved dynamic stiffness matrix Ψ in Eq. ( 30), the function f (λ) = det{Ψ(λ)} is not a polynomial, although it has non-negative real roots.Thus, the polynomial deflation has little or no possibility of working.
In these terms, alternatively we use the following modified version of Eq. ( 45), which characterizes what we call power deflation: where each N j power should be determined in a manner that, given two starting points λ 1n and λ 2n , the following condition is met: For calculation of the first eigenvalue (n = 1 ⇒ j = 0), we have F 0 (λ) = f (λ) and λo = 0. Consequently, before starting the iterations for each root, it is guaranteed that the first secant line favorably intersects the λ axis at a point located ahead the starting points.Then, from Eq. ( 50), we can obtain the following general expression for calculating the N j powers that are applied in Eq. ( 49): Latin American Journal of Solids and Structures 7(2010) 307 -333 If this expression results in N j < 1, N j = 1 is adopted because the corresponding deflation is not needed.
For the iterations of the first mode (n = 1 ⇒ j = 0), we calculate N 0 using Eq. ( 51) and apply Eq. (49) in order to obtain F 1 to use in the secant method.For the iterations of the second mode (n = 2 ⇒ j = 1), using F 1 , we calculate N 1 by Eq. ( 51) and apply Eq. ( 49) in order to obtain F 2 to use in the secant method; and thus, the other following modes are progressively calculated.
To illustrate the operation of the secant method with power deflation, or here simply called "the power secant method", we present in Fig. 4 the results for a simple clamped-clamped Euler-Bernoulli beam.In this case, we used ζ = 2 and appropriate starting points in order that convergence occurs for up to ten iterations.In the figure, we can visualize how the deflated functions converge for the eigenvalues of the problem, which are localized at the peaks of the function − log{|f (λ)|}.Obviously, the success of the method, in terms of the number of iterations that are needed to achieve convergence, strongly depends on the starting points and on the appropriate choice of the parameter ζ.
Depending on the problem, the parameter ζ can be previously admitted as a constant or varied, by trial and error, until the iterative process demonstrates appropriate convergence.In this last case, in the beginning of the iterations for each root, ζ should be re-estimated whenever necessary.
In Table 2, the method is applied to the ten possible distinct cases of boundary conditions for a simple Euler-Bernoulli beam, for which a total of six different frequency equations exist.

Latin American Journal of Solids and Structures 7(2010) 307 -333
Without the benefit of prior knowledge of the roots, the solution for each one of these equations uses arbitrary starting points as well as automatically determined value for the parameter ζ.
We note in Table 2 that, under this more general situation useful for computer programming, there is a reasonable increase in the number of iterations.We also observe that the final results precisely coincide with the known solutions of the problems [2].
In spite of the number of iterations, the power secant method is always faster than the bisection method.For comparison, lets take the example of the clamped-clamped beam in Table 2.The power secant method performed 168 evaluations of the frequency equation while a bisection method needs more than 1000 evaluations to reach the same accuracy.

COMPARISON EXAMPLES
The power secant method was implemented in the VIGENE program [3], which was developed on the MATLAB platform.This program is a part of the VIANDI computational system that is available for download at http://www.gmsie.usp.br .Section 7.1 presents a very simple example, which has analytical solution.Subsequently, sections 7.2 and 7.3 explore other simple examples, with results from the literature [7,8,12] also compared with additional solutions by FEM where, by dividing the members into hundred equal parts, the mesh was deemed appropriate for the comparisons.
Finally, in order to explore the entire VIGENE capability, section 7.4 presents a more complex model.

Simply supported Euler-Bernoulli beam under axial load
As shown in the appendix A, if shear deflection and rotatory inertia are disregarded, the simply supported beam has exact analytical solution, which does not require application of any iterative numerical processes.We can observe in Table 3 the perfect concordance between this analytical solution and the numerical solution of the VIGENE program.

Single Timoshenko beam for various boundary conditions
For comparison, from [7] and [12] we took the examples shown in Tables 4 and 5, respectively.A good concordance can be observed between VIGENE and ANSYS programs.

Building like frames comprising Timoshenko beam members
For the examples in Fig. 5 and 6 from [8], all the members have the same material and section properties, as indicated in Tables 6 and 7, respectively.An excellent agreement for the results of both examples can be observed, even when compared with FEM solutions by ANSYS.

Robotic handler comprising Timoshenko beam members
In order to explore all the potential of the present work, Fig. 7 depicts an actual example of a robotic handler where effects like rigid offset, end force release and skewed edge are employed.Table 8 gives section and material properties for all the members while Table 9 shows the obtained natural frequencies.The results of both VIGENE and ANSYS programs agree very well but it must be pointed out that the ANSYS model needed more than 240 beam elements to achieve the same accuracy obtained by the exact solution of the program VIGENE, which uses only four elements.

CONCLUSIONS
This work consolidates the possibility of using an iterative process for extracting eigenvalues, without the compulsory use of the Wittrick-Williams algorithm [6,10,11,14].Ref. [4] gives more details of the reasons why such an algorithm has operational difficulties when dealing with models that incorporate the special effects of end release and rigid offset.
For the examples here presented, the results in Table 10 show how an iterative determinant search method, based on the improved DSM of section 4.2, can be more computationally efficient.In terms of CPU time, the power secant method is always faster than a bracketing method based on the Wittrick-Williams algorithm.Therefore, as a well-known fact [1], the secant method is faster than any suitable determinant tracing, bracketing or bisection method.For CPU time comparison with FEM solutions, it must be considered that the FEM models, in general, cannot give accurate results unless fine meshes are employed.In the example of Fig 7, for instance, more than 240 beam elements were necessary in order to achieve a good accuracy.Having this in mind, Table 10 shows that the power secant method can be faster than the subspace iteration method [1] employed by the ANSYS program.When compared with FEM, however, the main advantage of the power secant method based on the improved DSM is the capability to find exact solutions that do not depend on the mesh refinement.
Finally, it must be emphasized that the Wittrick-Williams algorithm must use the classical formulation for the DSM while the power secant method must use the improved formulation as defined in section 4.2.As an additional advantage, the improved DSM preserves any local mode shape for which the nodal displacement vector Q o is null.Any other method based on the classical DSM cannot appropriately detect local modes since the classical DSM becomes ill-Latin American Journal of Solids and Structures 7(2010) 307 -333 conditioned.Moreover, using an inverse vector iteration method based on the improved DSM, it is possible to detect all the distinct mode shapes even when repeated natural frequencies exist.

FFigure 4
Figure 4 Application example of the power secant method (ζ = 2): calculation of the first three roots of the frequency equation f (λ) = cos(λ).cosh(λ) − 1, for a clamped-clamped single Euler-Bernoulli beam.The marks on the curves show the performed iterations.

Table 1
Maximum obtainable modal order for the example of Fig.2.

Table 2
Power secant method applied to single Euler-Bernoulli beams.>= 2 both starting points must be greater than the previous root c Roots obtained with six digits of tolerance a Boundary condition convention: P=Pinned; C=Clamped; F=Free; S=Sliding.bFor j
a Mode where maximum error occurs b Maximum relative error in percentage
a Mode where maximum error occurs b Maximum relative error in percentage

Table 8
Section and material properties for the robotic handler of Fig.7.

Table 9
Natural frequencies [Hz] for the robotic handler of Fig. 7. Mesh size = 1 mm; Element types: BEAM3 for members 3 & 4, BEAM44 for members 1 & 2 and BEAM54 for the top of members 1 & 2; Consistent mass matrix. a

Table 10 CPU
(Intel Core2 E7400 2.8 GHz) Time[s]for the examples of section 7.

Table a
Case b For all the cases, the first 20 natural frequencies were calculated; c Bracketing method based on the Wittrick-Williams algorithm which uses the classical DSM; a Previous table with natural frequency results; b d Iterative secant method using the improved DSM;e Subspace iteration method.