Valery N Pilipchuk

Application of Special Nonsmooth Temporal Transformations to Linear and Nonlinear Systems under Discontinuous and Impulsive Excitation

Keywords: Pulsed excitation, parametric loading, nonsmooth transformations, nonsmooth dynamics

Linear and nonlinear mechanical systems under periodic impulsive excitation are considered. Solutions of the differential equations of motion are represented in a special form which contains a standard pair of nonsmooth periodic functions and possesses a convenient structure. This form is also suitable in the case of excitation with a periodic series of discontinuities of the first kind (a stepwise excitation). The transformations are illustrated in a series of examples. An explicit form of analytical solutions has been obtained for periodic regimes. In the case of parametric impulsive excitation, it is shown that a nonequidistant distribution of the impulses with dipole-like temporal shifts may significantly effect the qualitative characteristics of the response. For example, the sequence of instability zones loses its different subsequences depending on the parameter of the shifts. It is shown that the method's applicability can be extended for nonperiodic regimes by involving the idea of averaging.

Introduction

Instantaneous impulses acting on a mechanical system can be taken into account either by imposing the special conditions of matching for the velocities and coordinates of the system at points of the impulses' location or by introducing Dirac's functions into the equations of motion. The first approach, in principal, deals with the differential equations for an unforced system [START_REF] Samoilenko | Differential Equations with Impulsive Excitation[END_REF]. However, these equations are treated separately in each of the intervals between impulses and, hence, instead of a single system, a whole sequence of systems must be considered. The second method gives a single set of equations over the whole time interval without any conditions of matching imposed on the variables. However, in this case, analysis can be carried out correctly within the framework of the theory of generalized functions (distributions) [START_REF] Richtmyer | Principles of Advanced Mathematical Physics[END_REF], which requires additional mathematical proofs in nonlinear cases [START_REF] Maslov | Asymptotic solution of equations with a small dispersion[END_REF][START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] (see this book also for a more complete list of publications regarding systems with impulses). Both of the above-mentioned ways of modeling are fruitfully employed for different quantitative and qualitative considerations of mechanical systems under pulsed excitations [START_REF] Samoilenko | Differential Equations with Impulsive Excitation[END_REF][START_REF] Hsu | Impulsive parametric excitation[END_REF][START_REF] Hsu | Application of the theory of impulsive parametric excitations problems[END_REF][START_REF] Faure | Percussions en mécanique non linéaire sur certaines solutions périodiques de phénomènes non linéaires excités par des percussions[END_REF][START_REF] Faure | Percussions en mécanique non linéaire: (I) Cas des interractions entre systèmes. (II) Théorie des percussions presque périodiques[END_REF][START_REF] Faure | Théorie des oscillations paramétriques cas des excitations par des percussions presque périodique sur un exemple de système couplé[END_REF]. A method is described below which enables one, on the one hand, to eliminate the singular terms in the equations and, on the other hand, to obtain solutions in the form of a single analytic expression over the whole time interval. Originally this method was proposed for a class of smooth vibrating systems. Here, the term 'smooth' means that the motion is described by means of smooth (differentiable) functions of time. However, a special feature of the technique is that the smooth solutions to be defined are expressed through a pair of nonsmooth functions such as a saw-tooth sine and rectangular cosine. Mathematical properties of the related transformation can be briefly listed as follows [START_REF] Pilipchuk | A transformation of vibrating systems based on a nonsmooth periodic pair of functions[END_REF]:

• an introduction of the nonsmooth (saw-tooth) time transforms the coordinates into elements of two-components hyperbolic numbers algebra1 ; • under certain conditions, the result of the differentiation of the elements remains in the algebra, and the result of integration remains in the algebra; • a slow timescale can be introduced into the transformation to describe the nonperiodic processes. Based on these points, the introduction of a nonsmooth independent variable into differential equations of motion can be justified. Different physical and mechanical applications and treatments of saw-tooth temporal transformation were discussed in [START_REF] Pilipchuk | Analytical study of vibrating systems with strong non-linearities by employing saw-tooth time transformations[END_REF][START_REF] Vakakis | Normal Modes and Localization in Non-linear Systems[END_REF]. It should be noted that the transformation itself has to be viewed as a preliminary step of investigation. It reduces an original system to a boundary-value problem with respect to a new independent variable, i.e., the saw-tooth bounded argument. Then, one can invoke the appropriate method to solve the boundary-value problem accordingly to the physical contents of the model. For example, special successive iterations processes and asymptotic techniques were used for different systems performing weakly and strongly nonlinear periodic motions, and modulated regimes as well [START_REF] Pilipchuk | Analytical study of vibrating systems with strong non-linearities by employing saw-tooth time transformations[END_REF][START_REF] Vakakis | Normal Modes and Localization in Non-linear Systems[END_REF]. The shooting method was applied to study a complicated set of periodic orbits close to the separatrix [START_REF] Pilipchuk | Study of a class of subharmonic motions using a nonsmooth temporal transformation (NSTT)[END_REF]. The transformation can also be realized with respect to a spatial independent variable when dealing with different spatially periodic structures [START_REF] Salenger | Discreteness effects in the forced dynamics of a string on a periodic array of non-linear supports[END_REF][START_REF] Pilipchuk | On one variant of nonsmooth transformations of variables for 1-D elastic systems of a periodic structure[END_REF][START_REF] Pilipchuk | Study of the oscillations of a nonlinearly supported string using a nonsmooth transformation[END_REF]. As was mentioned above, this paper focuses on a class of dynamical systems under periodic impulsive and discontinuous excitations. A generalized version of the transformation based on the sawtooth function with an asymmetric slope [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF] will also be considered. This enables us to take into account the dipole-like shifts of singularities of the external force and, hence, to involve a more general class of dynamical systems.

Some Reasons for Nonsmooth Transformations

This section presents some remarks concerning nonsmooth temporal transformations. Some of the basic relations will be reproduced for future references as well. A theoretical reason for the nonsmooth transformations follows from a simple general observation. Indeed, linear transformation of coordinates may significantly simplify a problem of linear dynamics. Weakly nonlinear transformations of coordinates can be viewed as a basic step in most of the analytical techniques in the theory of quasilinear systems. The theory of solitons gives a great example of a strongly nonlinear theory based on strongly nonlinear transformations of variables. Now one can conclude that dynamical systems with discontinuities can be simplified by means of the appropriate nonsmooth transformations of variables. For instance, in the case of systems with absolutely rigid constraints, special nonsmooth transformations have been proposed in [START_REF] Zhuravlev | Investigation of certain vibro-impact systems by the method of nonsmooth transformations[END_REF]. This realization deals with the positional variables (coordinates) of a vibro-impact2 system and enables one to eliminate the constraints that the vibro-impact systems must include by definition. The transformations instantaneously change the direction of the coordinate axes at the instance of impact (a geometrical treatment is given in Figure 1). As a result, the transformed system becomes free from its absolutely rigid barriers and the corresponding differential equations do not contain any 'impact terms'.

