A dynamic branch-switching method for parametrically excited systems

The branch-switching algorithm in static is applied to steady state dynamic problems. The governing ordinary differential equations are transformed to nonlinear algebraic equations by means of harmonic balance method using multiple frequency components. The frequency components of the (irrational) nonlinearity of oscillator are obtained by Fast Fourier Transform and Toeplitz Jacobian method (FFT/TJM). All singularities, folds, flips, period doubling and period bubbling, are computed accurately in an analytical manner. Coexisting solutions can be predicted without using initial condition search. The consistence of both stability criteria in time and frequency domains is discussed. A highly nonlinear parametrically excited system is given as example. All connected solution paths are predicted.


Introduction
A nonlinear dynamical system can be described accurately by a set of nonlinear ordinary differential equations in time variable when the spatial variables are eliminated by means of the finite element method.It is well known that a nonlinear dynamical system will have multiple solutions [9].The number of solutions is not known before hand and may change as the system parameters vary.There are qualitative and quantitative studies of the solutions of a dynamical system in its parameter space.For steady state numerical solutions, there are time integration methods [32], frequency methods [13,30,39] and time-frequency meth-ods [1].The time methods are classified by the ways of replacing the time derivatives by differences: implicit, explicit and Runge-Kutta methods.The time methods are simple to implement and can predict many kinds of responses (superharmonic, subharmonic and chaotic) reasonably well.However, at singularities, they do not converge and the trajectories look like chaotic.For higher order subharmonic solutions, the time step must be very small so that it takes exceedingly long time to reach the steady state and it will be very difficult to distinguish between round off errors and chaos.All possible combinations of initial conditions must be tried in a blanket search manner for a set of system parameters to find all possible steady state solutions.Some sensitive and/or unstable solutions can be missed easily.
Alternatively, the frequency methods replace the single time variable by many frequency parameters and the differential equations are transformed into algebraic equations, which are much easier to handle than the former.Newtonian method is very popular in solving nonlinear equations.The incremental harmonic balance (IHB) method [13] linearizes the ordinary differential equations before using the harmonic balance method.The resulting incremental algebraic equations are linear and therefore, the original IHB method can not study bifurcation and chaos [4].To find subharmonic responses the initial trial solution must be obtained by other means, e.g., time method or the singular perturbation method.Leung and Fung [16] apply harmonic balance to the ordinary differential equations first to obtain a set of nonlinear algebraic equations which are linearized by the Newtonian method afterwards.Higher order dynamic period doublings and chaos are predicted accurately and are proved to be consistent with those of Ueda [38] by analog simulation.The frequency methods rely on the possibility of expanding the nonlinearity in terms of the Fourier series.They can not handle irrational functions.Lewandowski [27,28] made the method in a systematic manner.
The first author and his colleagues have studied the frequency method extensively.Leung and Fung [17] Shock and Vibration 6 (1999) 183-196 ISSN 1070-9622 / $8.00 © 1999, IOS Press.All rights reserved extends it to nonlinear steady state vibration of frames by finite element method.They applied the method to study dynamic snap through [18], spinning structures [19], elasto-plastic systems [20], two-harmonic excitation [21], and dynamic sub-structuring [22].Nonlinear modal analysis was suggested [26].Fergusson and Leung [3] gave an explicit formula to find the harmonic balance matrix using matrix manipulation.Leung and Chui [5] studied coupled Duffing oscillators and the nonlinear vibration of von-Karman's square plate [14].Leung and Ge [23] discovered that the Jacobian matrix in the harmonic balance is actually a Toeplitz matrix which can be represented by a few vectors.Computational efficiency is greatly improved.Together with FFT technique, the method is called the FFT/TJM (Toeplitz Jacobian Matrix).The FFT/TJM method was applied to discontinuous oscillators [6], nonlinear multi-modes of Euler beams [24], double-period excitation [7], and construction of invariant torus [8].
The time-frequency methods [13] are essentially frequency methods except that the Fourier expansion of nonlinearity is replaced by the Fourier transform of the nonlinear functions which are evaluated in discrete time.One of the time-frequency methods is called the alternating frequency time (AFT) method [13], which has three major drawbacks.(i) The Fourier series is simply truncated at a prescribed number of terms and the leading terms are not necessarily balanced.(ii) The Jacobian is approximately given by difference quotient which is not accurate enough to locate the singularity point.(iii) A cut-off fraction is pre-set to determine the significant terms to be included in the computation to reduce the number of unknowns.The method is not conformed to the spirit of the Liapunov-Schmidt reduction and important nonlinear phenomena associated with the neglected terms can be missed easily.
The dynamic branch switching method presented is also a time-frequency method.As the system parameters vary, our aims are (i) to follow the existing solution branch; (ii) to find the singular points; and (iii) to switch to a new branch when the exiting branch becomes unstable at a singular point.Our method is different from the AFT method in that (i) the Fourier series is not simply truncated but rather averaged so that the retaining terms are balanced in Galerkin's sense; (ii) the Jacobian is also averaged analytically which is possible for piecewise continuous functions; (iii) branch switching strategy is adopted when the original path becomes unstable near the bifurcation points; and (iv) the stability of the solution is determined by Floquet-Liapunov theory.
Parametrically excited nonlinear systems have been studied extensively.Watt and Cartmell [41] studied an externally loaded parametric oscillator by multipletime-scale.Szabelski and Warminski [36] investigated a self-excited system with parametric and external excitation by harmonic balance using two terms approximation.Kahraman and Blankenship [12] used harmonic balance to find the interactions between commensurate parametric and forcing excitation with clearance.Sanchez and Nayfeh [35] investigated the global behavior of a biased nonlinear oscilltor under external and parametric excitations by varying the bias parameter and solved the resulting equations by analog computer, FFT and Poincare maps.We use the FFT/TJM method to study the parametrically excited system until chaos occurs.
An equation representing a parametrically excited buckled beam will be taken as an example.Period doubling and period bubbling (period doubling and then period halving) curves are determined accurately.Sudden change from deterministic to chaotic response within 0.5% change of a system parameter is predicted without difficulty.The complete response curve for a parametrically excited system is given for the first time.

