Paulo Flores
email: pflores@dem.uminho.pt

A parametric study on the dynamic response of planar multibody systems with multiple clearance joints

Keywords: Multibody dynamics, Clearance joints, Contact forces, Friction effect

A general methodology for dynamic modeling and analysis of multibody systems with multiple clearance joints is presented and discussed in this paper. The joint components that constitute a real joint are modeled as colliding bodies, being their behavior influenced by geometric and physical properties of the contacting surfaces. A continuous contact force model, based on the elastic Hertz theory together with a dissipative term, is used to evaluate the intra-joint contact forces. Furthermore, the incorporation of the friction phenomenon, based on the classical Coulomb's friction law, is also discussed. The suitable contact-impact force models are embedded into the dynamics of multibody systems methodologies. An elementary mechanical system is used to demonstrate the accuracy and efficiency of the presented approach, and to discuss the main assumptions and procedures adopted. Different test scenarios are considered with the purpose of performing a parametric study for quantifying the influence of the clearance size, input crank speed and number of clearance joints on the dynamic response of multibody systems with multiple clearance joints. Additionally, the total computation time consumed in each simulation is evaluated in order to test the computational accuracy and efficiency of the presented approach. From the main results obtained in this study, it can be drawn that clearance size and the operating conditions play a crucial role in predicting accurately the dynamic responses of multibody systems.

Introduction