Unfortunately, there is no direct way to use nonsmooth transformations of coordinates in the case of the external impacts, due to external force. Another question is how to apply the nonsmooth transformation of coordinates to a broad class of systems which, within a certain meaning, are close to vibro-impact, but do not include any absolutely rigid constrains. A related attempt was made in [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly supports[END_REF] when considering the localized vibro-impact modes in a linear string with concentrated masses between the two rigid walls. To deal with differential equations of motion over all time domains and under no constraint conditions, the rigid walls were replaced by a special potential field with a power function of a higher odd exponent, x m . In this case, the classical vibro-impact situation appears when m -→ ∞, as the rigid walls limit. Is it possible to construct a transformation of the coordinate x(t) -→ s(t) depending on m as x = f (s, m) and, giving nonsmooth transformation by [START_REF] Zhuravlev | Investigation of certain vibro-impact systems by the method of nonsmooth transformations[END_REF], when m -→ ∞? This question has not yet been answered. Probably, it is hardly possible before the solution of the original differential equation is known. However, starting from a study of the asymptotic m -→ ∞, it was preliminarily shown in [START_REF] Pilipchuk | The calculation of strongly nonlinear systems close to vibro-impact systems[END_REF] that for certain kind of systems, a nonsmooth transformation of time can be introduced as

τ = τ (t) : (-∞, ∞) t -→ τ ∈ [-1, 1],
(1)

x(t) = x(τ (t)) ≡ X(τ ), (2)
where x(t) is a periodic function, the period of which is normalized to be equal to four, and τ is a saw-tooth periodic function defined as

τ (t) = t, -1 ≤ t ≤ 1 -t + 2, 1 ≤ t ≤ 3 ; τ (t) ∀t = τ (t + 4). (3)
This function can also be expressed as τ (t) = (2/π ) arcsin sin(π t/2). A possible geometrical treatment of Equation ( 1) is given in Figure 2. It is important to note that transformation (1) is not invertible over all periods of motion T = 4 and, hence, representation ( 2) is not always correct. For any periodic motion, a representation has been obtained in a two-component form.

It is written in terms of both nonsmooth oscillating time, τ , and its first generalized derivative, e = dτ (t)/dt = τ (t) (Figure 3), as [START_REF] Pilipchuk | A transformation of vibrating systems based on a nonsmooth periodic pair of functions[END_REF] x

(t) = X(τ ) + Y (τ )e, ( 4
)
where the functions X and Y are either known (if the original function, x(t), is known), or supposed to be defined as a solution of the transformed equations. Expression (4) can be viewed as a special 'complex' value with real X and 'imaginary' Y components, where the basis 'imaginary' element e generates a circling group: e 2 = 1, e 3 = e, e 4 = 1, . . . . These complex numbers are elements of algebra without division by the basis elements {1, e}. This remark significantly simplifies all mathematical manipulations with the representation when transforming the differential equations of motion. For example,

(X + Y e) 2 = X 2 + Y 2 + 2XY e, (X + Y e) 3 = X 3 + 3XY 2 + (Y 3 + 3Y X 2 )e.
Here, the right-hand sides are of the structure

{• • •} + {• • •}e, as the element X + Y e itself is. A less trivial equality exp(X + Y e) = exp(X)[cosh(Y ) + sinh(Y )e]
(5)

will be taken into account below.

Taking into account that e 2 = 1, one obtains

1 X + Y e = X X 2 -Y 2 - Y X 2 -Y 2 e if X = ±Y.
The last example explains why division is not always possible.

In this paper, we also use a generalized version of representation (4) based on the saw-tooth function with an asymmetric slope:

τ (t; γ ) = t/(1 -γ ), -1 + γ ≤ t ≤ 1 -γ , (-t + 2)/(1 + γ ), 1 -γ ≤ t ≤ 3 + γ , e(t; γ ) = ∂τ (t; γ ) ∂t τ (t; γ ) ∀t = τ (t + 4; γ ) , (6)
where γ ∈ (-1, 1) is a parameter characterizing the inclination of the 'saw's teeth' as Figure 4 shows. The parameter of inclination γ is introduced in such a manner that the period T = 4 of the function and its amplitude value (= 1) do not depend on the inclination. In this case, the algebraic structure becomes more complicated. The above multiplication rule, e 2 = 1, should be replaced with a more complicated one:

e 2 = α(γ ) + β(γ )e; α = 1 1 -γ 2 , β = 2γ α. ( 7
)
In the symmetric case, one has α(0) = 1, β(0) = 0, and ( 7) is reduced to the old multiplication rule.

Properties of representation (4) based on the asymmetric saw-tooth function (6) were considered in previous work [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF]. Note that the result of differentiation of Equation ( 4) formally contains a periodic series of singularities at points = {t : τ (t; γ ) = ±1}. The singularities are discarded by imposing the condition

Y | t ∈ = Y | τ =±1 = 0 ( 8 )
that can be viewed as a necessary condition for the continuity for x(t). Indeed, the stepwise function e(t; γ ) has discontinuities at points . It means that the 'imaginary part' of representation (4), Y (τ )e, would be discontinuous at , if condition (8) was false. Now, let us assume that condition ( 8) is false and one has In this case, function x(t) represented by equality (4) has a periodic series of steps of amplitude 2p/(1γ 2 ) at points . An idea of this work is to employ these discontinuities in order to describe nonsmooth dynamical processes. Note that nonsmooth temporal transformation does not necessarily imply any absolutely rigid barriers. Instead of the coordinates, the parameter of time is transformed by means of the nonsmooth function. As a result, the coordinates acquire a special algebraic structure with a series of mathematically suitable properties. In the next sections, this kind of transformation will be applied to systems subjected to impulsive and discontinuous excitations. Transformation (1-4) instantaneously changes the direction and scale of time when the external impulses act.

Y | τ =±1 = p = 0. ( 9
)
Note that the imaginary basic element e of representation ( 4) is a periodic piece-wise constant step function of time t. In this connection, we note the paper by Liu [START_REF] Liu | Discontinuous and impulsive excitation[END_REF] where the solutions of quasilinear systems under the action of an impulsive or discontinuous force (with a singularity at the point t = a) were presented in the form x(t) = x 1 (t) + x 2 (t)H (ta), where H (t) is Heaviside's step function. An analogous form of the solution is used to describe moving discontinuities in the wave theory in [START_REF] Witham | Linear and Non-Linear Waves[END_REF]. It is known [4, p. 18] that the representation x(t) = y(t) + νH (ta) (ν = const.) transforms the nonlinear system ẋ = f (t, x) + νδ(ta) into a so-called Carathéodory equation with respect to y(t) (see Appendix 1). The Carathéodory form of equation admits a series of existence theorems. A principal difference is that components x 1 (t), x 2 (t) of the above representation depend on the original time parameter, whereas equality (4) provides the introduction of the oscillating time bounded as τ ∈ [-1, 1].

To this end, note that representations (1), (4) and the representation based on functions (6) as well, stay to transform a system on the manifold of periodic solutions. This may bring doubt regarding the practical meaning of the representations. That is why Sections 5 and 6 show how to extend the applicability of the representation. For example, the idea of averaging can be involved, if inserting slow varying parameters into the transformation. Let t • = εt be a slow time scale, where 0 < ε 1. In the case of damped or modulated vibrations, solutions can be sought in the form [START_REF] Pilipchuk | A transformation of vibrating systems based on a nonsmooth periodic pair of functions[END_REF] x

= X(τ, t • ) + Y (τ, t • )e,
where τ = τ (ωt), e = e(ωt) and the constant ω is assumed to be of order 1.

Of course, this representation does not lead to 'another method' of the two variables expansion. However, replacing the fast timescale by oscillating argument τ , one brings some special features into the two-variables expansion method. For instance, a system's fast time oscillation will be described with the boundary-value problem on the standard domain τ ∈ [-1, 1] and the slow time (averaged) equations follow from the boundary conditions at points τ = ±1 but are not a result of the elimination of the secular terms.

Systems with Periodic Pulsed Excitation

We will start from the simple case of a linear system to illustrate an application of the nonsmooth temporal transformation with no additional complications which may be brought about by technical details related to a concrete practical problem. Also, the example itself has a definite practical meaning.

PERIODIC SERIES OF IMPULSES

We will seek periodic solutions of the first-order differential equation

v + λv = μ ∞ k=-∞ [δ(t + 1 -4k) -δ(t -1 -4k)]; λ, μ = const. ( 10
)
For positive λ, the function to be determined, v(t), can be treated as a velocity of a point moving under periodic impulsive force and viscous damping. This equation has a simple form involving the generalized functions. However, the common ways of solving, generally speaking, require first the determination of a class of functions among which the solution should be sought. Then, one can apply either the generalized Fourier series/transforms procedure, or consider the equation between impulses matching pieces of the solution at the points of impulse. The saw-tooth temporal transformation gives us another way to get a closed-form analytical solution.

First, the right-hand side of Equation ( 10) can be expressed in terms of a generalized derivative of rectangular cosine (see Figure 3). As a result, Equation ( 10) is rewritten as

v + λv = μ 2 ė(t). (11)
Representing periodic solutions in the form of (4), v(t) = X(τ (t)) + Y (τ (t))e(t), and substituting this into Equation [START_REF] Pilipchuk | Analytical study of vibrating systems with strong non-linearities by employing saw-tooth time transformations[END_REF], gives

(Y + λX) + (X + λY )e + Y - μ 2 ė(t) = 0.
The first two terms in this expression present an element of structure (4). Both components must be separately equal to zero to provide a complete zero element. 3 This will give us a set of two first-order differential equations. The elimination of the last (underlined) term leads to the boundary conditions. As a result, one obtains the boundary-value problem as

Y + λX = 0, X + λY = 0 Y | τ =±1 = μ 2 . ( 12
)
So the original differential equation [START_REF] Pilipchuk | A transformation of vibrating systems based on a nonsmooth periodic pair of functions[END_REF] has been transformed into boundary-value problem [START_REF] Vakakis | Normal Modes and Localization in Non-linear Systems[END_REF] for two first-order differential equations. The basic advantages, however, is that the transformed system does not include δ-functions any more, and the new independent variable is bounded, |τ | ≤ 1.

Solving boundary-value problem [START_REF] Vakakis | Normal Modes and Localization in Non-linear Systems[END_REF], one obtains the periodic solution of the original equation as

v = X + Y e = μ 2 cosh λ (-sinh λτ + e cosh λτ ).
Taking into account equality [START_REF] Hsu | Impulsive parametric excitation[END_REF] and the expression e 2 = 1, finally gives

v = μ 2 cosh λ exp[-λτ (t)e(t)]e(t). ( 13
)
Figure 5 shows dependencies v(t) for different values of λ. Note that the discontinuous solution v(t) is described by a single expression of the closed form (13) using the two standard nonsmooth periodic functions τ (t) and e(t). These functions admit a clear mechanical treatment as the coordinate and velocity of the vibro-impact oscillator (see Figure 1) under the initial conditions x(0) = 0, ẋ(0) = 1 (see Section 7 for more details).

The boundary-value problem ( 12) is replaced by a more complicated one, when considering the impulsive excitation with a dipole-like shift, i.e. the right-hand side of Equation ( 11) is expressed by a derivative of the asymmetric saw-tooth function with a slope γ [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF]:

v + λv = p ∂e(ωt; γ ) ∂(ωt) = 2p 1 -γ 2 ∞ k=-∞ [δ(ωt + 1 -γ -4k) -δ(ωt -1 + γ -4k)]. (14)
In this case, representation for a periodic solution is taken as v = X(τ ) + Y (τ )e; τ = τ (ωt; γ ), e = e(ωt; γ ). [START_REF] Pilipchuk | On one variant of nonsmooth transformations of variables for 1-D elastic systems of a periodic structure[END_REF] Substituting expression [START_REF] Pilipchuk | On one variant of nonsmooth transformations of variables for 1-D elastic systems of a periodic structure[END_REF] into differential equation ( 14) and collecting separately the 'real' and the 'imaginary' components and singular terms, one obtains

(ωαY + λX) + [ω(X + βY ) + λY ]e + (ωY -p) ∂e(ωt; γ ) ∂(ωt) = 0.
The resulting boundary-value problem is

ω(X + βY ) = -λY, ωαY = -λX, ωY | τ =±1 = p. ( 16
)
Taking the boundary conditions into account, one obtains

Y = p ω cosh γ λ ω cosh λ ω τ cosh λ ω + sinh γ λ ω sinh λ ω τ sinh λ ω exp γ λ ω τ , X = - ωα λ Y . ( 17
)
The X-component is defined by differentiation according to the second equation in [START_REF] Pilipchuk | Study of the oscillations of a nonlinearly supported string using a nonsmooth transformation[END_REF].

HARMONIC OSCILLATOR UNDER PERIODIC IMPULSIVE EXCITATION

We will seek a steady-state periodic solution of the second-order differential equation

ẍ + ω 2 0 x = 2p ∞ k=-∞ [δ(ωt + 1 -4k) -δ(ωt -1 -4k)], (18)
where p, ω 0 and ω are constant parameters.

In the first step, Equation ( 18) is rewritten as

ẍ + ω 2 0 x = p de(ωt) d(ωt) . ( 19
)
We seek a periodic solution of the period T = 4/ω in the form

x(t) = X(τ (ωt)) + Y (τ (ωt))e(ωt). ( 20
)
Substituting this into Equation [START_REF] Vedenova | Normal oscillations of a string with concentrated masses on non-linearly supports[END_REF], under the condition of continuity for x(t) (Equation ( 8)), gives

(ω 2 X + ω 2 0 X) + (ω 2 Y + ω 2 0 Y )e + (ω 2 X -p) de(ωt) d(ωt) = 0.
The first two terms in this expression present an element of the algebraic structure considered. Both components must be separately equal to zero to provide a complete zero element. This will give us a set of two first-order differential equations. Elimination of the last (underlined) term leads to the boundary conditions. As a result, one obtains the boundary-value problem as

X + ( ω 0 ω ) 2 X = 0, Y + ( ω 0 ω ) 2 Y = 0, X | τ =±1 = p ω 2 , Y| τ =±1 = 0. ( 21
)
So the original differential equation [START_REF] Zhuravlev | Investigation of certain vibro-impact systems by the method of nonsmooth transformations[END_REF] has been transformed into boundary-value problem [START_REF] Liu | Discontinuous and impulsive excitation[END_REF] for two second-order differential equations. The basic advantage, however, is that the transformed system does not include δ-functions any more. By solving boundary-value problem [START_REF] Liu | Discontinuous and impulsive excitation[END_REF], one obtains the periodic solution of the original equation as

x = X(τ (ωt)) = p ωω 0 sin ω 0 ω τ (ωt) cos ω 0 ω , Y ≡ 0. ( 22
)
This is a continuous solution for all t, and all possible resonances are given by the roots of the equation cos(ω 0 /ω) = 0 as

ω = 2 π ω 0 n , n = 1, 3, 5, . . . , ( 23
)
where the factor 2/π appeared due to different normalization of the periods for trigonometric and saw-tooth sines.

In the presence of linear viscous damping, boundary-value problem [START_REF] Liu | Discontinuous and impulsive excitation[END_REF] becomes coupled throughout the terms Y (τ ) and X (τ ) in the first and second equations, respectively. However, the problem remains resolvable due to the linearity. Let the equation of motion be written as

ẍ + 2ζ ẋ + ω 2 0 x = p de(ωt) d(ωt) , ( 24
)
where ζ is the damping factor. The boundary-value problem for the X and Y components of the steady-state periodic solution has no discontinuous terms and is written as

X + 2 ζ ω Y + ω 0 ω 2 X = 0, Y + 2 ζ ω X + ω 0 ω 2 Y = 0, X | τ =±1 = p ω 2 , Y| τ =±1 = 0. ( 25
)
Solving the boundary-value problem in a standard way involving the characteristic equation, one finally obtains

x = X + Y e = p βω 2 (cos 2 β cosh 2 α + sin 2 β sinh 2 α) × [cosh α cos β cosh ατ sin βτ -sinh α sin β sinh ατ cos βτ + (sinh α cos βτ cosh ατ sin β -sinh ατ sin βτ cosh α cos β)e], ( 26
)
where τ = τ (ωt) and e = e(ωt); α = ζ /ω and β = ω 2 0ζ 2 /ω. Figure 6 illustrates the qualitatively different responses of the system when varying the external excitation frequency. In different proportions, the responses combine properties of the harmonic damped and nonsmooth motions generated by the impact loading. For example, when ω ω 0 and ω ζ , the system approaches the limit of a free particle under the periodic impulsive force. In this case, the boundary-value problem is reduced to

X = 0, Y = 0; X | τ =±1 = p ω 2 , Y| τ =±1 = 0 ( 27
)
and its solution gives the saw-tooth time shape of the motion, x = pτ (ωt)/ω 2 , which is approached by the time history record in Figure 6c. Now let us consider the N-degrees-of-freedom system Introducing the principal coordinates x j (t) as

M ÿ + Ky = p de(ωt) d(ωt) , ( 28
)
y = N j =1
x j (t)e j , [START_REF] Boettcher | Nonperturbative square-well approximation to a quantum theory[END_REF] one obtains a decoupled set of impulsively forced harmonic oscillators

ẍj + ω 2 j x j = p j de(ωt) d(ωt) , ( 30
)
where p j = e T j p. Making use of solution [START_REF] Witham | Linear and Non-Linear Waves[END_REF] for each of the oscillators [START_REF] Andrianov | Asymptotic solutions for nonlinear systems with high degrees of nonlinearity[END_REF] and taking into account expansion [START_REF] Boettcher | Nonperturbative square-well approximation to a quantum theory[END_REF], one obtains the steady-state response

y = N j =1 (e T j p)e j ωω j sin ω 0 ω j τ (ω j t) cos ω 0 ω j . ( 31
)
The resonance condition is

ω = 2 π ω j n ; n = 1, 3, 5, . . . ; j = 1, . . . , N.
These resonances can also be obtained in a regular way, by taking the generalized Fourier expansion of the right-hand side of Equation [START_REF] Andrianov | Asymptotic solutions for nonlinear systems with high degrees of nonlinearity[END_REF]. Note, however, that a convergence of the related Fourier expansions is slow, whereas solutions [START_REF] Witham | Linear and Non-Linear Waves[END_REF] and [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF] are not series and, hence, the convergence problem is of no importance for them.

NONLINEARIZABLE OSCILLATOR UNDER PERIODIC IMPULSIVE EXCITATION

Let us consider an essentially nonlinear oscillator with a power homogeneous characteristic under the periodic sequence of impulses

ẍ + x m = p de(ωt) d(ωt) , ( 32
)
where m is a positive odd number; p, ω = const. When m = 1, one has the impulsively forced harmonic oscillator considered in Section 3.2. For a large exponent, m 1, the system gives a smooth approximation of the simplest vibroimpact oscillator (Figure 1). Free vibration of the system admits an exact analytical solution in terms of complicated-enough special functions [START_REF] Lyapunov | Collection of Works[END_REF][START_REF] Rosenberg | The Ateb(h)-functions and their properties[END_REF][START_REF] Kamenkov | Stability and Vibrations of Nonlinear Systems[END_REF]. The oscillator was of interest to many authors [START_REF] Vakakis | Normal Modes and Localization in Non-linear Systems[END_REF][START_REF] Atkinson | On the superposition method for determining frequencies of nonlinear systems[END_REF][START_REF] Mickens | Construction of approximate analytical solutions to a new class of nonlinear oscillator equation[END_REF][START_REF] Manevitch | Method of Normal Vibrations for Essentially Non-Linear Systems[END_REF][START_REF] Boettcher | Nonperturbative square-well approximation to a quantum theory[END_REF][START_REF] Andrianov | Asymptotic solutions for nonlinear systems with high degrees of nonlinearity[END_REF]]. An approximate analytical solution for arbitrary exponent m was obtained in [START_REF] Pilipchuk | The calculation of strongly nonlinear systems close to vibro-impact systems[END_REF] in the power series form with respect to the saw-tooth function of time.