Basic formulation
Consider a scalar ODE in the form: where f is an analytical function, t is the time, x is the unknown response, Ω is the fundamental frequency of the linearized system, F and g are the amplitude and time function of external force, respectively.Over dot represents differentiation with respect to time () = d()/dt = Ωd()/dτ in which τ is a non-dimensional time factor, τ = Ω • t.Solution of Eq. ( 1) can be expressed in a Fourier series where a 0 , a i and b i are the unknown Fourier coefficient and L is the order of subharmonics.If one considers the solution precise to subharmonic order 3, then L = 3.It will be shown later that chaotic motion via period doubling occurs when the periodic solution bifurcates from L = 4 to L = 8 and bifurcates to higher order subharmonic at an accelerated rate.Set increments δx, δΩ, δF and expand Eq. ( 1) in Taylor series near x 0 , Ω 0 , F 0 which satisfy Eq. ( 1), where HOT denotes the higher order terms.The incremental equation is in which and r(τ ) is the residual given by The solution of the incremental equation satisfies the condition of harmonic balance, so that The neighbouring state is then reached by the Newtonian method when the residual of the approximate solution approaches to zero.The other step in the IHB is the formation of the Jacobi matrix for Newtonian iteration by the Galerkin method Truncating the Fourier series in Eq. ( 5) to a finite number of terms N , and substituting Eq. ( 4) into Eq.( 8), we have where where The implementation of Eqs ( 11) and ( 12) requires analytic calculation of the closed form integrals, which restricts the method from many applications.

Time discretization on IHB
An improvement is now presented using windowing and sampling techniques to discretize the continuous integration in Eqs (11) and (12).The Galerkin averaging is achieved most efficiently using the Fast Fourier Transformation (FFT).The results obtained are exactly the same as those found by Eqs (11) and (12).Compared with the method described in the previous section, a very small amount of analytical work is required before iteration.The alternative transformations from the frequency domain to the time domain are included in order to make use of the advantage of each domain.A fast transformation pair of real form based on the standard FFT algorithm [31] is developed in this scheme to bridge the response between the time domain and the frequency domain, the basic algorithm being provided as follows.The order of subharmonic is assumed to be L.
a) The backward transformation from a 0 , a i , b i to its time domain responses x k The complex form of Eq. ( 13) being where The standard inverse FFT algorithm can solve Eq. ( 14).
b) The forward transformation from x k to a 0 , a i , and b i x k , The first and second terms in a i and b i are corresponding to the complex coefficients c i and c M−i from the results of the standard forward FFT.It should be noted that the discrete Fourier transform implies periodicity in both the time and frequency domains.The correct application of the sampling theorem [31] for M 2N + 1, where M is a power of 2, must be satisfied where with The Jacobi matrix is rewritten as where The transformation formulas to find J cc oj , J cc ij , J sc ij are A kj , Similar discretizations are also available for J cs oj , J cs ij , and J ss ij .Therefore, we can find the solution {Q  14) and (15).It was found [23] that the Jacobian matrix thus obtained is actually a Toeplitz matrix which can be represented by a few vectors.Efficient computational method can be derived.