Over the last few decades, the Computer-Aided of Mechanical Systems emerged as an important scientific part of the Applied Mechanics field, with important applications in several branches of engineering [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF][START_REF] Haug | Computer-Aided Kinematics and Dynamics of Mechanical Systems -Volume I: Basic Methods[END_REF][START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Nikravesh | Planar Multibody Dynamics: Formulation, Programming, and Applications[END_REF][START_REF] Mccarthy | Geometric Design of Linkages: Systems and Controls[END_REF]. This has been made possible thanks to the impressive improvements of the computers at the level of both hardware and software, due to the creation of robust and accurate computational tools, and due to the demand for analysis of complex multibody mechanical systems. Decades ago, the design of machines and components was based on trial and error and knowledgeable craftsmanship. Later, algebraic methods for analysis eliminated part of the limitations of the trial and error and lead to documented methods used in the design of mechanical components [START_REF] Mccarthy | Geometric Design of Linkages: Systems and Controls[END_REF][START_REF] Flores | Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies[END_REF]. In today's industry there is little room for error and a great need for optimized and cost effective production of components and machines with high reliability and durability. Still, all areas of research are to some point based on the construction of models and, therefore, on assumptions and approximations. The development of faster data processing capabilities allowed more accurate and detailed methodologies and techniques inspired to improve the already existing formulations. The analysis of complex multibody systems is an area where in the past few decades faster data processing has lead to an increased research effort. This includes kinematic and dynamic analysis, synthesis and optimization of the motion forces in multibody systems [START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Nikravesh | Planar Multibody Dynamics: Formulation, Programming, and Applications[END_REF][START_REF] Mccarthy | Geometric Design of Linkages: Systems and Controls[END_REF][START_REF] Flores | Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies[END_REF].

All of the mechanical systems do not have perfect kinematic joints due to functional tolerances between the adjacent segments [START_REF] Flores | Dynamics of Multibody Systems with Spherical Clearance Joints[END_REF]. Either due to the loads carried by the mechanical devices or misalignments that are required for their operation, real joints must be lubricated [START_REF] Flores | Lubricated revolute joints in rigid multibody systems[END_REF] or include bushings [START_REF] Ambrósio | Sensitivity of a vehicle ride to the suspension bushing characteristics[END_REF], generally made with metals and polymers. By using rubber bushings a conventional mechanical joint is transformed into a joint with clearance allowing for the mobility of the over constrained system in which it is used [START_REF] Park | Effect of Steering-Housing Rubber Bushings on The Handling of a Vehicle[END_REF]. Therefore, appropriate tribological models must be devised in the framework of their application in general multibody systems [START_REF] Flores | A study on dynamics of mechanical systems including joints with clearance and lubrication[END_REF].

The characterization of the normal contact forces in the non lubricated joints is realized by using the continuous contact force model [START_REF] Lankarani | A Contact Force Model With Hysteresis Damping for Impact Analysis of Multibody Systems[END_REF] while their tangential forces are obtained by using appropriate friction force models [START_REF] Flores | Influence of the contact-impact force model on the dynamic response of multibody systems[END_REF]. Joints with rubber bushings are usually described by the methodology suggested by Park and Nikravesh [START_REF] Park | Effect of Steering-Housing Rubber Bushings on The Handling of a Vehicle[END_REF]. However, this is a linear model that does not include coupling between radial and axial or bending loading. A more advanced model for rubber bushings is realized by obtaining the static nonlinear response of the rubber bushings using finite elements and implementing the force-displacement curves as carpet plots in the computational code that supports the analysis [START_REF] Ambrósio | Improved bushing models for general multibody systems and vehicle dynamics[END_REF].

The serious consequences of the clearance joints on the behavior of the multibody systems have motivated various theoretical and experimental investigations over the last three decades [START_REF] Flores | Spatial rigid-multi-body systems with lubricated spherical clearance joints: modeling and simulation[END_REF][START_REF] Dubowsky | Dynamic Analysis of Mechanical Systems With Clearances, Part 1: Formulation of Dynamic Model[END_REF][START_REF] Dubowsky | Dynamic Analysis of Mechanical Systems With Clearances, Part 2: Dynamic Response[END_REF][START_REF] Dubowsky | Design and Analysis of Multilink Flexible Mechanism With Multiple Clearance Connections[END_REF][START_REF] Kakizaki | Modeling the Spatial Dynamics of Robotic Manipulators with Flexible Links and Joint Clearances[END_REF][START_REF] Bengisu | A Theoretical and Experimental Investigation of Contact Loss in the Clearances of a Four-Bar Mechanism[END_REF][START_REF] Bauchau | Modelling of Joints with Clearance in Flexible Multibody Systems[END_REF][START_REF] Innocenti | Kinematic Clearance Sensitivity Analysis of Spatial Structures With Revolute Joints[END_REF][START_REF] Wu | Interval Approach for the Modeling of Tolerances and Clearances in Mechanism Analysis[END_REF]. Furthermore, mechanical systems with clearance joints can have a predictable nonlinear dynamic response [START_REF] Farahanchi | Chaotic and Periodic Dynamics of a Slider-Crank Mechanism with Slider Clearance[END_REF][START_REF] Rhee | Dynamic Response of a Revolute Joint with Clearance[END_REF][START_REF] Flores | Dynamic behaviour of planar rigid multibody systems including revolute joints with clearance[END_REF]. This is an important feature for the design and control of these systems. Consequently, proper modeling of clearance joints in mechanical systems is required to achieve better understanding of the dynamic performance of the machines and mechanisms. This aspect gains paramount importance due to the demand for the proper design of the real joints in many industrial applications. The general purpose computational tools used for the design and analysis of multibody systems have a wide number of mechanical systems modeling features that require the description of rigid or flexible bodies for which geometry, mass, center of mass, moment of inertia, and other relevant properties are defined [START_REF] Ryan | ADAMS-Multibody System Analysis Software, Multibody Systems Handbook[END_REF][START_REF] Smith | DADS-Dynamic Analysis and Design System, Multibody Systems Handbook[END_REF]. The computational codes also provide a large library of kinematic joints that constrain relative degrees of freedom between connected bodies. The kinematic joints available in the commercial programs are represented as ideal joints, that is, there are no clearances or deformations in them.

Thus, modeling the dynamics of multibody systems with clearances and imperfections is a challenging issue in mechanical design and much work still remains to be done to achieve satisfactory modeling tools.

The primary objective of this work is to explore the mathematical modeling of real joints in the multibody systems. Furthermore, a parametric study of a planar mechanical system is performed in order to quantify the influence of the main parameters on the dynamic response of mechanical systems with multiple clearances joints. The selected parameters are the clearance size, the input crank speed, the number of joints modeled as clearance joints. During the numerical simulations the time of computation is evaluated in order to analyze the computational efficiency of the presented methodology. In general, the present research work is expected to be of great use in the analysis of multibody systems with clearance joints, namely in regard to the design and control tasks of such systems.

In a broad sense, there are two main approaches to model multibody systems with contact-impact events, namely the regularized models [START_REF] Gilardi | Literature survey of contact dynamics modeling[END_REF] and non-smooth dynamics formulation [START_REF] Flores | Modeling and analysis of rigid multibody systems with translational clearance joints based on the nonsmooth dynamics approach[END_REF]. The regularized approach is in fact a penalty method in which the contact forces and deformations are modeled by a set of spring-damper elements that represent the surface compliance of the contacting bodies [START_REF] Lankarani | A Contact Force Model With Hysteresis Damping for Impact Analysis of Multibody Systems[END_REF]. This approach is quite simple, very straightforward to implement and presents also a good computational efficiency. In addition, in the regularized approach, there are no impulses at the instant of contact, therefore, there is no need for impulsive dynamics calculations and the contact loss can be easily determined from position and velocity data [START_REF] Flores | Dynamic Analysis for Planar Multibody Mechanical Systems with Lubricated Joints[END_REF]. One of the disadvantages associated with this formulation is that it can introduce highly-frequency dynamics into the system, due to the presence of stiff springs in compliant surfaces. If these dynamics require the integration algorithm to take smaller steps, then the speed of simulation will be strongly penalized [START_REF] Flores | On the contact detection for contact analysis in multibody systems[END_REF]. Another drawback of the regularized approach is that the difficulty to chosen the contact parameters such as the equivalent stiffness or the degree of nonlinearity of the deformation, especially for complex contact scenarios [START_REF] Flores | Contact-impact analysis in multibody systems based on the nonsmooth dynamics approach[END_REF]. An alternative way to treat the contact-impact problems in multibody systems is to use the non-smooth dynamics approach, namely the Linear Complementary Problem (LCP) [START_REF] Pfeifer | Multibody dynamics with unilateral constraints[END_REF][START_REF] Anitescu | Time-stepping for three-dimensional rigid body dynamics[END_REF] and Differential Variational Inequality (DVI) [START_REF] Pang | Differential variational inequalities[END_REF][START_REF] Tasora | Large-scale parallel multi-body dynamics with frictional contact on the graphical processing unit[END_REF]. The complementarity formulations associated with the Moreau's time-stepping algorithm for contact modeling in multibody systems have used by many researchers [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementary systems in dynamics[END_REF][START_REF] Förg | Simulation of Unilateral Constrained Systems with Many Bodies[END_REF][START_REF] Pfeiffer | Numerical aspects of non-smooth multibody dynamics[END_REF]. Assuming that the contacting bodies are truly rigid, as opposed to locally deformable or penetrable bodies as in the penalty approaches, the complementarity formulations resolve the contact dynamics problem by using the unilateral constraints to compute contact impulses or forces to prevent penetration from occurring. Thus, at the core of the complementarity approach is an explicit formulation of the unilateral constraints between the contacting rigid bodies [START_REF] Brogliato | Numerical simulations of finite dimensional multibody nonsmooth mechanical systems[END_REF]. The basic idea of complementarity in unilateral multibody systems can be stated as for a unilateral contact either relative kinematics is zero and the corresponding constraint forces are zero, or vice versa. The product of these two groups of quantities is always zero. This leads to a complementarity problem and constitutes a rule which allows the treatment of multibody systems with unilateral constraints [START_REF] Glocker | Complementarity Problems in Multibody Systems with Planar Friction[END_REF][START_REF] Pang | Complementarity formulations and existence of solutions of dynamic multirigid-body contact problems with Coulomb friction[END_REF][START_REF] Trinkle | Dynamic multi-rigid-body systems with concurrent distributed contacts[END_REF][START_REF] Pfeiffer | The idea of complementarity in multibody dynamics[END_REF]. One of the first published works on the complementarity problems is due to Signorini [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF], who introduced an impenetrability condition in the form of a Linear Complementary Problem. Later, Moreau [START_REF] Moreau | Application of convex analysis to some problems of dry friction[END_REF] and Panagiotopoulos [START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF] also applied the concept of complementarity to study nonsmooth dynamic systems. Pfeiffer and Glocker [START_REF] Pfeifer | Multibody dynamics with unilateral constraints[END_REF] extended the developments of Moreau and Panagiotopoulos to multibody dynamics with unilateral contacts, being the complementarity considered of paramount importance. Indeed, complementarity problems proved to be a very useful way to formulating problems involving discontinuities [START_REF] Kwak | Complementarity Problem Formulation of Three-Dimensional Frictional Contact[END_REF][START_REF] Anitescu | Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementarity Problems[END_REF][START_REF] Brogliato | Some perspectives on the analysis and control of complementarity systems[END_REF][START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF]. In turn, the DVI has been recognized to be a powerful tool to deal with multiple contact problems in multibody dynamics. This approach has the advantage that it does not need the use of small time steps as in the case of penalty approaches, which means that simpler integrator schemes can be used such as the Euler method [START_REF] Pang | Differential variational inequalities[END_REF]. However, the algorithmic procedures that results from DVI approach is of great complexity. This formalism has been used with success by Tasora et al. [START_REF] Tasora | Large-scale parallel multi-body dynamics with frictional contact on the graphical processing unit[END_REF] to model and analyze multibody systems involving hundreds of thousand contacts. The DVI approaches are also interesting in the measure that they can easily deal with friction problems without need to modify Coulomb's friction law. In short, the different methods to deal with contact-impact events in multibody systems have inherently advantages and disadvantages for each particular application. To the belief of author, none of the formulations briefly described above can a-priori be said to be superior compared to other for all applications. It is a fact that a specific multibody problem might be easier to describe by one formulation, but this does not yield a general predominance of this formulation in all situations.

Modeling revolute clearance joints

In standard multibody models it is assumed that the connecting points of two bodies, linked by an ideal revolute joint, are coincident. The introduction of the clearance separates these two points. Figure 1 shows a typical connection with clearance joints found in planar multibody systems, where the clearance size is strongly exaggerated in order to illustrate the phenomena associated with the revolute joints with clearance, namely, the bouncing effect. In a revolute clearance joint, when contact exists between the journal and bearing, a contact-impact force is applied perpendicular to the plane of collision. The force is typically applied as a spring damper element. If this element is linear, the approach is known as the Kelvin Voigt model [START_REF] Timoshenko | Theory of Elasticity[END_REF]. If the relation is nonlinear the model is generally based on the Hertz contact law [START_REF] Hertz | On the contact of solids -On the contact of rigid elastic solids and on hardness[END_REF]. In the contact or following mode, the journal and the bearing are in permanent contact and a sliding motion related to each other is assumed to exist. In this mode the relative penetration depth varies along the circumference of the bearing. In practice, this mode is ended when the journal and bearing separate from each other and the journal enters in free flight mode.

In the free flight mode, the journal can move freely inside the bearing boundaries, that is, the journal and the bearing are not in contact, hence, and no reaction force develop at the joint. In the impact mode, which occurs at the termination of the free flight mode, impact forces are applied and removed in the system. This mode is characterized by a discontinuity in the kinematic and dynamic characteristics, and a significant exchange of momentum occurs between the two impacting bodies. At the termination of the impact mode, the journal can enter either a free flight or a following mode. Figure 3 shows a generic configuration of a revolute joint with clearance in a multibody system. In the dynamic simulation, the behavior of the revolute clearance joint is treated as an oblique eccentric impact between the journal and the bearing. The mechanics of this type of impact involves both the relative normal velocity and the relative tangential velocity [START_REF] Zukas | Impact Dynamics[END_REF]. When the impact occurs, an appropriate contact law must be applied and being the resulting forces introduced as generalized forces in the system's equations of motion [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF].

Regarding to Fig. 3, the relative penetration vector between journal and bearing walls can be defined as

2 B J d d - ⎛ ⎞ = -⎜ ⎟ ⎝ ⎠ δ e n ( 1
)
where e is the eccentricity vector, that is, the vector that defines the distance between journal and bearing centers, d B and d J are the bearing and journal diameters and n is the vector that defines the normal direction of the plane of collision.

The magnitude of the penetration depth can be evaluated as

= e -c δ (2)
in which c represents the radial clearance defined as

2 B J d -d c = (3)
Since the most suitable and sophisticated contact-impact force models are dependent on the contact velocities, it is important to evaluate these velocities in order to account for the dissipative effects during the contact-impact process [START_REF] Flores | A study on dynamics of mechanical systems including joints with clearance and lubrication[END_REF]. In particular, in the continuous force contact model it is necessary to calculate the relative velocity between impacting surfaces. The relative velocity vector between bearing and journal can be expressed by

( ) ( ) ( ) C C B J C C B J B J d dt ω ω - = + - - r r v n r r t (4)
where C B r and C J r are the position vectors of the contact points in the bearing and journal, n and t are the normal and tangential vectors, and ω B and ω J are the bearing and journal angular velocities. It should be highlighted that vector t is obtained by rotating the vector n in the counter clockwise direction by 90º.

The relative velocity between the contact points is projected onto the plane of collision and onto the normal plane of collision, yielding a relative tangential velocity, v T , and a relative normal velocity, v N . The normal relative velocity determines whether the contact bodies are approaching or separating. Similarly, the tangential relative velocity determines whether the contact bodies are sliding or sticking. The relative scalar velocities, normal and tangential to the plane of collision are found by projecting the relative impact velocity onto each one of these directions

T N v = v n (5) T T v = v t (6)
Although a revolute joint with clearance does not constrain any degree of freedom from the mechanical system, as the ideal joint does, it imposes some restrictions, limiting the journal to move within the bearing. Thus, when the clearance is present in a revolute joint, two kinematic constraints associated with an ideal joint are removed and two degrees of freedom are introduced instead. The dynamics of the joint is then controlled by contact-impact forces between the journal and bearing. Thus, whilst a perfect revolute joint in a mechanical system imposes kinematic constraints, a revolute clearance joint leads to force constraints. Therefore, joints with clearance can be defined as force-joints instead of kinematic joints.

When contact between the journal and bearing takes place, contact-impact forces act at the contact points. The contributions of these contact-impact forces to the external forces applied to the mechanical system are found by projecting them onto the X and Y directions. Since these forces do not act through the centre of mass of the bearing and journal bodies the moment components for each body need to be evaluated. The suitable numerical models for normal and tangential contact-impact forces are presented and discussed in the following section.

In short, in dynamic analysis, the deformation is known at every time step from the configuration of the system and the forces are evaluated based on the state variables. With the variation of the contact force during the contact period, the dynamic system's response is obtained by simply including updated forces into the equations of motion. Since the equations of motion are integrated over the period of contact, this approach results in a rather accurate response. Furthermore, this methodology accounts for the changes in the system's configuration during the contact periods. The methodology presented in the previous sections is quite simple and generic to apply in the dynamic modeling and analysis of any mechanical system. For details on this issue, the interested reader is referred to the work developed by Flores [START_REF] Flores | Dynamic Analysis of Mechanical Systems with Imperfect Kinematic Joints[END_REF].

Numerical models for contact forces

Impact, such as it happens in a revolute clearance joint, is one of the most common types of dynamic loading conditions that give rise to impulsive forces, which in turn excite higher vibration modes and affect the dynamic characteristics of the mechanical systems. Thus, for a revolute joint with clearance, the contact between the journal and bearing can be modeled by the well known Hertz contact law [START_REF] Hertz | On the contact of solids -On the contact of rigid elastic solids and on hardness[END_REF] n N F = Kδ [START_REF] Flores | Dynamics of Multibody Systems with Spherical Clearance Joints[END_REF] where K is the stiffness coefficient and δ is the relative penetration. The exponent n is set to 1.5.

The parameter K depends on the material and geometric properties of the contacting surfaces.

For instance, for two spheres in contact the generalized stiffness coefficient is a function of the radii of the spheres i and j and the material properties as [START_REF] Goldsmith | Impact -The Theory and Physical Behaviour of Colliding Solids[END_REF] 1 2 4 3( )

i j i j i j R R K R + R σ σ ⎡ ⎤ = ⎢ ⎥ + ⎢ ⎥ ⎣ ⎦ (8)
where the material parameters σ i and σ j are given by,

2 1 k k k E ν σ - = , (k=i,j) (9)
the quantities ν k and E k are the Poisson's ratio and the Young's modulus associated with each sphere, respectively.

Hertz contact law given by Eq. ( 7) is a pure elastic model, that is, it does not include any energy dissipation. Lankarani and Nikravesh [START_REF] Lankarani | A Contact Force Model With Hysteresis Damping for Impact Analysis of Multibody Systems[END_REF] extended the Hertz contact law to include energy loss due to internal damping as follows

4

n e N c F K δ δ δ - ⎡ ⎤ - = + ⎢ ⎥ ⎣ ⎦ (10)
where the stiffness coefficient K can be evaluated by Eqs. ( 8) and ( 9), c e is the restitution coefficient, δ is the relative penetration velocity and ( ) δ -is the initial impact velocity.

The Coulomb's friction law of sliding friction can represent the most fundamental and simplest model of friction between dry contacting surfaces. When sliding takes place, the Coulomb law states that the tangential friction force F T is proportional to the magnitude of the normal contact force, F N , at the contact point by introducing a coefficient of friction c f [START_REF] Greenwood | Principles of Dynamics[END_REF]. The Coulomb's friction law is independent of relative tangential velocity. In practice this is not true, because friction forces can depend on many parameters such as material properties, temperature, surfaces cleanliness, and velocity of sliding. Therefore, a continuous friction force-velocity relationship is desirable. Furthermore, the application of the original Coulomb's friction law in a general-purpose computational program may lead numerical difficulties because it is a highly nonlinear phenomenon that may involve switching between sliding and stiction conditions. In order to avoid such difficulties, a modified Coulomb law is used [58]

T T f d N T F c c F = - v v ( 11
)
where c f is the friction coefficient, F N is the normal force, v T is the relative tangential velocity and c d is a dynamic correction coefficient, which is expressed as

0 0 0 1 1 0 1 0 if if 1 if T T d T T v v v v c v v v v v v v ⎧ ≤ ⎪ - ⎪ = ≤ ≤ ⎨ - ⎪ ≥ ⎪ ⎩ (12)
where v 0 and v 1 are given tolerances for the tangential velocity [START_REF] Flores | Dynamic Analysis of Mechanical Systems with Imperfect Kinematic Joints[END_REF]. The dynamic correction factor c d prevents that the friction force changes direction in the presence of almost null values of the tangential velocity, which would be perceived by the integration algorithm as a dynamic response with high frequency contents, forcing it to reduce the time step size.

Multibody systems formulation

The equations of motion for a dynamic multibody system subjected to holonomic constraints can be stated in the form [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF]

T ⎡ ⎤ ⎧ ⎫ ⎧ ⎫ = ⎢ ⎥⎨ ⎬ ⎨ ⎬ ⎢ ⎥⎩ ⎭ ⎩ ⎭ ⎣ ⎦ q q M Φ q g Φ 0 λ γ (13)
with the reference frame placed at the center of mass for each body, M is the system mass matrix, q Φ is the Jacobian matrix of constraint equations, the vector q contains the generalized state accelerations, λ is the vector that contains the Lagrange multipliers, g is the vector of generalized forces and γ is the vector of quadratic velocity terms that is used to describe Coriolis and centrifugal terms in the equations of motion [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF]. Equation ( 13) is formed as a combination of the equations of motion and kinemtaic constraint equations, often referred to as a mixed set of differential and algebraic equations.

A set of initial conditions (positions and velocities) is required to start the dynamic simulation. The selection of the appropriate initial conditions plays a key role in the prediction of the dynamic performance of mechanical system [START_REF] Nikravesh | Initial condition correction in multibody dynamics[END_REF]. In the present work, the initial conditions are based on the results of kinematic simulation of mechanical system in which all the joints are assumed to be ideal, that is, without clearance. The subsequent initial conditions for each time step in the simulation are obtained in the usual manner from the final conditions of the previous time step [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF]. In order to stabilize or keep under control the constraints violation, Eq. ( 13) is solved using the Baumgarte stabilization technique [START_REF] Baumgarte | Stabilization of Constraints and Integrals of Motion in Dynamical Systems[END_REF]. In turn, the integration process is performed using a predictor-corrector algorithm with both variable step size and order [START_REF] Gear | Numerical Initial Value Problems in Ordinary Differential Equations[END_REF].

The use of numerical algorithms with automated adjust step size is particularly important in contact problems whose dynamic response is quite complex due to the suddenly change in kinematic configuration. In such events, the use of a constant time step is computationally inefficient and the system could be overlooked due to insufficient time resolution. Thus, automated time step size adaptability is therefore a crucial part of the dynamic solution procedure. Moreover, the abrupt configuration changes caused by rapidly variation of contact forces are source of stiffness, since the natural frequency of the system is widely spread. Thus, the time step size must be adjusted in order to capture the fast and low components of the system response [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF].

Results and discussion

This section contains extensive results obtained from computational simulations of a planar slider-crank mechanism with multiple clearance joints when subjected to different test scenarios in order to carry out a parametric study. This study takes into account the main functional parameters of the slider-crank mechanism, chiefly, clearance size, input crank speed and number of joints considered as clearance joints. The computational time consumed during numerical simulations is used as a measured of the computational efficiency of the proposed methodologies.

Description of the slider-crank mechanism

Figure 4 shows the configuration of the slider-crank mechanism selected, which consists of four rigid bodies that represent the crank, connecting rod, slider and ground, one ideal revolute joint and one ideal translational joint. Two revolute joints with clearance exist between the crank and connecting rod, and between the connecting rod and slider. The length and inertia characteristics of this mechanism are listed in Table 1. Due to the presence of the two non-ideal joints this system has five degrees of freedom. The acceleration due to gravity is taken as acting in the negative Y direction and the mechanism is defined as moving in a vertical plane. Tab. 1: Geometric and inertia properties of the slider-crank mechanism.

The crank, which is the driving link, rotates with a constant angular velocity of 5000 rpm.

The initial configuration of the mechanism is defined with the crank and the connecting rod collinear, and the journal and bearing centers coincident. Furthermore, the initial positions and velocities necessary to start the dynamic analysis are obtained from kinematic simulation of the slider-crank mechanism in which all the joints are considered to be ideal. Table 2 gives the main parameters used for the different models, required to characterize the problem and for the numerical methods, required to solve the system dynamics. The diametric clearance size of the non-ideal revolute joints is taken to be equal to 0.05 mm, which corresponds to the actual clearance size in a typical journal-bearing with nominal dimensions used in the present work. For example, for a journal-bearing in which the nominal diameter falls into the interval between 18 and 30 mm, the minimum and maximum recommended diametric clearance sizes are 0.02 and 0.06, respectively [START_REF] Ansy | Dimensional and Tolerancing[END_REF].

In what follows, several numerical results are presented and analyzed to demonstrate the computational implementation and efficiency of the proposed methodologies. The dynamic response of the slider-crank mechanism is quantified by plotting the slider velocity and acceleration, joint reaction force and the reaction moment that acts on the crank. In addition, the journal center orbit inside the bearing boundaries and the maps that relates the eccentricity and eccentricity velocity are plotted too. The global results are relative to two complete crank rotations after steady state has been reached. Furthermore, the results for this system are compared with the dynamic behavior of the slider mechanism with ideal joints only.

Influence of the clearance size: one clearance joint

In this section, the influence of the clearance size on the dynamic behavior of the slidercrank mechanism is investigated. In the simulation of the mechanism, only one joint is modeled as clearance joint, namely, the joint that connects the connecting rod and slider.

Figure 5 shows the results for the case in which the diametric clearance size is equal to 0.05 mm and crank speed of 5000 rpm, as it was presented in the previous section. Note that the results are plotted against those obtained for ideal joint, being reported for the two full crank rotations after steady-state has been reached. The results plotted are the slider velocity and acceleration, the reaction force that is produced in the clearance joint, the crank reaction moment that acts on the crank, the trajectories of the journal center relative to the bearing center and, finally, the Poincaré map. The eccentricity and eccentricity velocity in the clearance joint are the parameters selected to plot this map.

In Fig. 5(a) it is observed that the existence of a clearance joint does not influence the slider velocity in a significant manner. In sharp contrast, the slider acceleration is strongly affected by the contact-impact forces that takes place between journal and bearing contacting surfaces in the clearance joint. The slider acceleration is subjected to peaks caused by contact forces that are propagated through the rigid bodies of the mechanism, as perceived in Fig. 5(b) where the acceleration of the slider is displayed. The same phenomena can be observed in the curve of joint reaction force and crank moment represented by Figs. 5(c) and 5(d), respectively.

As far as the trajectory of the journal center relative to the bearing center is concerned, only one type of motion between the two bodies is observed, namely the permanent or continuous contact mode. The relative penetration depth between the journal and bearing is visible by the lines of the journal path that are plotted outside the clearance circle, as illustrated in Fig. 5(e). From the Poincaré map of Fig. 5(f) it can be concluded that the behavior of the slider-crank mechanism modeled with a clearance joint tends to be aperiodic, or even chaotic, due to the nonlinear nature of the contact-impact phenomena verified between the joint components. This conclusion can easily be explained based on the spread out nature of the Poincaré map. In fact, non periodic responses are extremely sensitive to initial conditions and are densely filled by orbits or points in the Poincaré map. A complicated looking phase in a Poincaré map is one indicator of chaotic motion. Quasi-periodic orbits fill up the Poincaré maps as the chaotic orbits, but they do so in a fully predictable manner since there is no sensitive dependence on the initial conditions [START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF]. are clearly exaggerated in order to study the influence of the variation of the clearance size on the system's dynamic response. Figure 6 shows that when the clearance size is increased the dynamic system behavior tends to be non periodic. This idea is supported by the higher peaks in the slider acceleration curves, which means that impacts followed by rebounds take place, instead of continuous or permanent contact between the journal and bearing walls. When the clearance is small, the system response tends to be closer to the ideal response meaning that the journal and bearing experiment a smaller number of impacts and the journal follows the bearing wall. Hence, the clearance joint behavior tends to be periodic instead of nonlinear or chaotic.

This observation can be drawn by examination of Fig. 7 relative to the joint reaction force. The level of peaks force increase with the increasing the clearance size, owing to the higher number of impacts that take place. This conclusion can be confirmed in the journal center trajectories plots and in the respective Poincaré maps displayed in Figs. 8 and9, respectively. Observing Fig. 8 it can be drawn that for larger clearances there are some periods in which journal is in freeflight motion inside the bearing boundaries. In fact, for clearance size equal 0.2 and 0.5 mm the three types of journal motion are observed, namely, the free-flight motion, the impact mode and the permanent or continuous contact mode. From the Poincaré maps analysis, the slider-crank mechanism behavior can easily be characterized, and it is possible to distinguish between periodic, quasi-periodic and chaotic responses, when the clearance size increases. This chaotic response suggests that impacts followed by some rebounds take place. Figure 9(a) shows a quasiperiodic motion, because the orbits fill up the Poincaré maps in a fully predictable manner, thus, there is no sensitive dependence on the initial conditions. It is clear that when the clearance is reduced the dynamic response tends to be periodic or regular, which indicates that the journal follows the bearing wall. In fact, in mechanical systems nonlinearities arise from intermittent motion, clearance joints, friction effect, and contact forces, among others. The relation between the clearance size and the type of motion observed is clearly identified from plots in Figs. 8 and9.

Influence of the input crank speed: one clearance joint

The influence of the input crank speed on the dynamic behavior of the slider-crank mechanism is investigated in this section. Again, in the simulation of the mechanism, only one joint is modeled as clearance joint, that is, the joint that connects the connecting rod and slider. Figure 10 illustrates the global results for the simulation in which the diametric clearance size is equal to 0.05 mm and the input crank speed is of 50 rpm. The effect of the clearance joint is evident on the system's response, namely, in high peaks on the slider acceleration, joint reaction force and crank moment. These peaks are associated with the contact-impact forces developed in clearance joint. Since in this simulation the input crank speed is relatively low (50 rpm) the journal is almost always in permanent contact with bearing wall, as illustrated in Fig.

10(e). Hence, the system's behavior tends to be quasi-periodic as it is visible in the respective Poincaré map of Fig. 10(f), where a strange attractor is evident and is also a strong indicator for this quasi-periodic response.

The influence of the input crank speed on the dynamic response of the slider-crank mechanism is illustrated in Figs. 11 to 14. Four different values for input crank speed were selected, chiefly, 50, 500, 1000 and 5000 rpm. These different input crank speeds imply the use of different scales for results plotted. By and large, the influence of the clearance joint on the output results is similar for all crank speed simulations, that is, the slider acceleration and joint reaction force curves present peaks due to collisions between the journal and bearing. From Fig. 13 it can be observed that the relative deformation between the journal and bearing increases with the crank speed, because the impacts are severer. In addition, for higher values for crank speed, the system's behavior tends to be nonlinear, as it is quite visible in Poincaré maps of Fig. Figure 15 shows the global results for the case in which the joint that connects the crank and connecting rod is modeled as clearance joint with 0.05 mm for diametric clearance size, being the system operating at 5000 rpm. In general, the results, as expected, are dependent on the contact-impact phenomena that take place between the journal and bearing. From Fig. 15(e) it can be observed that the journal and bearing are in permanent or continuous contact mode. The system's response is also chaotic, as illustrated in Fig. 15(f).

Concluding Remarks

A general methodology for modeling and simulating multiple clearance joints in rigid multibody systems was presented and discussed throughout this work. The equations of motion that govern the dynamic response of the general multibody systems incorporate the contact forces due to the collisions of the bodies that constitute the clearance joints. A continuous contact force model provides the intra-joint contact-impact forces that develop during the normal operations of the mechanisms. A suitable model for revolute clearance joints is embedded into the multibody systems methodology. An elementary slider-crank mechanism was used to demonstrate the methodology proposed and to perform a parametric study. The main functional parameters were used for this study, namely, the diametric clearance size, the input crank speed and the number of joints modeled as clearance joints. In addition, the computation time was evaluated in each simulation in order to test the numerical efficiency. In general, the presented methodology is simple, general and cheap from computational view point. This particular issue has not been addressed before in the literature.

From the numerical simulations performed, it can be concluded that there are considerable changes in joint reaction forces and on the crank moments due to the effect of clearances. These changes are mainly due to the peaks observed in the slider acceleration curves. It is also important to note that the magnifications of forces and moments in the slider-crank mechanism as affected by revolute clearance joints essentially depend on the clearance size and input crank speed. From the main results obtained in this work, it can be also observed that multibody systems with clearance joints are well known as nonlinear dynamic systems that, under certain conditions, exhibit a chaotic response. From the results presented above, it is found that the dynamics of the revolute clearance joint in mechanical systems is sensitive to the clearance size value operating conditions. With a small change in one of these parameters the response of the system can shift from chaotic to periodic and vice-versa.

Finally, it should be highlighted that the results presented in this paper represent an upper bound of the joint reaction forces and crank moments due to the existence of clearance joints, since the elasticity of bodies and lubrication action in the joints were not included in the analysis.

These effects tend to reduce the values of these joint reaction forces and crank moments.

Moreover, lubrication could also change the entire behavior of the mechanism.

Figure 1 .

 1 Figure 1. Typical connection with clearance joints found in planar multibody systems.

Figure 2

 2 Figure 2 depicts a revolute joint with clearance, that is, the so-called journal bearing, where the difference in diameters between the bearing and the journal defines the diametric clearance. Several published research works focused on the different modes of motion of the journal inside the bearing boundaries. Most of these consider a three mode model for predicting the dynamical response of articulated systems with revolute clearance joints [55-61]. The three different modes of journal motion inside the bearing are the contact or following mode, the free flight mode, and the impact mode, which are illustrated in Fig. 2.

Figure 2 .

 2 Figure 2. Types of journal motion inside the bearing boundaries.

Figure 3 .

 3 Figure 3. Generic configuration of a revolute joint with clearance in a multibody system.

Figure 4 .

 4 Figure 4. Slider-crank mechanism with two revolute joints with clearance.

Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .

 56789 Figure 5. (a) Slider velocity; (b) Slider acceleration; (c) Joint reaction force at the clearance joint; (d) Crank moment required to maintain its angular velocity constant; (e) Journal center trajectory relative to the bearing center; (f) Poincaré map.

Figure 10 .

 10 Figure 10. (a) Slider velocity; (b) Slider acceleration; (c) Joint reaction force at the clearance joint; (d) Crank moment required to maintain its angular velocity constant; (e) Journal center trajectory relative to the bearing center; (f) Poincaré map.

Figure 11 .Figure 12 .Figure 13 .Figure 14 . 5 . 4 .

 1112131454 Figure 11. Slider acceleration for different crank speed values: (a) n=50 rpm; (b) n=500 rpm; (c) n=1000 rpm; (d) n=5000 rpm.

Figures 16 to 19

 19 Figures 16 to 19 illustrate the results obtained when the mechanism is modeled with one, two and three clearance joints, being the global system's response similar, namely in what concerns with peaks on the slider acceleration and joint reaction force. It have to be noted that the journal trajectories and Poincaré maps have similar look for all the cases simulated, as illustrated in Figs. 18 and 19.

Figure 15 .Figure 16 .

 1516 Figure 15. (a) Slider velocity; (b) Slider acceleration; (c) Joint reaction force at the clearance joint; (d) Crank moment required to maintain its angular velocity constant; (e) Journal center trajectory relative to the bearing center; (f) Poincaré map.

Figure 17 .Figure 18 .Figure 19 .Figure 20 .

 17181920 Figure 17. Joint reaction force for different number of clearance joints: (a) One clearance joint between connecting rod and slider; (b) One clearance joint between crank and connecting rod; (c) One clearance joint between crank and connecting rod and one clearance joint between connecting rod and slider; (d) All three revolute joints have clearance.

  Tab. 2: Parameters used in the dynamic simulations of the slider-crank mechanism with clearance joints.

	Nominal bearing diameter	20.0 mm	Poisson's ratio	0.3
	Standard diametric clearance	0.05 mm	Integrator algorithm	Gear
	Coefficient of restitution	0.9	Baumgarte -α, β	5
	Coefficient of friction	0.01	Integration step size	1×10 -6 s
	Young's modulus	207 GPa	Total time of simulation	10.0 s

Acknowledgements

This work is supported by the Portuguese Foundation for the Science and Technology under the research project BIOJOINTS (PTDC/EME-PME/099764/2008). The author would also like to acknowledge the anonymous reviewers for their insightful comments and suggestions on an earlier draft of this paper.


Fernando B Duarte
email: fduarte@mat.estv.ipv.pt

J Tenreiro

Machado • Gonçalo

Monteiro Duarte
email: monteiro.duarte@gmail.com

J A Tenreiro Machado

Monteiro Gonçalo

Duarte

Dynamics of the Dow Jones and the NASDAQ Stock Indexes

Keywords: Pseudo-Phase Plane, Fourier Transform, Power Law, Fractional Calculus

The goal of this study is the analysis of the dynamical properties of financial data series from worldwide stock market indices. We analyze the Dow Jones Industrial Average ( ∧ DJI) and the NASDAQ Composite ( ∧ IXIC) indexes at a daily time horizon. The methods and algorithms that have been explored for description of physical phenomena become an effective background, and even inspiration, for very productive methods used in the analysis of economical data. We start by applying the classical concepts of signal analysis, Fourier transform, and methods of fractional calculus. In a second phase we adopt a pseudo phase plane approach.

Introduction

The study of fractional order systems has received considerable attention, due to the fact that many physical systems are well characterized by fractional models. The importance of fractional order mathematical models is that it can be used to make a more accurate description and it can even give a deeper insight into the processes underlying long range memory behavior [START_REF] Tenreiro | A probabilistic interpretation of the fractiional-order differentiation[END_REF][START_REF] Lima | Fractional dynamics in mechanical manipulation[END_REF][START_REF] Jiang | Bubble diagnosis and prediction of the 2005-2007 and 2008-2009 chinese stock market bubbles[END_REF]. It seems that there are many distinct analogies between the dynamics of complex physical and economical or even social systems. The methods and algorithms that have been explored for description of physical phenomena become an effective background and inspiration for very productive methods used in the analysis of economical data [START_REF] Plerou | Econophysics: financial time series from a statistical physics point of view[END_REF][START_REF] Nigmatullin | Universal distribution function for the strongly-correlated fluctuations: General way for description of different random sequences[END_REF][START_REF] Rui | A data-reconstructed fractional volatility model[END_REF][START_REF] Mendes | A fractional calculus interpretation of the fractional volatility model[END_REF].

In this paper we study the Dow Jones and the NASDAQ indexes at a daily time horizon at the closing.

People on Wall Street found it difficult to analyze the daily jumble of up-a-quarter and down-an-eighth, or whether stocks generally were rising, falling or staying even. Charles Dow a journalist, neither financier nor broker, devised his stock average to make sense of this confusion. He began in 1884 with eleven stocks, most of them railroads. Railroads were among the biggest and sturdiest companies in America at that time, which is why they dominated Dow's first average. Few stocks of industrial companies were publicly traded, and those were considered highly speculative. On May 26, 1896, he introduced the industrial average. In October of that year, Dow's original average shed the last of its non-railroad stocks and became the twenty-stock railroad average. To complete this line of history, the utilities average came along in 1929 and the railroad average was renamed the transportation average in 1970. Nowadays, of course, there are plenty of indicators to tell investors what the stock market is doing. The Dow Jones Industrial Average is in sync with other major market barometers. That is true despite the difference in computation methods; the Dow is unweighted while almost all other indexes weight their stocks by market capitalization, which is price times shares outstanding. It is also true despite the fewer number of stocks in the Dow [START_REF] Dow | the History of the Dow Jones Averages[END_REF].

The NASDAQ (National Association of Securities Dealers Automated Quotation) Stock Market, founded in 1971, was the world's first electronic stock market. The purpose of its founding was to popularize the OTC (over-the-counter) securities market which, up to that point, had been relatively unknown and unused by many stock players. With its first day of trading on February 8, 1971, the NASDAQ system displayed quotes for over 2,500 over-the-counter stocks. The NASDAQ stock market is full of technology stocks of up-and-coming companies, some with expensive stock prices and some for just pennies. It continues to be America's most popular market in a day and age when technology still seems to be the wave of the future [2].

Bearing these ideas in mind this paper is organized as follows: sections 2 and 3 present some fundamental concepts, and the dynamical analysis, respectively. The existence of a power law relationship typical of systems with fractional calculus is shown. Finally, section 4 draws the main conclusions and addresses perspectives towards future developments.

Fundamental concepts

In this section we present a review of fundamental concepts involved in the experiments. The technique used to determine the fractional behavior of the Dow Jones [3] and NASDAQ [4] index signals is based on the slope of their trendlines in the frequency spectrum. Additionally, the pseudo phase space (PPS ) is obtained using the method of the time delays [START_REF] Abarbanel | Tsimring: The analysis of observed chaotic data in physical systems[END_REF][START_REF] Feeny | Fractional derivatives applied to phase-space reconstructions[END_REF].

Fourier transforms

In order to examine the behavior of the signal spectrum, a power law trendline is superimposed to the Fourier Transform (FT) signal. To study the signal spectrum, we approximate the modulus of the FT amplitude through power functions

|F{x(t)}| ≈ p ω q , {p, q} ∈ IR ( 1
)
where F is the Fourier operator, x(t) is the index and t is time, ω is the frequency, p a positive constant that depends on the signal amplitude and q is the trendline slope [START_REF] Gabaix | A unified econophysics explanation for the power-law exponents of stock market activity[END_REF]. According to the values of q the signals can exhibit an integer or fractional order behavior.

The standard Fourier transform describing the data in the 'Fourier domain' is precise in frequency, but not in time. Small changes in the signal at one location cause change in the Fourier domain globally. It is of interest to have transformed domains that are simultaneously precise, both in time and frequency domains.

The Windowed Fourier transforms (WFT) are important in providing simultaneous insight in time and frequency behavior of the signal. A window function is a function that is zero-valued outside of some chosen interval. When a signal is multiplied by a window function, the product is also zero-valued outside the interval: all that is left is the "view" through the window. In this paper, among the several window functions, we adopt the Gaussian window for a sliding-window Fourier transform. The coefficients of a Gaussian window are computed from the following equation:

W (t) = e -1 2 (α t T /2 ) 2 (2)
where T is the window length, -T 2 ≤ t ≤ T 2 and α is the reciprocal of the standard deviation. The width of the window is inversely related to the value of α: a larger value of α produces a narrower window.

Pseudo Phase Plane

An essential problem in nonlinear time series analysis is to determine whether or not a given time series is a deterministic signal from a low-dimensional dynamical system. If it is, then further questions of interest are: what is the dimension of the phase space supporting the data set? Is the data set chaotic? The key to answering these questions is embodied in the method of phase space reconstruction, which has been rigorously justified by the embedding theorems. Takens' embedding theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF][START_REF] Packard | Geometry from a time series[END_REF] asserts that if a time series is one component of an attractor that can be represented by a smooth d-dimensional manifold (with d an integer) then the topological properties of the signal are equivalent to the topological properties of the embedding formed by the m-dimensional phase space vectors

y(t) = [s(t), s(t + τ ), s(t + 2τ ), • • • , s(t + (m -1)τ )] (3)
whenever m > 2d + 1. In equation (3) τ is called the delay time and d is the embedding dimension. Different choices of t and m yield different reconstructed trajectories. The vector y(t) can be plotted in a d-dimensional space forming a curve in the PPS. If d = 2 we get a two-dimensional time delay space where the signal {s(t), s(t + τ )} is related to the model {s(t), ṡ(t)}. In this case the PPS it is called pseudo phase plane (PPP ). We expect, with the PPP of the signal, to draw conclusions about the system dynamics [START_REF] Feeny | Fractional derivatives applied to phase-space reconstructions[END_REF].

There have been various proposals for choosing an optimal delay time, τ , for topological properties based on the behavior of the autocorrelation function. These include the earliest time t at which the autocorrelation drops to a fraction of its initial value or has a point of inflection. These definitions seek to find times where linear correlations between different points in the time series are negligible, but they do not rule out the possibility of more general correlations. Some argue that a better value for τ is the value that corresponds to the first local minimum of the mutual information. The mutual information is a measure of how much information can be predicted about one time series point giving full information about the other. The values of τ at which the mutual information has a local minimum are equivalent to the values of τ at which the logarithm of the correlation sum has a local minimum. It is not clear which method, if any, is superior for all topological properties. However, optimal values based on the behavior of the autocorrelation function is the easiest to compute and is adopted in our experiments [START_REF] Francisco Martins De Lima | Análise dinâmica de vibrações em manipuladores robóticos[END_REF].

Dynamics of financial indexes

In this section we study numerically the signals corresponding to the Dow Jones and the NASDAQ indexes. For both signals we study the fractional behavior and the PPP reconstruction.

Figures 1 and2 depict the time evolution of the two indexes with the well-know noisy of chaotic-like characteristics.

Fourier Analysis

Figure 3 shows the amplitude of the FT for the Dow Jones (|F{x D (t)}|) signal index. A trendline is calculated and it is superimposed over the signal. For the Dow Jones index the slope yields q = -0.79, revealing a fractional order behavior. Figure 4 shows the amplitude of the FT for the NASDAQ |F{x N (t)}| signal index. A trendline is also calculated and superimposed over the signal. For the NASDAQ index the slope yields q = -0.83 and reveals, also, a fractional order behavior. In fact, in both cases we get a fractional noise in between white and pink noise, corresponding to a considerable volatility.

We now adopt the windowed Fourier transform and we consider α = 2.5, T = 366 days (one year) and that two consecutive windows are superimposed over a range of window length of β = 50%.

For the Dow Jones index we considerer a total of the 110 windows centered at t = 183λ D days for λ D = 1, 2, • • • , 110 and for the NASDAQ index we adopt a total of 51 windows centered at t = 183λ N days for λ N = 1, 2, • • • , 51. For each of the above FT a power trendline is calculated and superimposed over the signal. Table 1 shows the values of the slope q for both indexes where a fractional order behavior is clearly observed.

We know that the lower the negative slope the higher the attenuation of the high frequencies and, consequently the smoother the time evolution of signal. Therefore, the Fourier transform analysis suggests that the NASDAQ index is more volatile than the Dow Jones index. This is in accordance with reality. In fact, the NASDAQ index tends to have a more variable quotation -it usually outpaces the Dow Jones index both to the upside and the downside. This is because the NASDAQ index is heavily weighted to technology stocks, containing more speculative (i.e., high risk/high reward) companies than the Dow Jones index, which is constituted by large and stable companies. It is interesting to note, in the window Fourier transform, that the NASDAQ volatility has been increasing in the last years. However, this is not surprising if we have two aspects in mind: i) volatility is normally seen in the market as a proxy for the investors' emotions and fears, and ii) in the last decade we have witnessed two major crisis and, therefore, investors have been quite anxious.