We will seek a periodic (T = 4/ω) solution of Equation ( 32) in the form

x = X(τ (ωt)), ( 33
)
where the Y -component of the representation is supposed to be zero, based on the final result of the previous subsection. Note that the equation considered is nonlinear and the two components of the transformed system are coupled. However, the system still admits a family of solutions on which Y ≡ 0. Substituting Equation (33) into Equation ( 32), one obtains

(ω 2 X + X m ) + (ω 2 X -p) de(ωt) d(ωt) = 0 (34)
and the equivalent boundary-value problem

ω 2 X + X m = 0, ( 35
)
ω 2 X | τ =±1 = p, ( 36
)
where the boundary condition stays for the elimination of the periodic singular term from the transformed Equation (34). In the case of free vibration (p = 0), the boundary condition becomes homogeneous and the parameter ω is treated as oscillator's natural frequency which was previously unknown and should be coupled with a parameter of amplitude.

Due to the symmetry provided by the oddness of m, the system admits a family of solutions where X(-τ ) ≡ -X(τ ). Taking into account this remark, we will seek an odd with respect to τ solution and take care of only one point, τ = 1, of the boundary condition (36).

In this example, the difference between free and forced vibrations is reduced to different values of the right-hand side of the boundary condition (36). A scheme of successive approximation [START_REF] Pilipchuk | The calculation of strongly nonlinear systems close to vibro-impact systems[END_REF] gives an approximate solution of Equation (35) in the power series form

X = A τ - ητ m+2 (m + 2) + mη 2 τ 2m+3 2(m + 2)(2m + 3) - m(m 2 + m -1)η 3 τ 3m+4 2(m + 2) 2 (2m + 3)(3m + 4) + • • • ,( 37
)
where τ = τ (ωt) and A is arbitrary constant, and

η = A m-1 (m + 1)ω 2 .
The solution is in power series form. However, it is periodic with respect to the original parameter of time, t.

The substitution of solution (37) into the boundary condition (36) for τ = 1 gives the algebraic equation

A 1 -η + mη 2 2(m + 2) - m(m 2 + m -1)η 3 2(m + 2) 2 (2m + 3) + • • • = p ω 2 . ( 38
)
Equation ( 38) remains for the elimination of the periodic singular terms4 from the transformed system (34).

If p = 0, Equation (38) couples the parameters A and ω associated with the amplitude and frequency of the free vibration. When p = 0, it gives a nonlinear frequency response for the system under external impulsive excitation (Figure 7).

A point of interest is the limit of the vibro-impact oscillator when m -→ ∞ (see Section 1). For a sufficiently large exponent m = 9, the phase trajectory's projection on the plane xv and the time history record for the velocity show the discontinuities of the velocity v created by external impulses when the coordinate x reaches its amplitude values (see Figures 8 and9). In the vibro-impact limit, when m -→ ∞ the amplitude values of the coordinate are defined by the rigid barriers appearing at points x = ±1. It means that both external and internal impacts must act at the same time instances to provide the periodic regime considered.

To this end, note that starting from the boundary-value problem in the form ω 2 X +f (X) = 0; ω 2 X | τ =±1 = p and using successive approximations, one can consider a more general case of the odd characteristic f (x) represented by Maclaurin's series.

Parametric Excitation

In this section, two cases of parametric excitation will be considered based on a relatively simple linear model. The time shapes of a piecewise-constant excitation and an impulsive one will be described by means of the functions e(ωt; γ ) and ∂e(ωt; γ )/∂(ωt), respectively. Both systems can be considered as particular cases of the more general example represented in [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF]. However, the parameters selected in [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF] correspond, in principal, to the case of impulsive excitation only (see Section 4.2). Note that both of the systems are piecewise linear and can be solved in a regular manner [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF]. These examples have been chosen to illustrate the applicability of the saw-tooth transformation of time to systems under parametric discontinuous excitations. Finally, there are two reason to use the transformation as a preliminary step in analysis. Firstly, the saw-tooth transformation automatically gives the conditions of matching at the points of discontinuity of the excitation. Secondly, the solution is represented as a closed-form expression, and this is an important feature in many cases.

PIECEWISE-CONSTANT EXCITATION

Let us consider a linear oscillator under the periodic piecewise-constant excitation

ẍ(t) + ω 2 0 [1 + εe(ωt; γ )]x = 0, (39)
where ω 0 , ω, γ are constants, and the parameter ε is not necessarily small. We will seek a periodic solution with the period of excitation T = 4/ω in the form [START_REF] Pilipchuk | On one variant of nonsmooth transformations of variables for 1-D elastic systems of a periodic structure[END_REF]. The coordinate of the system x(t) and velocity ẋ(t) must be continuous functions of time under the excitation considered. It means that neither velocity ẋ(t) nor acceleration ẍ(t) can include any δ-functions. The first derivative of representation [START_REF] Pilipchuk | On one variant of nonsmooth transformations of variables for 1-D elastic systems of a periodic structure[END_REF], due to expression [START_REF] Faure | Percussions en mécanique non linéaire sur certaines solutions périodiques de phénomènes non linéaires excités par des percussions[END_REF], is

ẋ(t) = αY + (X + βY )e + Y ∂e(ωt; γ ) ∂(ωt) ω. ( 40
)
The last term consisting of δ-functions can be excluded by imposing the boundary condition for Y -component of the solution as

Y | τ =±1 = 0. ( 41
)
Under this condition, the second derivative is

ẍ(t) = ω 2 [α(X + βY )] + ω 2 [βX + (α + β 2 )Y ]e
+ ω 2 (X + βY ) ∂e(ωt; γ ) ∂(ωt) , ( 42
)
where the multiplication rule for basic element [START_REF] Faure | Percussions en mécanique non linéaire sur certaines solutions périodiques de phénomènes non linéaires excités par des percussions[END_REF] has been taken into account once again. Analogously to the previous step, one eliminates the underlined term by means of the conditions