Double pivoting arc-length method
Note that in Eq. ( 10), the number of incremental is more than the number of equations available due to the presence of δΩ and δF .One of the increments must be fixed so that the other increment can be solved from Eq. (10).Therefore, we can rewrite Eq. (10) to where G ∼ is an augmented matrix with order of 2N + 1 by 2N + 2 and δλ is called the control increment or the active increment, which is either δΩ or δF .If the active increment vanishes, an improved solution with diminishing residual vector {R ∼ } is guaranteed by the nature of the Newtonian method.The process is repeated until the magnitude of the residual vector {R ∼ } is acceptably small and a solution (an equilibrium state) is obtained.This process is called iteration.On the other hand, if the value of the active increment is augmented artificially from a known equilibrium state to find a new solution, the process is called an augmentation.An alternative application of augmentation and iteration is an useful strategy in the IHB scheme.The method mentioned above is successful in predicting the dynamic response curves in general.However, for typical nonlinear stability phenomena, this strategy is not sufficient since the Jacobian is in the neighborhood of turning points and bifurcation points.Furthermore, when the system damping is very small, it may also fail to converge near the sharp peaks of the response curve.This imperfection is usually improved by appending an additional equation to the solution scheme.A well-known example of solving strategy is the socalled arc-length method as widely discussed in nonlinear FEM field by Risk [34], Crisfield [2], Wagner and Wriggers [40].The basic idea of such methods is to append an additional equation to the original set of Eq. ( 10) constraining the active increments by a given arc-length.Although this method has been successful in nonlinear static FEM, it needs some modifications when applies to dynamic problems [29].A double pivoting arc-length method is applied here to avoid the singular Jacobian at turning points and bifurcation points.Denote X n+1 = λ and enlarge the vector {X ∼ } from order 2N + 1 to order 2N + 2 to include λ.Here Eq. ( 24) can be written as where {δ We pivot the active increment on the fastest changed variable [25] in the last augmentation step.Assume that the active increment in the current iteration is δX k .The incremental linear Eqs (10) can then be solved for the 2N + 1 unknowns β i , where β i are scale factors whose function is similar to directional derivatives, The remaining unknown derivative δx k /δl is solved from the constraint Eq. ( 27), which can be rewritten as The subjective Jacobian {G ∼ k } of Eq. ( 25) will be regular for all turning points since we have chosen the most fastest changed variable as active at each augmentation step.The solution of the iteration are expressed as: where {X ∼ k } 0 is the initial vector for iteration.

Stability points and secondary branches
It is well known that the physical response curve will switch to the new bifurcated branch when the stability of original branch is lost, or one can only get an unstable solution with its amplitude growing monotonically [5].The cascade diagram, which represents a sequence of period doubling, is the most celebrated example of stability loss.Although the algorithm of IHB can take practically all the harmonic terms into consideration, it is shown by numerous computation experiments that subharmonic responses can never be obtained automatically without using branch-switching algorithm.The emphasis in the following discussion will focus on the computational aspects.Our objective is to develop a numerical scheme which is consistent with the context of our IHB scheme mentioned above for bifurcation problems of period doubling and period halving.
A point is a stability point when where {J ∼ } is the enlarged Jacobian taking the double period response into consideration.The calculation of the double period branch switching requires additional consideration.Near a stability point, the dimension of equations is enlarged so that subharmonic terms are accommodated.And then the eigenvalues of the Jacobian are used to find the umber of possible branches near the bifurcation point.The number of branches is determined by the number of zero eigenvalues λ i of Eq. ( 30).
The eigenvectors Φ i indicate the directions of the secondary branches associated with λ i = 0.The branch switching is performed by substituting the factors β i in the Eq. ( 28) with the scaled eigenvector giving the tangent vector along the secondary branch in which {X ∼ i } 0 is the nearest bifurcation points and superscripts I, II represent the solution branches.The pivoting of the control increment depends on the character of bifurcation points.For instance, in case of symmetry breaking the most effective pivoting could be on the first element in vector {δX ∼ i } while for period double bifurcation, pivoting on the one of the coefficients of the associated subharmonic terms will successfully lead to the solution on the secondary branch.The value of {δX ∼ k } can be approximated as This value is critical for a successful branch switching.

Floquet-Liapunov theory
When a solution is obtained, its stability characteristic can be checked by the Floquet-Liapunov theory.
An approximate method based on this theory was developed by Hsu [11], and Hsu and Cheng [10] who discretized the state transition matrix by a series of step function and uses a single pass numerical integration to a small perturbation.The basic idea is presented as follows.When the solution is perturbed by δx, the autonomous linear incremental equation is where in which M (t), C(t), K(t) are time dependent functions.Rewrite Eq. ( 34) as where {u Discretize the time domain and replace [A ∼ (t)] with step function: where the sampling period ϑ = T /M and M is the number of discrete points in the defined period.The perturbed system can be approximated by Then the stability of Eq. ( 34) is checked by evaluating the eigenvalues of the transition matrix [Q ∼ ] at the end of the specified period from time t = MT to time t = (M + 1)T .The two eigenvalues are so called Floquet multipliers.If the absolute magnitudes of the Floquet multipliers are less than unity, the solution is stable.For any small perturbation δx, its oscillation is diminishing.The explicit form of [Q ∼ ] can be written as The two eigenvalues are given by where For dissipative systems, the solutions are stable when the moduli of the complex eigenvalues are within the unit circle.When the system parameter changes, instability occurs if one of its eigenvalues leaving the unit circle.That a positive real eigenvalue leaves the unit circle while the other remains inside predicts a limit point in folded solution curves.The solution on the folded part is shown by a shift of the equilibrium position.When a negative real eigenvalue leaves the unit circle through −1, the point is known as flip bifurcation point and an unstable one period solution will bifurcate into a stable double period solution.This phenomenon is most commonly met at bifurcation points and happened in the system we considered below.When a pair of complex conjugate eigenvalue leaves the unit circle, it is called Hopf bifurcation which has been discussed by many authors.

Numerical example
As an example, we consider a class of problem whose excitation appears as a coefficient in the equation of motion.There is a number of earlier studies on chaos in parametrically excited systems related to the forced motion of buckled beams [33,37].The resulting equation can be written in terms of x as Fig. 1, Fig. 1.A parametrically excited oscillator.
For the numerical solution of Eq. ( 40), we have chosen c = 0.2, ω 0 = 0.5 and α =1.From the view of parametric study, this choice corresponds to a continuous change of the amplitude of external force.The linearized incremental equation in time domain is obtained as: where r is the unbalanced residual in the current iteration which will vanish when a solution is reached.The related incremental equation in frequency domain are then formulated through FFT/TJM as the same as that presented by Eq. (10).The stability characteristic can be studied by Floquet theory, with which is periodic with period T = 2πL.The number of harmonic terms and residual tolerance used in each iteration depends on the required order of the subharmonic solutions.The total number of terms must be increased when large unbalanced residual terms are detected.We observe later that the sequence of period double bifurcation occurs at closely neighbouring points.The evaluation of bifurcation and branch switching is possible only if the unbalanced residual can be kept small.in particular, for period 4π orbits, we have used 2N MAX + 1 = 57 terms and for period 8π orbit 2N MAX + 1 = 113 terms, which are able to keep the residual less than 10 −6 and 10 −8 , receptively.

Results and discussion
For small values of q (q < 0.28637), the solutions are attracted only by the stable fixed point x = 0, x = 0 with two Floquet multipliers less than unit in period 2π.It should be noted that our method eliminates the interference of the initial conditions and connect all possible fixed point of the equation with one pseudo arc-length curve.So no additional efforts are required to predict the three coexisting solutions of Eq. ( 42), say at q = 0.6.Figure 2 shows the folded curve, there are three solutions, two stables and one unstable, in some area.At point a, the original period 1(2π) solution loses stability with one of its Floquet multipliers leaving the unit circle through −1.The stable period 2(4π) solutions could be obtained by using dynamic branch switching.The period 1 solution regains its stability at point b where the period 2 solutions merge with period 1.A period bubbling is completed.At point c we have a positive real Floquet multiplier exceeding unit which relates to a fold or a saddle node bifurcation.The solution is then attracted by the unstable fixed point x = π, x = 0, tracing along the upper unstable part of the curve with period 2π.The curve folds back again at point d where the period 1 solution is attracted the fixed point x = −π, x = 0 and all Floquet multipliers stay in the unit cycle.Period doubling cascade happens when q reaches 0.2869544.We combine the arc-length and branch-switching algorithm to trace the stable branches.Figure 3 shows the period 2π stable branch bifurcating into period 4π.At point f the branches bifurcates again to 8π, and at point g the stable solution period bifurcates into 16π (Fig. 4).Further period double bifurcations on the stable branches are              If the determinant of the Jacobian in Eq. ( 10) takes the double period into consideration, it is called the enlarged determinant.The zero points of the enlarged determinant are correspondent to the points where the Floquet multipliers pass through the unit cycle.This consistent relation in the determination of solution stability is very clear from theories and also verified by our numerous computational practices.The benefit from the enlarged determinant check is that additional calculation to obtain Floquet multipliers of solution is avoided.So based on this we can use the inverse iteration to find the position of local bifurcation point.We

Conclusion
A numerical method for determining the steady state forced vibration of a highly non-linear system has been presented.The method in which an alternative FFT/TJM algorithm is used, is comparatively easy to formulate the harmonic balance equations as compared to that of IHB.An effective algorithm to search for the specified stable branches and characterise the stability of solutions has been described.The newly developed method was applied to a forced parametric excitation system.A number of well-known nonlinear phenomena such as fold and period double bifurcation is obtained.The results generated here are coincided with that of direct numerical integration.The computer time is greatly reduced and the response diagrams are neatly constructed.