Peaks in the Fourier transforms' charts for both indexes are easily identifiable. They occur to the frequencies corresponding to 1 day (ω 1 = 7.27E -5rad/sec), 1.5 days (ω 2 = 4.85E -5rad/sec) and 3 days (ω 3 = 2.42E -5rad/sec). This suggests the existence of a short-range periodicity in the stock market. In fact, ω 3 corresponds simply to the sampling frequency and ω 2 , a cross-effect between ω 1 and ω 3 . However, ω 3 reveals that we have an half-week period limit cycle.

Pseudo Phase Plane Analysis

In order to study the PPP properties of the indexes x D (t) and x N (t), we adopt the earliest time t at which the autocorrelation ρ has a point of inflection. This value is the delay time, τ , used for the PPP construction [START_REF] Lima | Pseudo phase plane, delay and fractional dynamics[END_REF].

Figures 7 and8 depict the autocorrelation ρ versus τ , for the Dow Jones and NASDAQ indexes, respectively. For the Dow Jones index the first local minimum occurs for the time lag τ D = 2330 days, while for the NASDAQ index it occurs for τ N = 2174 days.

Figures 9 and 10 depict the PPP values for the Dow Jones and NASDAQ indexes, for the chosen time lag. The PPP charts reveal that we can subdivide them into several different partitions. Based on a visual analysis of the pattern of the Dow Jones PPP chart we decided to consider five partitions as illustrated in Figure 11. For the NASDAQ index's PPP chart we consider six partitions as shown in Figure 12.