(X + βY )| τ =±1 = 0. ( 43
)
Substituting expressions [START_REF] Pilipchuk | On one variant of nonsmooth transformations of variables for 1-D elastic systems of a periodic structure[END_REF] and (42) into equation of motion (39) and taking into account Equation [START_REF] Faure | Percussions en mécanique non linéaire sur certaines solutions périodiques de phénomènes non linéaires excités par des percussions[END_REF] gives an element of the algebraic structure

{• • •} + {• • •}e.
Equating each of the two components of the element to zero, one obtains the set of differential equations for X(τ ) and Y (τ ) in matrix form, as

α αβ β α + β 2 X Y + r 2 1 αε ε 1 + βε X Y = 0, ( 44
)
where r = ω 0 /ω. We put

X Y = B 1 μ exp(λτ ); B, μ, λ = const. ( 45
)
From the differential equations for X and Y (Equation ( 44)), subject to the condition B = 0, we obtain a characteristic equation which has the following two pairs of roots:

λ 2 = [-(1 -γ )ε -(1 -γ ) 2 ]r 2 ≡ ±k 2 , λ 2 = [(1 + γ )ε -(1 + γ ) 2 ]r 2 ≡ ±l 2 .
We now consider the case of negative pairs. This means that the coefficient of x in Equation (39) is positive at any time t.

The general solution of a system of equations subjected to this condition will have the form

X = B 1 sin kτ + B 2 cos kτ + B 3 sin lτ + B 4 cos lτ, Y = μ 1 (B 1 sin kτ + B 2 cos kτ ) + μ 2 (B 3 sin lτ + B 4 cos lτ ), ( 46
)
μ 1 = - 1 α αk 2 -r 2 βk 2 -εr 2 , μ 2 = - 1 α αl 2 -r 2 βl 2 -εr 2 .
In view of the homogeneity of the boundary conditions, the requirement that the solution should be nontrivial leads to the relation

[μ 1 (1 + βμ 2 )l cos k sin l -μ 2 (1 + βμ 1 )k cos l sin k] × [μ 1 (1 + βμ 2 )l cos l sin k -μ 2 (1 + βμ 1 )k cos k sin l] = 0.
Figure 10 shows a family of curves in which this expression is equal to zero and, hence, periodic solutions exist with related instability zones. Figure 11 shows different time history records when crossing a region covered by two neighboring curves. An increasing of the amplitude in Figures 11b and11e indicates the parametric instability of the system, when the point (r, ε) is located inside the instability regions shown in Figure 10.

PARAMETRIC IMPULSIVE EXCITATION

Let us consider the case of a system with a single degree of freedom under the parametric impulsive excitation of the time shape given by a second derivative of the asymmetric sawtooth function, i.e. let us write the equation of motion as (compare with Equation (39)) This case was considered in [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF] based on the saw-tooth transformation of time. It was shown that the periodic solutions of the period T = 4/ω exist when parameters of the system satisfy the relation

ẍ + ω 2 0 1 + ε ∂e(ωt; γ ) ∂(ωt) x = 0. ( 47
)
p 2 = 2r 2 (1 -γ 2 ) 2 sin 2 2r cos 4r -cos 4γ r , ( 48
)
where r = ω 0 /ω and p = εr 2 . The dependence of p on r for fixed γ has a branched zone-like structure characteristic of similar problems (see Figure 12, where, in view of the symmetry, only the upper half-plane is shown). A curious effect is observed when the parameter γ which characterizes the slope of the 'saw' is varied. The question concerns the transverse 'collapse' of the sequences of zones which are regularly situated across each of the s zones, where s is an integer which depends on the parameter γ (the form of this relationship will be established later). For instance, when γ = 1/5, each fifth zone is missing and, when γ = 1/2, each second zone is missing as Figures 6a and6b, respectively, show. An explanation of this phenomenon was given in [START_REF] Pilipchuk | Calculation of mechanical systems with pulsed excitation[END_REF].

On the Use of the Method of Averaging

The above technique is applicable for analyzing periodic solutions, which are only achieved under certain initial conditions. Solutions of the Cauchy problem for an arbitrary initial condition are, generally speaking, nonperiodic and can be obtained by involving the idea of averaging. An appropriate formalization is given by the two-variables expansion technique [START_REF] Nayfeh | Perturbation Methods[END_REF].

Let us consider the first-order differential equation

v + εv 3 = p de(t) dt (49)
under the initial condition

v| t =0 = v • . ( 50
)
This Cauchy problem describes the one-dimensional motion of a material point under a periodic series of impulses, nonlinear damping and no elastic forces. Note that any standard formulation of the idea of averaging is hard to apply directly to Equation (49), because the right-hand side of the equation is neither smooth nor even a continuous function of time.

To proceed correctly, one should specify the class of solutions to be defined, because the equation combines the nonlinearity and singular functions (see Section 1 and Appendix 1). A preliminary analysis shows that a role in the solutions can be played by functions with discontinuities of the first kind at points = {t : τ (t) = ±1}. Solutions of the problem (49), (50) are sought in the form

v = X(τ, t • ) + Y (τ, t • )e; τ = τ (t), e = e(t), (51)
where the notation t • = εt indicates the 'slow time' scale. Hence, the solution is assumed to be dependent on two time variables -fast oscillating time τ and slow time t • . Substituting Equation (51) into Equation (49), we obtain the expression

∂Y ∂τ + ε ∂X ∂t • + X 3 + 3XY 2 + ∂X ∂τ + ε ∂Y ∂t • + 3Y X 2 + Y 3 e + (Y -p) ė = 0.
This finally leads to the set of partial differential equations

∂X ∂τ = -ε ∂Y ∂t • + 3Y X 2 + Y 3 , ∂Y ∂τ = -ε ∂X ∂t • + X 3 + 3XY 2 (52)
with the boundary condition

Y | τ =±1 = p. ( 53
)
The asymptotic (ε → 0) solution of problem ( 52), (53) can be found using the standard procedure of the method of two-variables expansion in a power series form, as

X = ∞ i=0 ε i X i (τ, t • ), Y = ∞ i=0 ε i Y i (τ, t • ). ( 54
)
Substituting Equation (54) into Equations ( 52), (53) and matching coefficients of the same powers of ε, gives a zero-order problem

ε 0 : ∂X 0 ∂τ = 0, ∂Y 0 ∂τ = 0; Y 0 | τ =±1 = p.
The solution satisfying the boundary condition is

X 0 = A 0 (t • ), Y 0 = p, ( 55
)
where A 0 is an arbitrary function of the slow time.

In the first order of ε, taking into account Equation (55) gives

ε 1 : ∂X 1 ∂τ = -(3pA 2 0 + p 3 ), ∂Y 1 ∂τ = - dA 0 dt • + A 3 0 + 3A 0 p 2 ; Y 1 | τ =±1 = 0.
Integrating the equations gives

X 1 = -(3pA 2 0 + p 3 )τ + A 1 (t • ), Y 1 = - dA 0 dt • + A 3 0 + 3A 0 p 2 τ + C 1 (t • ). ( 56
)
Figure 13. The closed-form analytical approach for the velocity of a nonlinearly damped material point of the periodic impulsive excitation.

Boundary conditions for Y 1 give C 1 (t • ) ≡ 0 and the leading-order averaged differential equation for

A 0 (t • ) is dA 0 dt • + A 3 0 + 3A 0 p 2 = 0. (57)
So, in the first order of ε, one has Y 1 ≡ 0, and the function A 1 (t • ) in the expression for X 1 remains unknown. Equation (57) admits an elementary integration which gives a solution in the following explicit form

A 0 = - 1 3p 2 + const • exp(6p 2 t • ) -1/2 .
Finally, a leading-order asymptotic solution of the Cauchy problem can be written as (Figure 13)

v = ± - 1 3p 2 + 1 3p 2 + (x • -p) 2 exp(6p 2 εt) -1/2 + pe(t) + O(ε), ( 58
)
where the plus or minus signs are taken as v •p > 0 or v •p < 0, respectively. The first term in expression (58) is a transient response of the system, and it tends to zero when the time tends to infinity. It can be seen that the impulsive external excitation strongly affects the time scale of the transient response, because that timescale is multiplied by the impulse parameter squared, p 2 .

Weakly Nonlinear Resonances in the Duffing Oscillator under the Periodic Impulsive Force

Let us consider the Duffing oscillator under the periodic impulsive excitation

ẍ + ω 2 0 x = -ε x 3 -p de(ωt) d(ωt) (59)
in the neighborhood of linear resonance condition [START_REF] Lyapunov | Collection of Works[END_REF]. Introducing the detuning parameter σ into the differential equation (59) as

ω 2 0 = ( nπ 2 ω) 2 + εσ (60) gives ẍ + nπ 2 ω 2 x = -ε σ x + x 3 -p de(ωt) d(ωt) . ( 61
)
We shall seek the solution in the form (51), where τ = τ (ωt), e = e(ωt). Assuming the continuity condition for x(t),

Y | τ =±1 = 0 (62)
gives the first derivative of Equation ( 51) as

ẋ = ω ∂Y ∂τ + ε ∂X ∂t 0 + ω ∂X ∂τ + ε ∂Y ∂t 0 e. ( 63
)
Substituting Equations ( 51) and (63) into differential equation (61) gives

ω 2 ∂X ∂τ τ =±1 = εp, ( 64
)
ω 2 ∂ 2 X ∂τ 2 + 2εω ∂ 2 Y ∂τ ∂t 0 + ε 2 ∂ 2 X ∂t 02 + nπ 2 ω 2 X = -ε(σ X + X 3 + 3Y 2 X), ( 65
)
ω 2 ∂ 2 Y ∂τ 2 + 2εω ∂ 2 X ∂τ ∂t 0 + ε 2 ∂ 2 Y ∂t 02 + nπ 2 ω 2 Y = -ε(σ Y + Y 3 + 3X 2 Y ). ( 66
)
As a result of the nonsmooth transformation of the argument, the boundary-value problem (62), (64-66) does not include any singular functions, and one can correctly use the averaging technique.

We seek an asymptotic solution in a power series form as

X = X 0 (τ, t • ) + εX 1 (τ, t • ) + ε 2 X 2 (τ, t • ) + • • • , Y = Y 0 (τ, t • ) + εY 1 (τ, t • ) + ε 2 Y 2 (τ, t • ) + • • • .
Substituting this into Equations ( 62) and (64-66), and matching coefficients of the same degrees of ε, gives a series of boundary value problems with respect to the X and Y components of the solution.

In the first step, one has the homogenous boundary value problem

∂ 2 X 0 ∂τ 2 + nπ 2 2 X 0 = 0, ∂X 0 ∂τ τ =±1 = 0, ∂ 2 Y 0 ∂τ 2 + nπ 2 2 Y 0 = 0, Y 0 | τ =±1 = 0 (67)
and its nonzero solution is written as