For each of these partitions, power law x(t) = ax(tτ ) b and linear x(t) = mx(tτ ) + n trendlines are calculated. Tables 3 and4 depict the values of the trendlines parameters {a, b} and {m, n} for both indexes.

For all of the PPP partitions we superimpose (Figures 11 and12), over the temporal data, the corresponding values of the power law and the linear trendlines mappings over the original data. Moreover, the corresponding errors, e P (t) = x(t)a[x(tτ )] b , and e L (t) = x(t)mx(tτ )n and histograms are obtained and, for the each error type, the corresponding values for the arithmetic mean and the standard deviation are calculated. Tables 5 and6 show the values of the arithmetic mean and the standard deviation {μ, σ} of the errors e P and e L , for the Dow Jones and NASDAQ indexes, respectively. 1 Slope q for the windowed Fourier transform, for the NASDAQ and the Dow Jones indexes. x

Fig. 11 The five partitions and corresponding power law and linear law trendlines for the Dow Jones index's PPP.

(t) = x(t) -a[x(t -τ )] b and e L (t) = x(t) -mx(t -τ ) -n.
Figures 13 and 14 depict the PPP partitions and the trendlines approximation for the Dow Jones and the NASDAQ indexes, respectively, for both types of trendline mappings. Figures 15 and16 depict the relative error's charts for the Dow Jones and the NASDAQ indexes, respectively.

These charts demonstrate that we have relationships between distinct time periods. This observation is consistent with the fractional order long-range memory effect depicted by the Fourier transform. At the present time it is neither clear what is the required number of partition to characterize completely the PPP, nor the relations between the trendline parameters and the Fourier transform details. Moreover, further research towards having the "best" trendline is needed. Nevertheless, it is clear that long term memory relations exist and that more research efforts are still needed.

Conclusions

The Dow Jones and The NASDAQ indices were studied using several techniques usually adopted in dynamical systems. The Fourier spectrum of the Dow Jones and NASDAQ indexes was approximated by trendlines. Based on the slope of the trendlines the fractional order behavior was evidenced. To provide simultaneous insight in time and frequency behavior of the signal we study the spectrum signals using a sliding-window gaussian Fourier transform.

For the PPP reconstruction an alternative technique based on the time series was also adopted. It was calculated the time lag for each index signal using the first minimum value of the autocorrelation. The tests suggest that the time lag obtained for the minimum autocorrelation values leads to patterns of relationship between several time partitions. Fig. [START_REF] Rui | A data-reconstructed fractional volatility model[END_REF] Temporal data and corresponding power law and linear law mapping for the NASDAQ index.

Figure 5

 5 Figure3shows the amplitude of the FT for the Dow Jones (|F{x D (t)}|) signal index. A trendline is calculated and it is superimposed over the signal. For the Dow Jones index the slope yields q = -0.79, revealing a fractional order behavior. Figure4shows the amplitude of the FT for the NASDAQ |F{x N (t)}| signal index. A trendline is also calculated and superimposed over the signal. For the NASDAQ index the slope yields q = -0.83 and reveals, also, a fractional order behavior. In fact, in both cases we get a fractional noise in between white and pink noise, corresponding to a considerable volatility.We now adopt the windowed Fourier transform and we consider α = 2.5, T = 366 days (one year) and that two consecutive windows are superimposed over a range of window length of β = 50%.For the Dow Jones index we considerer a total of the 110 windows centered at t = 183λ D days for λ D = 1, 2, • • • , 110 and for the NASDAQ index we adopt a total of 51 windows centered at t = 183λ N days for λ N = 1, 2, • • • , 51.Figure5depicts the amplitude of a sliding-window Fourier Transform, |Fw{x D (t)}|, centered at t = 183λ D for λ D = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100}, for the Dow Jones index, while figure6shows the amplitude of a sliding-window Fourier Transform, |Fw{x N (t)}|, centered at t = 183λ N for λ N = {5,10, 15, 20, 25, 30, 35, 40, 45, 50}.