X 0 = A 0 (t • ) sin nπ 2 τ, Y 0 = D 0 (t • ) cos nπ 2 τ, ( 68
)
where A 0 (t • ) and D 0 (t • ) are arbitrary functions. Collecting all terms of order ε gives the boundary-value problem as

ω 2 ∂ 2 X 1 ∂τ 2 + nπ 2 2 X 1 = -2ω ∂ 2 Y 0 ∂τ ∂t 0 + σ X 0 + X 3 0 + 3Y 2 0 X 0 , ω 2 ∂X 1 ∂τ τ =±1 = p, ω 2 ∂ 2 Y 1 ∂τ 2 + nπ 2 2 Y 1 = -2ω ∂ 2 X 0 ∂τ ∂t 0 + σ Y 0 + Y 3 0 + 3X 2 0 Y 0 , Y 1 | τ =±1 = 0. (69)
The right-hand sides of Equation ( 69) depend implicitly on the slow time throughout the still unknown functions A 0 (t • ), D 0 (t • ), and its derivatives A 0 (t • ) and D 0 (t • ).

The solution of the differential equations (69) can be represented as

X 1 = A 1 (t 0 ) sin nπ 2 τ + B 1 (t • ) cos nπ 2 τ + 1 ω D 0 + 1 nπ ω σ A 0 + 3 4 A 3 0 + 3 4 A 0 D 2 0 τ cos nπ 2 τ + 2 nπ ω D 0 + 1 nπ ω σ A 0 - 9 16 A 3 0 - 21 16 A 0 D 2 0 sin nπ 2 τ - 1 8n 2 π 2 ω 2 (A 3 0 -3A 0 D 2 0 ) sin 3nπ 2 τ, Y 1 = C 1 (t 0 ) sin nπ 2 τ + D 1 (t • ) cos nπ 2 τ + 1 ω A 0 - 1 nπ ω σ D 0 + 3 4 A 2 0 D 0 + 3 4 D 3 0 τ sin nπ 2 τ + 1 4n 2 π 2 ω 2 (D 3 0 -3D 0 A 2 0 ) sin nπ τ sin nπ 2 τ.
The boundary condition for X 1 (see remark at the end of Section 2) gives

B 1 (t • ) ≡ 0 and D 0 = 1 nπ ω 2(-1) (n-1)/2 p + σ A 0 + 3 4 A 0 (A 2 0 + D 2 0 ) . ( 70
)
The boundary condition for Y 1 gives

C 1 (t • ) ≡ 0 and A 0 = - 1 nπ ω σ D 0 + 3 4 D 0 (A 2 0 + D 2 0 ) . ( 71
)
In this step, functions A 1 (t • ) and D 1 (t • ) remain unknown. Equations ( 70) and (71) allow us to determine the slow functions A 0 (t • ) and D 0 (t • ) and to complete the previous step. These equations can be investigated in a regular manner for the nonlinear vibration theory, introducing polar coordinates.

In the case when no external impulses are acting, i.e. p = 0, Equations ( 70) and ( 71) stay for elimination of the so-called 'secular terms' in the solution. If p = 0, secular terms occur in the solution. These terms, however, are only secular with respect to 'oscillating time' τ and, hence, the solution remains periodic with respect to original time, t.

Conclusion

It has been shown that the nonsmooth transformation of the argument is applicable to linear and nonlinear dynamical systems under periodic impulsive and discontinuous excitations. A standard pair of nonsmooth periodic functions plays an important role in the proposed method. In the symmetric case, γ = 0, the pair of functions which is used {τ (t), e(t)} is associated in a natural way with the motion of a free point mass between two fixed arresting devices, that is, with one of the two simplest mechanical oscillators. These systems are a vibro-impact system with a single impact pair and a harmonic oscillator. It is known that the latter generates a pair of trigonometric functions {sin t, cos t} possessing a number of convenient mathematical properties. From the mechanical point of view, the functions {sin t, cos t} are associated with a free rotating disc with angular velocity ω = 1. A 'complicated' linear regime of an elastic structure (say, a free mass-spring model of N degrees of freedom) usually involves a set of harmonics {sin ω k t, cos ω k t}, k = 1, . . . , N. It means that the regime can be represented by a set of N free rotating discs with angular velocities ω k , i.e. it appears that the oscillation of an elastic structure is equivalent to the rigid-bodies motion.

In the same manner, the functions {τ (t), e(t)} are associated to the motion of a free rigid body with absolutely rigid constraints. It has been shown that a more complicated nonsmooth (or a quasi nonsmooth) dynamical regime can be represented by this rigid-body motion and its different deformations, according to the power forms, τ 3 , τ 5 , . . . , or in another simple way. Hence, the introduction the pair {τ (t), e(t)} is not a purely mathematical manipulation. It has a certain physical meaning.

Another essential physical feature of the transformation is that the oscillating time, τ , is introduced in such a manner that it does not affect the metrical properties of the time, namely, the metric remains constant after the transformation: [dτ (t)] 2 = dt 2 . This defines a special role of the saw-tooth transformation among other periodic transformations of time. Due to this property, the transformed differential equations of motion remain within the framework of Newton's formulation (see Appendix 2).

Acknowledgement

This work was supported in part by the National Science Foundation of the U.S.A. Grant No. CMS-9634223.

Appendix 1 Carathéodory Equation

We will present some remarks concerning nonlinear differential equations with distributions involved as summands. We will also compare the two different representation for solutions based on nonsmooth functions.

Let us consider first-order (nonlinear) differential equation, ẋ = f (t, x), where x is a scalar or a vector, and the right-hand side satisfies the Carathéodory conditions [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]. Namely, in the domain D of the (t, x) space, the function f (t, x) is defined and continuous in x for almost all t; the function f (t, x) is measurable in t for each x; |f (t, x)| ≤ m(t), and the function m(t) is summable (at each finite interval if t is not bounded in D).

The conditions listed above mean that the restriction taken for the right-hand side's dependence on t is less strong than the classical existence theorem requires, i.e. the function f (t, x) is not necessarily continuous in t. It becomes possible if the right-hand side of the equivalent integral equation, x(t)x(t 0 ) = t t 0 f (ξ, x(ξ)) dξ is understood in the sense of Lebesque. Under the Carathéodory conditions, a series of existence theorems holds [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF].

Let the right-hand side of the equation additively involve a δ-impulse and the equation is written as

ẋ = f (t, x) + νδ(t -a), ( 72
)
where a and v are constant parameters (ν is a scalar or a vector). Many mechanical and physical problems can be represented in the form (72). EXAMPLE 1. A parametrically excited Duffing oscillator subjected to external impact at moment t = a is described with the equation

ü + ω 2 0 (1 + ε cos 2t)u -βu 3 = 2pδ(t -a). ( 73
)
Equation ( 73) is equivalent to Equation (72) when

x = u v , f (t, x) = v βu 3 -ω 2 0 (1 + ε cos 2t)u , ν = 0 2p .
In this case, the variable v is the velocity. EXAMPLE 2. After a special scaling of the variables, equation

d 4 W dξ 4 + γ (ξ)W = qδ(t -a) ( 74
)
describes a beam's centerline, W = W (ξ ), resting on an elastic foundation of variable stiffness, γ (ξ) and loaded by a transverse force localized at point ξ = a. Equation (74) can be represented in the form of (72) when t ≡ ξ ,

x = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ W dW dξ d 2 W dξ 2 d 3 W dξ 3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , f (ξ, x) = ⎡ ⎢ ⎢ ⎢ ⎣ 0 1 0 0 0 0 1 0 0 0 0 1 -γ (ξ) 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ x and ν = ⎡ ⎢ ⎢ ⎢ ⎣ 0 0 0 q ⎤ ⎥ ⎥ ⎥ ⎦ .
The right-hand side of Equation (72) no longer satisfies the Carathéodory conditions. However, a change of variables

x(t) -→ y(t) : x(t) = y(t) + vH (t -a) ( 75
)
brings the differential equation back into the form