Fig. 1

 1 Fig. 1 The temporal evolution of the daily closing value of Dow Jones index, from Oct, to Jun, 2009.

Fig. 2

 2 Fig. 2 The temporal evolution of the daily closing value of NASDAQ index, from Feb, to Jun, 2009.

Fig. 3

 3 Fig. 3 |F{x D (t)}| and the trendline for the Dow Jones signal index.

Fig. 4

 4 Fig. 4 |F{x N (t)}| and the trendline for the NASDAQ signal index.

Fig. 5

 5 Fig. 5 Windowed Fourier transform |Fw{x D (t)}| for Dow Jones index for α = 2.5, T = 366 days and β = 50% centered at t = 183λ D for λ D = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100}.

Fig. 6

 6 Fig. 6 Windowed Fourier transform |Fw{x N (t)}| for NASDAQ index for α = 2.5, T = 366 days and β = 50%, centered at t = 183λ N for λ N = {5, 10, 15, 20, 25, 30, 35, 40, 45, 50}.

Fig. 7

 7 Fig. 7 Autocorrelation ρ vs time lag τ , for the Dow Jones index.

Fig. 8 9 9Fig. 9

 899 Fig. 8 Autocorrelation ρ vs time lag τ , for the NASDAQ index.

Fig. 10

 10 Fig. 10 Pseudo phase plane for Nasdaq index and τ N = 2174.

Fig. 12

 12 Fig.[START_REF] Gabaix | A unified econophysics explanation for the power-law exponents of stock market activity[END_REF] The six partitions and corresponding power law and linear law trendlines for the NASDAQ index's PPP.

5 Fig. 13

 513 Fig.[START_REF] Francisco Martins De Lima | Análise dinâmica de vibrações em manipuladores robóticos[END_REF] Temporal data and corresponding power law and linear law mapping for the Dow Jones index.

Table 2

 2 Date of the Begin and the End time instants in each PPP partition, for the NASDAQ and the Dow Jones indexes.

		Dow Jones		NASDAQ
	Partition	Begin		End	Begin	End
	1	21/12/1932	15/09/1976	13/08/1971	29/11/1990
	2	20/09/1976	10/09/1987	30/11/1990	02/08/1991
	3	20/10/1987	16/03/1994	05/08/1991	21/04/1995
	4	17/03/1994	13/05/1998	24/04/1995	05/04/1999
	5	26/01/1988	13/03/2000	04/02/2000	13/07/2000
	6					22/04/1999	03/02/2000
			Dow Jones	NASDAQ
	Partition	a		b	a	b
	1	22.0338		0.5627	0.2823	1.415
	2	6.2494		0.8838	0.7073	1.384
	3	1308.2781	0.2482	8579377001	-2.342
	4	691.1133	0.3202	171.6524	0.3527
	5	3100.4451	0.1421	191.0417	0.2529
	6					6156.9266	-0.1193