ẏ = f (t, y + vH (t -a)) ≡ f (t, y), ( 76
)
where the Dirac function has been eliminated due to equality Ḣ (ta) = δ(ta) and, as a result, the right-hand side, f (t, y), satisfies the Carathéodory conditions. Note that due to H (ta), the right-hand side of the transformed equation is generally piecewise continuous in t, even if the function f (t, x) at the right-hand side of the original equation ( 72 (77)

We also introduce the old notation for the first-time derivative of the nonsmooth argument as ṡ(t) = d|t -a|/d(ta) = e(t). One has e(t) = -1 for ta < 0 and e(t) = 1 for ta > 0 and hence, e 2 = 1 at least for almost all t. where the expressions

R f = 1 2 [f (a + s, X + Y ) + f (a -s, X -Y )], I f = 1 2 [f (a + s, X + Y ) -f (a -s, X -Y )]
have been obtained analogously to Equation (80). Equality (81) gives a set of equations and a boundary condition

X + I f (s, X, Y ) = 0, Y + R f (s, X, Y ) = 0, Y | s=0 = v 2 . ( 82
)
As a result of the transformation, the Dirac function has been eliminated from the system. Equations (82) compared with Equation (76) have a more complicated form. However, one should pay attention, for its smoothness properties may become better. In fact, the step function H (ta) has not been involved in Equation (82). Another (sometimes suitable) point is that the new time variable s is half-limited, 0 ≤ s < ∞.

Appendix 2

Let us show that the saw-tooth time τ (t) is associated with the group properties of the oscillator d 2 x dt 2 + f (x) = 0 (83) and, hence, plays a special role among all possible periodic transformations of time.

In this particular case, a one-component version of the saw-tooth transformation of time can be realized. Alternatively, one may introduce an oscillating time parameter, say g, in a more general manner as g = g(t) and x = X(g(t)), where g(t) is a periodic function, not necessarily saw-tooth. One could consider trigonometric sine, elliptic sine, and other specially selected functions depending on the problem investigated. As a result, one obtains a new system in the form dg dt This transformation becomes physically reasonable if the variable coefficients appear to be functions of the new time parameter g over all periods of the oscillation. For example, let g = sin t. In this case, the transformation t -→ g, x -→ X is completed as

(1 -g 2 ) d 2 X dg 2 -g dX dg + f (X) = 0. ( 85
)
In this example, one undesirable result of the transformation is that the system has lost its original Newton form given by Equation (83) and, hence, the dynamical properties of the system cannot be treated any more in a regular manner.

To preserve Newton's form of the system, one should set dg dt in Equation (84). The first equality is satisfied by g = α + t and g = βt, where α and β are arbitrary constants. It means that the function g(t), to be periodic, must consist of different pieces of linear functions with positive and negative slopes of the same absolute value, i.e. it has a piecewise linear saw-tooth form. The simplest (monofrequency) case of the saw-tooth function is given by τ (t) (Equation ( 3)), as Figure 3 shows. In this case, the second equality of Equation ( 86) holds for all t under the condition (dX/dτ )| τ =±1 = 0, and the new equation has the same form as the original one: (d 2 X/dτ 2 ) + f (X) = 0.

Figure 1 .

 1 Figure 1. Geometrical treatment of the nonsmooth transformation of a positional variable on the simplest vibro-impact oscillator. In this case, x ∈ [-1, 1] is the original coordinate of the system, while s ∈ (-∞, ∞) is the distance passed by the system. The figure illustrates reflection from the right barrier.

Figure 2 .

 2 Figure 2. Geometrical treatment of the nonsmooth transformation of time. The new time oscillates in the bounded domain τ ∈ [-1, 1] with the constant velocity reflecting at points τ = ±1, while the original time is unbounded, t ∈ (-∞, ∞).

Figure 3 .

 3 Figure 3. The saw-tooth sine and its first two generalized derivatives.

Figure 4 .

 4 Figure 4. The saw-tooth function with asymmetric slope and its first two generalized derivatives. Parameter γ defines the teeth's inclination; it has been taken as γ = 0.7.

Figure 5 .

 5 Figure 5. Family of periodic solutions of Equation (11) for μ = 0.2 and λ = 0.4, 0.7, 1.0, 1.3, 1.6, 1.9.

Figure 6 .

 6 Figure 6. Different responses of the damped harmonic oscillator under the periodic impulsive excitation for p = 0.1, ζ = 0.5, ω 0 = 4 and (a) ω = 0.2 (large distance between the impulses); (b) ω = (2/π)ω 0 = 2.5465 (principal resonance).

Figure 6 .

 6 Figure 6. Continued. (c) ω = 6 (small distance between the impulses).

Figure 7 .

 7 Figure 7. Frequency response of the 'quasi-impact' oscillator under periodic impulsive excitation for p = 2, m = 3 (dashed line) and m = 9 (solid line). In standard normalization of the linear vibration theory, the frequency is (π/2)ω. The neighborhood of zero frequency is not shown.

Figure 8 .

 8 Figure 8. Projection of the impulsively forced 'quasi-impact' oscillator's orbit on the plane 'coordinate-velocity', xv. The velocity experiences momentary steps (see also the next figure), when the system's coordinate reaches its amplitude values. The periodic regime corresponds to the parameters p = 2, m = 9, ω = 3 and A = 1.79354.

Figure 9 .

 9 Figure 9. The velocity time record shows sudden vertical 'steps' caused by external impulses.

Figure 10 .

 10 Figure 10. Instability zones of the linear oscillator under the parametric piecewise constant force described by the first derivative of the asymmetric saw-tooth function for the parameter γ = 0.7. A numerical testing of the self-intersections of the instability zone boundaries is presented by the series of plots in Figure 11.

Figure 11 .

 11 Figure 11. Different time history records when taking the system's parameters gradually around the point of self-intersection of the instability zone (point P in Figure 10). (a) r = 3.20, ε = 0.42 (on the left of the instability zone, under point P ), (b) r = 3.24, ε = 0.42 (inside the instability zone, under point P ).

Figure 11 .

 11 Figure 11. Continued. (c) r = 3.28, ε = 0.42 (on the right of the instability zone, under point P ), (d) r = 3.24, ε = 0.52 (on the left of the instability zone, above point P ).

Figure 11 .

 11 Figure 11. Continued. (e) r = 3.28, ε = 0.52 (inside the instability zone, above point P ), (f) r = 3.32, ε = 0.52 (on the right of the instability zone, above point P ).

Figure 12 .

 12 Figure 12. 'Collapse' of the instability zones of the harmonic oscillator under the periodic impulsive excitation with the dipole-like shift: (a) γ = 1/5, each fifth zone is missing, (b) γ = 1/2, each second zone is missing.

  ) is continuous. 'Single-Tooth' Transformation of Time Let us transform Equation (72) by making use of a nonsmooth argument (a 'single-tooth' version of the saw-tooth transformation of time). Let the new argument be |t -a| = s = s(t).

LEMMA 1 .

 1 An arbitrary function x(t) can be represented in formx(t) = X(s) + Y (s)e. (78)Proof. A transformation inverted to (77) is expressed ast = a + se. (79)Substituting Equation (79) into Equation (78) and taking e = -1 and e = 1, one obtainsx(as) = X(s) -Y (s), x(a + s) = X(s) + Y (s) a + s) + x(as)], Y (s) = 1 2 [x(a + s)x(as)]. (80)So when X(s) and Y (s) are selected as Equation (80), the identity (78) holds for any x(t).Let a solution of Equation (72) be represented in the form of (78). Substituting Equations (78) and (79) into Equation (72) and taking into account the expression ė(t) = 2δ(t -a), one obtainsY + R f (s, X, Y ) + [X + I f (s, X, Y )]e + (2Yv)δ(ta) = 0,(81)

Complex numbers X + Y e, where e

= 1 (see the next section).[START_REF] Richtmyer | Principles of Advanced Mathematical Physics[END_REF] In terms of the related field of the applied mechanics, see V. I. Babitsky, Theory of Vibroimpact Systems, Moscow, Nauka, 1978 [in Russian].

Note that elements 1 and e are linearly independent on the time domain considered.

Recall that the quasilinear techniques give the corresponding relations after the elimination of the secular terms.


T M Cameron

J H Griffin

Mem Asme

An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems

A method is proposed for analyzing the steady-state response of nonlinear dynamic systems. The method iterates to obtain the discrete Fourier transform of the system response, returning to the time domain at each iteration to take advantage of the ease in evaluating nonlinearities there-rather than analytically describing the nonlinear terms in the frequency domain. The updated estimates of the nonlinear terms are transformed back into the frequency domain in order to continue iterating on the frequency spectrum of the steady-state response. The method is demonstrated by solving a problem with friction damping in which the excitation has multiple discrete frequencies.

Introduction

A variety of methods are presently available for calculating the response of nonlinear dynamic systems, (see, for example, Iwan (1974) for analytical approaches, and [START_REF] Press | Numerical Recipes[END_REF] for numerical approaches). The method introduced in this paper offers advantages over most other methods in its generality, efficiency, accuracy, and ease-of-use for calculating the steady-state response of forced nonlinear systems. The method involves an iterative scheme that alternates from the frequency to the time domain in order to take advantage of the convenience in performing certain tasks in each domain. Consequently, this Alternating Frequency /Time (AFT) method easily provides either time trajectory or frequency spectrum solutions.

The purpose of this paper is to introduce the AFT method and illustrate its utility with a friction damping problem: the response of a single degree-of-freedom (sdf) system with an elastic/perfectly-plastic friction damper subject to three discrete frequencies of excitation is calculated.

Alternative Spectral Methods