Table 3

 3 Parameters {a, b} of the power law trendline for the PPP partitions, for the NASDAQ and the Dow Jones indexes.

Table 4

 4 Parameters {m, n} of the linear law trendline for the PPP partitions, for the NAS-DAQ and the Dow Jones indexes.

		Dow Jones	NASDAQ
	Partition	m	n	m	n
	1	0.8878	235.13	3.6635	-152.48
	2	2.3964	389.03	9.8351	-1080.1
	3	0.7549	7235.3	-7.7537	7417.8
	4	0.618	7400.4	0.5304	1449.4
	5	0.2465	8985	0.0969	1167.0
	6			-0.0875	2636

Table 5

 5 Arithmetic mean and standard deviation {μ, σ} of the error histograms for the Dow Jones index: e

		Power law, e P	Linear law, e L
	Partition	μ	σ	μ	σ
	1	17.5949	265.2153	-0.0117	276.8396
	2	1.5065	506.9526	-0.0164	509.4117
	3	3.6146	517.0247	0.0118	555.1289
	4	0.1030	115.6174	0.0379	113.8387
	5	0.0052	143.5438	0.1396	143.4068
	6	-0.0002	110.92	-0.0202	110.9426

P (t) = x(t)a[x(tτ )] b and e L (t) = x(t)mx(tτ )n.

Table 6

 6 Arithmetic mean and standard deviation {μ, σ} of the error histograms for the NASDAQ index: e P

NODY9680_source.tex; 24/02/2010; 10:54 p. 7

NODY9680_source.tex; 24/02/2010; 10:54 p.[START_REF] Francisco Martins De Lima | Análise dinâmica de vibrações em manipuladores robóticos[END_REF]