One of the most popular methods for approximating the frequency response of nonlinear systems is known as the Harmonic Balance (HB) method in the vibration literature (e.g., [START_REF] Nayfeh | Nonlinear Oscillations[END_REF], and as the Describing Function (DF) approach in the control literature (e.g., Gelb and Vander Velde, 1968). This is essentially a spectral method in which the solution is assumed to have one dominant frequency component-or at most two. 1 In situations where it can be applied, and trusted to give accurate solutions (where a single harmonic description of the nonlinearity is adequate), the HB method is by far the most computationally efficient method for obtaining steady-state solutions to nonlinear dynamics problems. It has been used extensively to study friction damping (e.g., [START_REF] Griffin | Friction Damping of Resonant Stresses in Gas Turbine Engine Airfoils[END_REF]Sinha and[START_REF] Sinha | Friction Damping of Flutter in Gas Turbine Engine Airfoils[END_REF][START_REF] Sinha | Effects of Static Friction on the Forced Response of Frictionally Damped Turbine Blades[END_REF][START_REF] Dowell | The Behavior of a Linear, Damped Modal System with a Nonlinear Spring-Mass-Dry Friction Damper System Attached[END_REF][START_REF] Jezequel | Structural Damping by Slip in Joints[END_REF]. Generally, however, this method is limited to harmonically excited systems and it becomes increasingly inaccurate as higher harmonics in the expansion of the nonlinear force become significant. A fair amount of work is also required to expand the nonlinear force, although DF's for a large number of nonlinearities have been tabulated (Gelb and VanderVelde, 1968, Appendix B).

An Incremental Harmonic Balance (IHB) method which accommodates multiple harmonic components in the description of the nonlinear force was introduced by [START_REF] Lau | A Variable Parameter Incremental Method for Dynamic Instability of Linear and Nonlinear Elastic Systems[END_REF] The method was originally applied to continuous (cubic) nonlinearities [START_REF] Lau | A Variable Parameter Incremental Method for Dynamic Instability of Linear and Nonlinear Elastic Systems[END_REF][START_REF] Lau | Incremental Harmonic Balance Method with Multiple Time Scales for Aperiodic Vibration of Nonlinear Systems[END_REF], but [START_REF] Pierre | Multi-Harmonic Analysis of Dry Friction Damped Systems Using an Incremental Harmonic Balance Method[END_REF] extended the method to analyze the rigid/perfectly-plastic friction damper problem in which the nonlinear force is discontinuous. The IHB method appears to handle strong nonlinearities well, and provides accurate solutions over a much broader range than the HB method since it accommodates multiple harmonics. The method is ideally suited to parametric studies since it can step from a known state of vibration to a neighboring state which corresponds to an incremental change in one of the governing parameters of the ----r:rhe DF approach has been extended to allow two input frequencies in the case of odd, memoryless (nonhysteretic) nonlinearities [START_REF] Gelb | Multiple-Input Describing Functions and Nonlinear System Design[END_REF], Chapter 5 and Appendix D), but the nonlinear force and the system response are described only in terms of the two input frequencies.

system. However, there are also significant disadvantages to the IHB method.

Computationally the IHB method is no more efficient than numerical integration, at least for the discontinuous friction damping problem (cf., Pierre et a!., 1985, "Conclusions"), and the computational effort grows geometrically with the number of harmonics considered. As currently implemented, the frequencies are integer multiples of the excitation frequency, which will hide the subharmonic resonances that may occur in nonlinear systems. Another disadvantage of the IHB method is the amount of analytical work required. Functions in the equation are expanded in a first order Taylor series (which may be difficult to obtain for discontinuous functions) and substituted into the original equation to obtain a linear "incrementalized" equation. The incrementalized equation is solved by the Galerkin method, which requires further expansions and integrations. Finally, although the approach is general, the expansions are specialized, so each new dynamic system requires its own formulation and implementation.

Other numerical methods, using techniques similar to those in the AFT method, describe functions in both the frequency and time domains. In computational fluid dynamics, "Spectral'' and "Pseudo-Spectral" methods (cf., Zang and Hussaini, 1985) use Fourier transform techniques to expand functions into (trigonometric or Chebyshev) polynomials. This permits derivatives and other operations on those functions to be performed with higher order accuracy than finite difference methods achieve for similar effort. [START_REF] Cesari | Functional Analysis and Periodic Solutions of Nonlinear Differential Equations[END_REF] appears to be the first to obtain steady-state solutions to nonlinear dynamic equations by using a trigonometric polynomial expansion, which is equivalent to a truncated Fourier series or a discrete Fourier transform. [START_REF] Urabe | Galerkin's Procedure for Nonlinear Periodic Systems[END_REF] expands Cesari's results and [START_REF] Urabe | Numerical Computation of Nonlinear Forced Oscillations by Galerkin's Procedure[END_REF] shows how the Fourier coefficients may be computed at each iteration with a discrete Fourier transform, once the nonlinear functions have been evaluated in the time domain. Urabe identifies this approach as an application of the Galerkin method. [START_REF] Ling | Fast Galerkin Method and Its Application to Determine Periodic Solutions of Non-Linear Oscillators[END_REF] use a fast Fourier transform (FFT) to implement Urabe's approach, calling it the "Fast Galerkin" (FG) method. Although there are fundamental similarities in the concepts underlying the FG and AFT methods, there are differences in how the methods are formulated and implemented. These differences are discussed later.

3 An Alternating Frequency /Time (AFT) Method 3.1 Fundamental Concepts. Frequency domain methods have long been used to analyze linear systems because the frequency content of a signal is often of greater interest than the time history and because convolution-a convenient solution method-is more efficiently accomplished in the frequency domain . In general, nonlinear forces are nonlinear functions of the past and present states of the system. Since they are not linear functions of the states or explicit functions of time, they cannot be transformed directly into the frequency domain. This makes nonlinear systems difficult to analyze in the frequency domain, as evidenced by the complicated expansions required in the frequency domain methods discussed previously. Nonlinear forces are easy to evaluate in the time domain, however, since the state history of the system is readily available. This is evidenced by the ease of implementation and greater complexity allowed by numerical integration. Th. e AFT method employs the advantages of both domains by performing convolution (as well as filtering and convergence testing) in the frequency domain and by switching to the time domain to evaluate the nonlinear forces as functions of the state trajectory. The cost involved in this technique is the time required to switch back and forth between the frequency arid time domains. A fast Fourier transform 2 algorithm makes this cost reasonable, especially in obtaining steady-state solutions to periodically excited systems.

Consider a nonlinear dynamic system in which the terms containing nonlinearities are grouped on the right hand side of the equation with the external forcing term. The equation of motion for this system may be represented by:

L(x(t)J =f(t) +Nj(X,X) (1)
where L is the differential operator representing the linear terms, f(t) is the external forcing term, Nj( •) represents the nonlinear force terms, and x(t) is the desired system response. These are vector functions for multi-degree-of-freedom (mdf) systems, or scalar functions for sdf systems. Taking the Fourier transform (FT) of both sides yields:

A(w)•X(w)=F(w)+1j(w,X(w)) (2)
where A(w), the transform of L, is a (matrix) function of the frequency, w, F(w) is the transform of f(t), 1j(w,X(w)) is the transform of Nj(•), and X(w) is the transform of x(t). In the transform domain the nonlinear forces depend only on w and X(w) because derivatives of x(t) transform into linearly scaled forms of X(w):3

x(t) #X(w)

x(t)#-jwX(w) x(t) # -w 2 X(w)
(" #" denotes a Fourier Transform pair).

(3)

Since an FFT will be used to compute and represent the frequency dependent functions, ( 2) is required to hold at the discrete frequencies, wk, introduced by the FFT. Thus, (2) becomes a system of nonlinear algebraic equations4 : (6)

A(w)•X(w) =F(w) + 1j (w,X(w)) (4)
Here, AT is the sampling period and N is the number of samples in the FFT of each time domain function. In general, A is a matrix whose components are functions of the independent variable (wk), while X, F, and 7l are vector functions of wk. For example, at a given frequency, say w 4 , there will be X 1 (w 4 ), ••• , Xm(w 4 ) if there are m degrees-of-freedom. Therefore, the length of X(w) is the product of the number of frequencies in the FFT (N/2+ 1) times the number of degreesof-freedom (m). To minimize the bandwidth of L, the equations should be ordered so that X(w) appears as {X 1 (w 0 ),

X2(wo), . .. 'Xm (wo), Xl(wl), . . . 'Xm (wl), ... 'Xl(wNn), • • • Xm(wNn )J.
Various iterative methods are available for solving nonlinear algebraic equations, like (4). One problem facing any method is evaluation of the nonlinear term 1j(w, X (w)). To add it directly to F (w) this term must be known explicitly in terms of w, as 7l (w), but this form is difficult to determine in ~Fast Fourier transform" refers to a class of algorithms used to calculate the discrete Fourier transform (OFT-cf., Jackson, 1986, Chapter 7;and Oppenheim and Schafer, 1975, Chapter 3). The notation "FFT" will be used here to refer to the OFT as well as its implementation. The authors use an FFT subroutine from [START_REF] Press | Numerical Recipes[END_REF] . the frequency domain. However, if N 1 ( o) can be expressed as an explicit function of time, N 1 ( t), then it can be transformed to provide 77(w). To find N 1 (t), the most recent estimate of X(w) is inverse transformed to provide an estimate of x(t). Equation (3)-or the analogous FFT relationships-may be used to obtain other state variables so that N/ o) may be evaluated as N 1 (t), and transformed to 77(w). For nonlinearities that involve products of time domain functions this avoids frequency domain convolution. For historydependent nonlinearities, this procedure allows the nonlinearity to be evaluated directly from the state history as long as the time period used in the analysis is longer than the ''term of memory" of the nonlinearity. With short term memory nonlinearities, values of the nonlinear terms at the beginning of the period may be obtained from periodicity requirements. Another problem for nonlinear equation solvers is obtaining an initial estimate of the solution. An initial estimate (designated .XO(w)) may be obtained from the linear system, ignoring 77( o ), or from a known solution corresponding to a closely related state. For example, consider a sdf system with f(t) =A cos(w.t) where a parametric study of the effect of A on the system response is desired. Once the solution, X(w), is obtained for one value of A, it may be used as the initial estimate for a neighboring case with a slightly different value of A.

3.2 Implementation Considerations. The authors presently use two methods to solve (4). The first method, direct iteration, is an iterative implementation of (4):

A(w)o_x{O(w) =F(w) + 17(w,XU-1 (w)). ( 7) A(w) 5 and F(w) are fixed by the system and the external forcing function, respectively, once the w/s are selected. New generations of XUl (w) may be found by evaluating 77( • ), as described previously, adding it to F(w), and multiplying by A -1 • Since A is a constant matrix representing the linear part of the system, it may be inverted once and saved. A flow chart of the direct iteration approach is shown in Fig. 1.

The second approach used by the authors for solving (4) is the Newton-Raphson method (cf., [START_REF] Stoer | Introduction to Numerical Analysis[END_REF], with a difference quotient approximation of the Jacobian. For the example problem in this paper, Newton-Raphson is more robust than direct ithation (it converges at times that direct iteration fails) and it converges in fewer iterations. However, when both methods converge, computation time is greater for Newton-Raphson than for direct iteration because of the effort required to recompute the Jacobian at each iteration. Newton-Raphson is also more difficult to implement. In some cases Newton-Raphson fails to converge, so more robust methods (cf., [START_REF] Wolfe | The Secant Method for Simultaneous Nonlinear Equations[END_REF][START_REF] Broyden | A New Method of Solving Nonlinear Simultaneous Equations[END_REF][START_REF] Stoer | Introduction to Numerical Analysis[END_REF][START_REF] Shacham | Numerical Solution of Constrained Nonlinear Algebraic Equations[END_REF]) may be required if a solution is desired for all regions of the parameter space.

---sp;ofessor Alok Sinha of Pennsylvania State University has pointed out that in the absence of damping A may be singular. This may be remedied by adding a linear damping term to both sides of any equation at fault. The damping term may be lumped with the nonlinear term(s) on the right-hand side of an equation.

Equation (4) represents m(1 + N/2) equations (for m degrees-of-freedom) which for most methods is unwieldy, except for small values of N. The computational effort of the direct iteration method, O[mN(m + log(N)] multiplications per iteration, is much less affected by "N' than Newton-Raphson, which has O[(mN) 2 •log(N)] multiplications per iteration. In the approach used here, rather than limiting N, which would greatly reduce the number of frequencies considered in the solution, the direct iteration approach is used first to estimate which frequency components are significant. Then, if Newton-Raphson is employed, equation ( 4) is implemented only for those values of k for which X(') (wk) is a significant component in the last estimate of the solution vector obtained from the direct iteration method, (_x{O (w 0 ), • • • , _x(l) (wN 12)). This is discussed further in Section 5. This approach, by solving only the significant frequency components, keeps the number of equations and unknowns at a minimum while still allowing a full spectrum of frequencies to be considered in the solution.

There are two types of error that may corrupt FFT values. Aliasing error arises when the FT is nonzero for frequencies above the maximum frequency considered by the FFT (wmax = 1r I !J.T). Power at the higher frequencies is aliased or folded back and added to the FFT values for w < 1r/!J.T.

Leakage error occurs when the frequency values of the FFT (k•!J.w, or wk) do not coincide precisely with each frequency in the FT that has nonzero power. Power at a frequency not represented in the FFT is leaked to frequencies that are represented. Leakage can only be eliminated if the sampled time signal is periodic.

For the following sample problem aliasing is insignificant (as verified by taking smaller time steps), and both the excitation on the system and the system response are periodic, so no leakage will occur. Most nonlinear forces are periodic when subject to a periodic input, so the assumption of a periodic response to a periodic excitation is typically valid. Discussions of the DFT, FFT algorithms, aliasing, and leakage may be found in Jackson (1986), [START_REF] Oppenheim | [END_REF], and Press et a!. ( 1986).

Sample Problem

In many mechanical applications the most significant source of vibrational damping is provided by friction. For example, friction damping occurs naturally at joints in fabricated structures, friction damping is used to control the earthquake response of buildings, and friction dampers are used in high speed turbomachinery to reduce resonant stresses in blades. An elastic/perfectly-plastic friction model may also be used to analyze the elasto-plastic behavior of materials. Interest in these applications has motivated considerable research into how friction affects the dynamic response of structures (cf., Plunkett's review of friction damping research, 1980).

In the problem considered the nonlinear element consists of a spring in series with a friction constraint (Fig. 2). The friction contact is assumed to exhibit simple Coulomb-type behavior, i.e., when slip occurs the force is constant in magnitude and directed to oppose relative motion. Since the damper force in this model depends on the history of the response, it is easiest to illustrate the damper force by example. In Fig. 2 (b), for instance, the damper force exhibits bilinear hysteresis, which occurs when the input to the damper is harmonic. For more complicated motion the damper force could take on any value within a similar hysteresis loop, where the upper and lower limits of the input, which represent the extreme amplitudes of motion, need not be symmetric about zero. In the sample problem the equations have been normalized so that the friction force, 11-C 1 , is unity and does not appear explicitly. The system used to illustrate the AFT method has only one degree-of-freedom. From the formulation in Section 3.1 it is clear that the method easily extends to multiple degrees-offreedom (mdt). The computational savings of the AFT method over numerical integration should be even more dramatic for mdf systems since the number of frequencies to consider depends only on the excitation and the nonlinearities, not on the number of degrees-of-freedom.

An Elastic/Perfectly-Plastic Friction Damper.

The sample problem was previously unsolved, although it can be solved easily by numerical integration. It would be either difficult or impossible to solve by other methods. The accuracy of the solution is judged by comparison to numerical integration. The comparison is made in the time domain to show the accuracy of the AFT method in describing both the amplitude and phase of the steady-state response.

Multiple discrete frequencies of excitation arise in high speed turbomachinery applications because as the turbine blades rotate they pass through a nonuniform, circumferentially distributed pressure field. This introduces excitation frequencies that are integer multiples of the rotation speed of the turbomachine. It is also frequently necessary to model frictio•n dampers as having compliance, rather than being perfectly rigid. These physical considerations motivate the sdf model with an elastic/perfectly-plastic friction damper subjected to three discrete frequencies of excitation (Fig. 3). The nondimensionalized equation of motion for this system is: where f(t) =0.75 sin (0.5 1)+0.5 sin(0.75t)+0.375 sin(t), .1.

x+2tx+(1-E)x=f(t) -Ey(x,x) (8)
the fraction of critical viscous damping, is 0.01 (1 percent), and E, a measure of the relative spring stiffnesses, is 0.25. The distance the mass must displace for slip of the friction damper to occur is unity. The nonlinear force on the mass is equal to the stiffness, E, of the elastic friction damper times the displacement of the friction contact point, y( • ). The contact slips when the magnitude of the friction force equals E; the resulting force resists the motion. The specific excitation was chosen to provide an interesting response.

One period of the steady-state response of this system is shown in Fig. 4. The AFT and numerical integration solutions compare almost perfectly in both amplitude and phase. Since slip occurs at an amplitude of unity, there is considerable participation of the friction damper in this response. The frequency spectrum of the response (FFT) is shown in Fig. 5. This problem is used in the following section to make accuracy and efficiency comparisons between the AFT method and fourthorder Runge-Kutta integration.

Convergence, Error, and Computational Efficiency.

In the absence of aliasing and leakage the FFT of a function is equivalent to a truncated Fourier series representation. Therefore, the AFT method should exhibit the convergence and error characteristics of Fourier series. Given a periodic excitation, with period II, and assuming the time sampling is over one period, the maximum frequency in the FFT representation is: Since the highest frequency is proportional to the number of samples per period, N, a continuous function which exhibits w-2 convergence, like A(w)-1 for second order systems, will exhibit N-2 convergence in its FFT components. For the sample problem, the nonlinear force on the system, N 1 (t), is continuous, but has a discontinuous derivative, so 11(w) will also converge O(N-2 ). Therefore, X(w) is expected to converge O(N-4 ), which is confirmed by Fig. 5. (F(w) has infinite order convergence so it does not affect the convergence of X(w).)

In order to compare convergence of the time domain errors, as N increases, for the AFT and Runge-Kutta methods the error is assumed to be of the form:

EN=XN-Xex =aNc (10) where xN is the value of x( t) at a specified time using N sample points, Xex is an extrapolated "exact" answer, EN is the estimated error, N is the number of sample points per period, "a" is a constant, and "c" is the constant that gives the order of convergence of the error. Using three values of N gives three equations that allow the parameters Xex• a, and c to be evaluated. The consistency of the error convergence is checked with two of the original three values of N, plus and an additional, fourth value. These permit another calculation of the parameters in (10).

For the AFT method, the O(N-4 ) convergence of the FFT should produce similar convergence in the time domain error. This is the case if the FFT components are known precisely. However, changing N affects truncation and aliasing errors, so for different values of N the FFT components for the same frequency vary slightly. Consequently, additional errors appear in the time domain. Using N = 64, 128, and 256, in single precision, the AFT time domain error for equation ( 10) was found to be O(N-3 • 1 ). Since the sampling period, !::..T, is II/N, this is equivalent to O(!::..Tl• 1 ) convergence in the error. The error was consistent when estimated with N= 128, 256 and 512.

When the functioli.s in a problem are infinitely differentiable a fourth-order Runge-Kutta (R-K) integration has error O(!::..T'), or O(N-4 ). When a function has discontinuous derivatives, like Ey(t) (the nonlinear force in problem 4.1), that error estimate is invalid. Since the R-K method steps progressively forward in time, accumulated roundoff error can also dominate truncation error. In fact, for this problem, no consistent error convergence (using (10)) was observed for R-K in single precision, with N = 64, 128, 256, 512, and 1024. Using the same values of N in double precision R-K simulations, the error was O(!::..TL 6 ), or O(N-1. 6 ). The error results for both R-K and AFT methods are plotted in Fig. 6.

It is impossible to make a general comparison of computational efficiency between the AFT and R-K methods. The efficiencies of both depend strongly on system parameters.

Depending on damping, for example, the system may require any number of periods of oscillation before steady-state is reached. Another concern for R-K is that with complicated wave forms, like the solution to this problem (Fig. 4), it is very difficult to determine when steady-state has been achieved. The AFT method tests for convergence in the frequency domain where it is easy to discern. However, for a sample comparison, the following test was performed. This sample problem-a sdf system with three input frequencies and an elastic/perfectly-plastic friction damper-was analyzed through the 30th period of oscillation. The 30th period was selected as representative of an average number of cycles to reach steady-state. In order to achieve the same accuracy (i.e., similar error in Fig. 6) the AFT method was implemented in single precision with 256 points per period, and the R-K method was implemented in double precision with 512 points per period. In this case the AFT method took 21.7 seconds 6 and Runge-Kutta took 389 seconds on a Macintosh II, a factor of 18 difference.

Features of the AFT Method

Certain fundamental concepts underlying the AFT method are also employed in the FG (Fast Galerkin) method introduced by [START_REF] Ling | Fast Galerkin Method and Its Application to Determine Periodic Solutions of Non-Linear Oscillators[END_REF]. Ling and Wu implement Urabe's Galerkin method with an FFT, making the method computationally efficient. In Urabe's approach the system of nonlinear equations is put into state variable form, a solution is asssumed in the form of a truncated Fourier series, and the coefficients of the Fourier series are solved for iteratively-with nonlinear terms evaluated in the time domain. Since the coefficients of the Fourier series are equivalent to the components of the DFT, the AFT and FG methods are similar in this respect. There are, however, differences in how the methods formulate and implement solutions of equations. These differences highlight some of the features of AFT method.

The first difference between the AFT and Galer kin methods lies in the form of the equations to be solved. The AFT method permits the equations to be of any order, while Urabe's formulation, which Ling and Wu repeat, assumes the equations are in first order state variable form. The AFT method only finds the linearly independent functions, thereby reducing the number of equations to be solved.

Another feature of the AFT formulation is the general manner in which it handles nonlinearities. In the problems considered by Urabe, as well as Ling and Wu, the nonlinearities are expressed analytically. In many applications, such as the elastic-plastic damper problem, the nonlinearities are hysteretic and analytical expressions may be difficult or impossible to obtain. In the AFT approach the user is provided with the full state history in order to compute the nonlinear forces in a user supplied subprogram, in a manner similar to that used in standard time integration codes.

Perhaps the most significant computational advantage of the AFT method is its technique for reducing the number of equations. As discussed previously, increasing the number of frequencies included in the solution can significantly increase ~would also be possible to solve this problem using Runge-Kutta time integration in conjunction with a two point boundary value formulation. In this approach the displacement and velocity of the mass at the beginning of the period are unknowns whose values must be determined so that the solution satisfies the periodicity requirements. Such an approach requires an iterative method for solving the resulting nonlinear constraint equations (for a discussion of this approach see, for example, Press et al., 1986, Chap. 16). This procedure is less standard than time integration and, consequently, was not investigated here. However, it could be significantly more efficient than time integration and more competitive with the computational efficiency of the AFT method. Again, the outcome of such a comparison would depend on the specific physical system considered, as well as on the complexity of the frequency content of the excitation.

the computational effort, but in many problems it may be necessary to consider both low and high frequencies. The AFT implementation permits the .user to include a relatively large number of frequencies and to set a "cutoff" fraction. While using direct iteration (which is computationally efficient), the algorithm scans all the frequency components to estimate which components lie above the cutoff fraction relative to the largest component. These components are considered ~'signifi cant.'' If it becomes necessary to employ an intensive equation solver, like Newton-Raphson, only the equations corresponding to the significant components are included. This approach also makes it easy for the user to examine tradeoffs between higher accuracy (low cutoff) and lower computation time (high cutoff).

Concluding Remarks

In the absence of aliasing, leakage, and roundoff errors the AFT method converges uniformly to the exact answer as Noo. Even when it cannot be eliminated, means are available for limiting the error so that good approximate solutions may be obtained. For example, in many problems the power in a signal decays asymptotically to zero as the frequency increases, so aliasing cannot be eliminated (without analog filtering). However, it is still possible to obtain exact solutions to within the precision of the computer by selecting il T so that power at frequencies exceeding (1r/ ilT) cannot numerically affect the solution. In the sample problem the error in the AFT method dropped off much more quickly than the error in numerical integration.

Solutions generated by the AFT method are conveniently presented in either the frequency or time domain, or both. The AFT method also lends itself to parametric and sensitivity studies because the solution for one set of parameters supplies a good estimate of the solution for a "neighboring" set of parameters, which aids in convergence of the iterations.

Relative to other frequency domain methods, the AFT method is amenable to more general, complicated nonlinearities and excitations. By switching to the time domain to evaluate the nonlinear forces, the AFT method can use the same description of the nonlinear forces used for numerical integration. The only constraint on the excitation is that its Fourier transform can be approximated using the FFT approach.

FFFTFig. 1

 1 Fig. 1 Direct iteration implementation of the AFT method

  Fig. 2 Elastic/perfectly-plastic friction damper; (a) graphic representa• lion and (b) sample output for harmonic input f(t) i

  Fig. 4 Amplitude and phase comparison between AFT method and fourth-order Runge•Kutta solutions lor three frequency excitation 103

Fig. 6

 6 Fig.6Error convergence comparison between AFT and fourth-order R-K methods

The Fourier Transform definition used here ~ X(w) = j: 00 x(t)fiwt dt and the inverse transform is x(t) = 1/ 2,-j_ 00 X(w)e-Jwtdw . Analogous definitions are employed in the FFT.

 4 The "w" notation in (4) is for conceptual purposes. For digital implementation (4) is cast in stricter vector notation-A(k)•X(k) =F(k) +11(k)-where the kth component of each vector represents a sampled version of the corresponding continuous function at the frequency "'k (ignoring scale factors and errors).

Acknowledgment

This work was supported by the National Science Foundation, grant number MSM 8703887.