Open access peer-reviewed chapter

Appell-Gibbs Approach in Dynamics of Non-Holonomic Systems

Written By

Jiří Náprstek and Cyril Fischer

Submitted: 15 November 2017 Reviewed: 06 March 2018 Published: 18 July 2018

DOI: 10.5772/intechopen.76258

From the Edited Volume

Nonlinear Systems - Modeling, Estimation, and Stability

Edited by Mahmut Reyhanoglu

Chapter metrics overview

1,326 Chapter Downloads

View Full Metrics


Hamiltonian functional and relevant Lagrange’s equations are popular tools in the investigation of dynamic systems. Various generalizations enable to extend the class of problems concerned slightly beyond conventional limits of Hamiltonian system. This strategy is very effective, particularly concerning two-dimensional (2D) and simpler three-dimensional (3D) systems. However, the governing differential systems of most non-holonomic 3D systems suffer from inadequate complexity, when deduced using this way. Any analytical investigation of such a governing system is rather impossible and its physical interpretation can be multivalent. For easier analysis, particularly of systems with non-holonomic constraints, the Appell-Gibbs approach seems to be more effective providing more transparent governing systems. In general, the Appell-Gibbs approach follows from the Gaussian fifth form of the basic principle of dynamics. In this chapter, both Lagrangian and Appell-Gibbs procedures are shortly characterized and later their effectiveness compared on a particular dynamic system of a ball moving inside a spherical cavity under external excitation. Strengths and shortcomings of both procedures are evaluated with respect to applications.


  • Appell-Gibbs function
  • Lagrangian approach
  • non-holonomic systems
  • engineering applications

1. Introduction

The energy contained in a dynamic system is given by a scalar potential Et. It is a function of time and system response components (displacement, velocity, and acceleration vectors). Moreover, Et is a function of system parameters, position in a field of forces (potential or not), internal sources of energy and of the system evolution including a residual energy. The total energy of the system increases or decreases accordingly with external excitation and dissipation of energy. The form of energy contained within the system can have a deterministic or stochastic character and similarly also excitation and dissipation.

Considering the mechanical energy only, the total energy increase/decrease of the system with respect to time should be in equilibrium with the energy supplies and energy losses due to dissipation. This relation can be outlined by the following equilibrium:


where Pt is power supply (excitation energy per unity time) and St the specific dissipation of energy also per unity time (supposed to be independent on accelerations x¨). Functions Pt,St can dispose in special cases with a superior potential, which, however, cannot be incorporated into the potential part of total energy. Eq. (1) has a scalar character.

The energy is a primary value characterizing the system state and its evolution in time. The function Et and external influences are a background for the derivation of a governing differential system characterizing the system response with respect to initial and boundary conditions. The governing differential system is then deduced from the equivalence of Eq. (1) type using an adequate variational principle. It claims that the form of the system response corresponds with the minimum of energy spent among all admissible shapes of the system reaction. Take a note that many important settings of external forces and dissipation mechanisms do not admit the formulation by means of potentials. In such cases, they should be incorporated separately into the governing differential system using complementary principles and theorems, for example, virtual works, and so on.

We can find in monographs, for example, [1, 2, 3, 4, 5] and many others, various formulations of potentials Et and functions Pt,St combining the system parameters (physical and geometric) and the system response vectors x-displacements, ẋ-velocities, and x¨-accelerations. They can be selected in individual cases with respect to physical or geometric complexity of the system, components of the response, which are to be found, deterministic or stochastic character of the system and its excitation, and so on.


2. Basic considerations

Approaches commonly applied to construct mathematical models of dynamic systems with multiple degrees of freedoms (MDOF) follow mostly from principles symbolically outlined by Eq. (1). The equation of this type can be deduced using, for instance, a procedure of virtual displacements. They balance the energy flow in every step and subsequently applied minimization steps try to select such response trajectories, which represent a minimum of energy consumption among all admissible shapes. Let us get briefly through Lagrangian and Appell-Gibbs procedures in order to compare their basic properties. Later, we recognize that most of these properties can be regarded as positive or negative in dependence on a particular problem. Therefore, the solution method should be selected in every particular case very sensitively.

Let us remember that the aim of this chapter is a comparison of Lagrangian and Appell-Gibbs approaches effectiveness to process dynamic systems in holonomic and non-holonomic settings and to help estimate which one is more suitable to be employed in a particular case. Despite that the most important features of non-holonomic systems themselves are briefly treated as well, but for thorough evaluation of their properties, special literature should be addressed. Except five monographs cited in introductory section containing a large number of additional relevant references, a vast number of papers have been published concerning the investigation of various properties of non-holonomic systems.

Motion of an MDOF system with n degrees of freedom can be described by a system of n differential equations and l constraints:


Vector X represents external forces, while λ are unknown multipliers. The summation in Eq. (2a) characterizes influence of constraints (holonomic and non-holonomic) related with constraints (Eq. (2b)). These constraints reduce the number of the original degrees of freedom from n to k=nl. The system (Eq. (2)) includes n+l differential equations for x and λ unknown functions t, which can be determined, provided x,ẋ are given in an initial point t0. If the system (Eq. (2b)) is fully integrable, it provides l functions fr=frxt,r=1,..,l and constraints can be formulated as fr=frxt=cr. They are exclusively of a geometric character and the system is holonomic. Corresponding constraints are formulated in displacements only. In principle, l components of x can be eliminated and then remains to analyze the system with nl unknowns. Then, it can be considered λ0, and the second part on the right side of Eq. (2a) vanishes. The system with holonomic constraints takes the form:


However, frankly speaking, such an operation is possible rather exceptionally. In general, the full form of Eq. (2) should be treated, despite the system is holonomic. If some (or all) of constraints (Eq. (2b)) are not integrable, then the system is non-holonomic. In practice, we encounter these cases when the formulation of constraints includes velocities (more often velocities only).

We should remember that the non-holonomic constraints introduced in Eq. (2b) represent the most simple version of such constraints, as they are linear and in velocity. Many applications, for example, robotics, wind engineering, automotive systems, plasma physics, and so on, present more complicated types of non-holonomic constraints. Notifications to nonlinear non-holonomic constraints in velocity are given in elderly monographs [2, 3]. Later, many papers have appeared presenting results of systematic research at this field originating from particular physical or engineering problems, for example, [6, 7, 8], where higher derivatives of velocities in non-holonomic constraints are discussed. These attributes have been reflecting also in pure mathematical studies with respect to control theory and systems with a delayed feedback, see, for example, series [9, 10, 11] dealing with generalized Lagrange-d’Alembert-Poincaré equations or other studies devoted to non-holonomic reduction and related problems, see, for example, [12] and many others.

Let us realize now that the virtual work of every constraint force should vanish in the meaning as follows:


Therefore, we have with respect to Eq. (2a):


This equation holds for any arbitrary virtual displacements and represents a generalization of the principle of virtual works in statics and of the d’Alembert principle. The important issue is that it does not include any reactions of constraints. It has been well investigated in the study. For many details, see monographs, for example, [1, 3] and many others.

Let us consider that velocities in constraints Eq. (2b) are increased by virtual increments δẋ, so that they read


Deducting from Eq. (6), the initial state (Eq. (2b)) holds


Virtual increments of velocities δẋ fit into constraints requested for constraints (Eq. (4)) and, consequently, in Eq. (5) the δxs can be replaced by δẋs:


We revisit Eq. (2b) and perform differentiation with respect to t:


where d/dt represents the operator /t+i=1nẋi/xi. Considering two possible movements of the system in identical initial state and velocities in time t, but with different accelerations x¨ and x¨+δx¨, then we obtain with respect to Eq. (9):


It means that virtual accelerations δx¨ satisfy constraints requested similarly like virtual displacements or velocities following Eqs. (4) or (7). Therefore, we can write:


Some authors call relations (Eqs. (5), (8), and (11)) as the first, second, and third form of the system equation, see, for example, [3] and others.

Let us multiply each equation in the system (Eq. (2a)) by velocities ẋs. Summing them together, one obtains


which can be rewritten in the form


where T,V are kinetic and potential energy, respectively, and X˜s are forces, which cannot be included into the potential energy V. In other words, relation Eq. (13) indicates that the change of the full energy (kinetic and potential) is equivalent to power (work on velocities) of all forces X˜s, which do not contribute to the potential energy V. Relation Eq. (13) corresponds to equilibrium condition Eq. (1), where functions of excitation and dissipation Pt,St correspond with the influence of non-potential forces X˜s.


3. Lagrange’s equations

The original coordinates x should be replaced with respect to Lagrangian coordinates q. The reason is that they represent the most inherent coordinates respecting the real movement of the system and configuration of external forces. Let us write basic coordinates as functions of Lagrangian ones:


It can be easily shown


We reconsider Eq. (5) where the virtual displacement δx is replaced by


The last Eq. (16) can be modified using Eqs. (15), which implies


This equation can be rewritten now in the form:


where it has been denoted:


Inspecting the polynomial L, we recognize that it consists of the polynomial of the second and first degrees of components q̇ (coefficients are still functions of displacements q and time t) and the absolute part without any velocity components q̇. We can now assign the first part to the kinetic energy T, while the part without velocities to the potential energy V. So that L can be understood as the Lagrange function as usually defined


provided the dynamic system studied is holonomic and no constraints are applied. In such a case, all variations δqs are independent and Eq. (18) can be fulfilled only if every coefficient in curly brackets vanishes individually. Consequently, we obtain Lagrange’s equations in the form:


where Qs are generalized external forces as functions of q and t. These forces are basically linear transforms of original forces Xr, see Eq. (19). If holonomic constraints are inserted, then the number of remaining degrees of freedom is lower k=nl. Nevertheless, if there is possibility to define the system after elimination of inactive DOFs, then we can consider formally k=n again and Eq. (21) remains in force.

Let us suppose now that the system includes l non-holonomic constraints and those holonomic, which cannot be eliminated. Whatever is the reason for that, it still holds k+l=n. These constraints are described by constraints in Lagrange’s coordinates (analogous with Eq. (2b)) as follows:


This time, the variations δqs are not fully independent and only those components, which satisfy conditions:


can be regarded as independent.

In such a case, the right side of Eq. (21) should be completed:


To the system, Eq. (24) should be attached l constraints Eq. (22). So that, finally we have the system of n+l equations with unknowns q and λ. Multipliers λ are linearly related with forces in constraints. In particular cases, multipliers λ can be physically interpreted, for instance, they can have a meaning of reactions of a body moving along a given trajectory. Very knowledgeable explanation about manipulation and interpretation of Lagrange’s multipliers from the viewpoint of a general theory as well as of employment in particular cases can be found in the monograph concerning non-holonomic systems, see [2]. For additional information and a large overview of additional literature resources, see [5].

The real dynamic system is always influenced by energy dissipation. Some simple models can be introduced using Rayleigh function R, see, for example, [3]. This way is typically applicable, if linear viscous damping is considered and the Rayleigh function has a quadratic form in velocities q̇s. We can include this factor symbolically into Eq. (24), which reads now:


Hence, the completed system Eqs. (22) and (25) with n+l unknowns can be considered. However, we should be aware that this supplement is rather intuitive and does not follow from any rigorous derivation, although in practice it is widely and successfully used. Nevertheless, comparison of this system with the general relation Eq. (1) introduced in Section 1 is obvious. Take a note that more sophisticated versions of Lagrange’s equations have been developed inspired by physical problems; see, for instance, generalized Lagrange-d’Alembert-Poincaré equations discussed in [11].

Let us add that many details (internal mechanisms and inclusion into governing system) concerning more sophisticated models of damping can be found in monographs of the rational dynamics, for example, [1, 4]. See also papers oriented to practical aspects of the damping either of natural, for example, rheological, aeroelastic origin, or intentionally included in order to achieve the highest damping effectiveness, for instance [13].


4. Appell-Gibbs function and equation system

Although the Appell-Gibbs approach is not referred so often in the study as the Lagrangian procedure, there are some monographs treating the analytical dynamics, for example, [1, 3], where detailed features of this method are explained. Moreover, journal papers can be found where special aspects of the Appell-Gibbs approach are discussed. A close relation of the fifth Gaussian form and the Gibbs equations from the viewpoint of Dynamics is studied, for example, [14, 15], important remarks for application are concerned in [16, 17], as well as possibilities of extension for systems with time-dependent masses [18] are indicated.

Let us briefly outline principal steps leading to the Appell-Gibbs differential system with respect to essentials ascertained and introduced in Section 2. We should be aware that generalized external forces Qs, introduced in Eq. (19), follow in principle only k degrees of freedom, which remained free; thereafter, l constraints have been applied and the original number n of DOFs has been reduced to k=nl,0<ln. However, due to complicated relations inside the dynamic system, this fact is rather impossible to be employed in basic coordinates xs,s=1,..,n and Lagrange’s coordinates qs,s=1,..,n should be addressed, as we have also seen in previous Section 3. Nevertheless, it is worthy to involve only such coordinates qs, which correspond to k remained DOFs. It can be easily expressed in Lagrange’s coordinates, unlike basic coordinates xs. So that, as the first step, we reformulate some expressions of Section 2 concerning the transform from basic to Lagrange’s coordinates.

Velocities ẋr,r=1,..,n should be evaluated with respect to the fact that coordinates xr are functions of all Lagrange’s coordinates qs,s=1,..,k and time t, see Eq. (14):


which also implies


Differentiation of Eq. (26) with respect to time gives


The incremented acceleration vector, when keeping velocities and displacements, can be formulated as follows:


Deducting Eq. (28) from Eq. (29), one obtains


Hence, it can be written


where Qs are identical generalized forces, as they have been defined in Eq. (19). With reference to Eq. (11), we can reformulate this equation as follows:


This relation will be used later, see Eq. (36).

As a principal step of this section, we define now the Gibbs function G concentrating “acceleration energy” included in all n DOFs as follows:


When we pass from basic to Lagrange’s coordinates, only k active coordinates remain in force and so the expression Eq. (33) can be rewritten:


Expressions Eqs. (33) and (34) differ only in terms independent from accelerations.

Let us introduce the function H:


and evaluate its virtual increment:


The last parenthesis vanishes due to the relation Eq. (32). Therefore, if δx¨0, then the function δH is always positive:


which implies that accelerations q¨s,s=1,..,k should lead to a minimum of the function H, which means:


The energy dissipation terms Rx,Ry,Rz should be added to the right side of Eq. (38). At this moment, the conformity of Eq. (38) with the equivalence Eq. (1) is well pronounced, similar like in the previous section. The system Eq. (38) should be completed by geometric constraints:


Equations (38) and (39) are the Gibbs-Appell differential system including n equations, which can be written in the normal form and hence it is suitable to be immediately investigated using common methods.

The differential system (Eqs. (38) and (39)) represents the simplest and in the same time the most general form of equations of the dynamic system movement. The form of this system is very simple, and it can be used with the same effectiveness to the investigation of holonomic as well as non-holonomic systems, as the constraints can represent non-holonomic but also holonomic type of constraints. Unlike the Lagrangian approach, the non-holonomic or non-eliminable constraints do not augment the number of differential equations.

Procedure of the Appell-Gibbs equations employment in particular cases is obvious, looking back at this section. In the first step, the “so called kinetic energy of accelerations” 12r=1Nmrx¨r2 is composed using n acceleration components of the vector x¨. It represents the Appell-Gibbs function G. In a general case, this function includes also all coordinates x and velocities ẋ. Nevertheless, it is important that G in Lagrange’s coordinates contains only k selected components of accelerations q¨. Anyway, all n components of q̇ and q are still included as a result of a transformation from basic to Lagrange’s coordinates.

It is worthy to remind that the differentiation outlined in Eq. (38) is very easy in a particular case. Indeed, let us realize that G can be symbolically expressed as a sum of quadratic function of accelerations q¨s,s=1,..,kG2, linear function of these components G1 and function without accelerations G0. Differentiating G2, one obtains the relevant acceleration component in a linear form, which will be moved onto the left side together with a coefficient, which can be a function of all velocities and displacements q̇s,qs,s=1,..,n. Differentiation of G1 leads to acceleration-free coefficients and G0 can be omitted leading to zeroes. Sometimes, the so-called reduced Appell-Gibbs function G is defined where G0 is a priori omitted.

In the second step, the work of k given forces Q on k virtual displacements q is carried out. It has the form s=1kQsδqs. We substitute now back into Eqs. (38) and add l=nk geometric constraints following Eqs. (39). So we obtain k+l=n differential equations for n components of the vector q¨t. Take a note that no unknown multipliers λ emerge here, which on the other hand increases the number of unknowns in a Lagrangian approach.

The procedure working with accelerations instead with velocities provides much simpler governing differential system. Unlike velocities, the acceleration components in the Appell-Gibbs function are included only in a few parts of energy expression. Therefore, all parts including only velocity and displacement components disappear during the differentiation of the Appell-Gibbs function with respect to q¨r,r=1,..,k, and therefore they can be considered beforehand as unimportant.

Investigating problems with rotations, we work with Lagrange’s coordinates ω, which represent in fact velocities. So that by solving the abovementioned differential system, the displacements and velocities ω emerge as results. Rotations themselves remain unattended. May be, it is a forfeit for a relative simplicity of the governing system in comparison with the Lagrangian approach. However, this shortcoming is mostly apparent only. The main part of the result represents usually displacement components, which are obtained without restrictions. Together with velocities ω, they represent a full set of information needed to get through the shape of trajectories of the system response including rotation (illustrative example will be presented later in Section 6). If detailed rotations (not only velocities) are still needed, a subsequent integration can be performed independently using differential relations between rotation velocity vector ω and (for instance) Euler angles, see monographs [1, 3, 4] and others. They provide a detailed description of time history of a body orientation as a function of time t. This step can be useful, for instance, when a detailed animation is needed for presentation purposes.


5. Planar movement of a ball in a spherical cavity

5.1. Engineering motivation

Passive vibration absorbers of various types are very widely used in civil engineering. TV towers, masts, and other slender structures exposed to wind excitation are usually equipped by such devices. Conventional passive absorbers are of the pendulum type. Although they are very effective and reliable, they have several disadvantages limiting their application.

These shortcomings can be avoided using the absorber of the ball type. The basic principle comes out of a rolling movement of a metallic ball of a radius r inside of a rubber-coated cavity of a radius R>r. This system is closed in an airtight case, see, for instance, Figure 1. First papers dealing with the theory and practical aspects of ball absorbers have been published during the last decade, see [13, 19].

Figure 1.

Dynamic scheme of (a) spherical pendulum absorber, (b) ball absorber, and (c) ball absorber during testing in a dynamic laboratory, see [19].

5.2. Planar layout of the system, Lagrangian procedure

The version, when the ball is forced to move solely in a vertical plane, has been thoroughly studied using Lagrangian approach in [20, 21] and other detailed papers dealing not only with theoretical aspects but also with experimental verification in the laboratory and in situ examining absorbers installed on real structures.

The cavity is fixed to a vibrating structure. Their dynamic character is represented by a linear single degree of freedom (SDOF) system represented by a mass M. Inside of the cavity, the ball m in a vertical plane is moving, that is two degrees of freedom (TDOF) system should be investigated, as it is outlined in Figure 2. It follows from geometric relations:


where ρr=Rr. It holds for vertical or horizontal components of a displacement and velocity of the internal ball center:


Figure 2.

Basic scheme of a system.

Kinetic energy of a moving system of the ball m and the cavity M can be written in a form:


where m/κ=m+J/r2κ=5/7, while the potential energy is given by an expression:


The damping should be introduced in a form of a simple Rayleigh function:


m,M – mass of the ball m, mass of the cavity, M representing the protected structure;

J – inertia moment of the ball m;

bu,bφ – damping coefficients (logarithmic decrements, linear viscous damping);

Expressions Eqs. (42), (43), and (44) should be put into Lagrange’s equations of the second type, see Eqs. (24) or (25) and monographs, for example, [1, 3, 4] and others:


which give the governing equations of the system:


Equation (46) describes 2D movement of a ball absorber under excitation by the force Pt at any arbitrary deviation amplitudes including incidental transition through a limit cycle toward an open regime.

5.3. Illustration of some planar system features

Analysis of the governing system (Eqs. (46)) has been done in a couple of papers, for example, [20, 21]. Investigation has been carried out using the harmonic or multi-harmonic balance method, see, for example, [22, 23], respectively.

The system is auto-parametric, see, for example, [24] and other resources. Very rich overview of a theoretical basis of auto-parametric systems can be found in [25]. Expecting a single mode response, the Harmonic balance-based methods are applicable. Following approximate expressions for excitation and response can be written (cf., e.g., [22]):


Having four new variables α=αt,β=βt,γ=γt,δ=δt instead of two original unknowns φt,ζt, two additional conditions can be freely chosen:


After substituting Eqs. (47) and (48) into Eqs. (46) and substituting the sinϖ and cosφ functions by two terms of Taylor expansion, the harmonic balance procedure gives the differential system for unknown amplitudes Z=αβγδT, see, for example, [21, 23]:


System (49) for amplitudes Zt is meaningful if they are functions of a “slow time,” in other words, if their changes within one period 2π/ω are small or vanishing and individual steps of the harmonic balance operation are acceptable. The matrix M and the right-hand side vector F have the following form:



where A0=α2+β28,Aα=3α2+β28,Aβ=α2+3β28.

Let us consider stationary response of the system. In this case, the derivatives dZ/dt vanish and the right-hand side has to vanish too. Eq. (49) degenerates to the form of


Thus, to identify the stationary solutions, the zero solution points of F, depending on the excitation frequency and amplitude, should be traced. In the same time, the signum and the zero points of the Jacobian detJF have to be checked. The negative value of the Jacobian for a particular point indicates that the corresponding solution is stable, whereas when Jacobian vanishes a bifurcation could occur.

The curve Fαβγδω=0, projected into the planes ωR or ωS (for S2=γ2+δ2), forms the resonance curves known from the analysis of linear oscillators. However, correspondence of this curve to the original Eq. (46) is limited to the case of stationary response. It is necessary to remind that limits of stationarity of the response cannot be determined from properties of Eq. (52) itself. The complete Eq. (49) has to be taken into account for this purpose.

With respect to actual experiences regarding passive vibration absorbers and some interesting properties of system (46), the following reference input data have been introduced:


Utilizing Eqs. (52) and (51), the nonlinear resonance curves describing the stationary response of system (46) can be obtained. A set of such curves for excitation amplitudes p0=0.25,0.5,1,1.5,2.5 is shown in Figure 3. It is obvious for the first view the nonlinear character manifesting oneself by a dependence of a position of extreme points on an amplitude of excitation force. This effect is visible predominantly in a neighborhood of a conventional “linear” natural frequency of the absorber, although also the second natural frequency corresponding to the original natural frequency of the structure is affected. The resonance curves are typical for a system with “softening” nonlinearities.

Figure 3.

Nonlinear resonance curves describing the stationary response of the system for excitation amplitudes p0=0.25,0.5,1,1.5,2.5. Stable branches are shown as solid blue curves, unstable parts are indicated as the red dashed curves. Amplitudes, see Eq. (47), R=α2+β2 are shown in the left part of the figure, amplitudes S=γ2+δ2 are on the right.


6. Spatial version of the system, Appell-Gibbs procedure

6.1. Gibbs function

The spatial version of the ball absorber on the basis of rational dynamics has been widely investigated by authors of this chapter, see, for example, [26, 27]. Lagrangian approach and Appell-Gibbs procedures have been discussed in these papers combining analytical and numerical methods. Some important issues will be roughly outlined and for details see cited papers.

Unlike the planar version discussed in the previous section, the Appell-Gibbs approach is used to formulate the governing nonlinear differential system. The authors tried to formulate the spatial version using the Lagrangian procedure as well, see [28]. Although the governing system of the respective holonomic system has been successfully assembled, the further analysis appeared very cumbersome, and therefore, it has been given up to follow this way. Thus, the Appell-Gibbs approach is used to formulate the governing system. Its structure is much more transparent and represents a wider option of analytical-numerical investigation of detailed properties of the ball trajectories within the cavity.

With respect to Sections 2 and 4, the first step represents to construct the Appell-Gibbs function (often referred to as an energy acceleration function) defined as follows:


where M is the mass of the ball, J is central inertia moment of the ball with respect to point G, ω the angular velocity vector of the ball with respect to its center G, uG the displacement of the ball center with respect to absolute origin O, C contact point of the ball and cavity, A moving origin related with the cavity in its bottom point, see Figure 4. Coordinates x=xyz are Cartesian coordinates with origin in the point O. Hence, it holds:


where uA is the displacement of the moving origin A with respect to absolute origin O, uC the displacement of the contact point C with respect to moving origin A, un the displacement of the ball center G with respect to contact point C, n the cavity normal unit vector in point C.

Figure 4.

Ball rotation vector in moving coordinates.

Geometry of the cavity (radius R) with respect to moving origin A is given by equation:


where xA=xAyAzA are Cartesian coordinates with origin in the moving origin A.

Using Pfaff theorem and adopting a conjecture of non-sliding contact between the ball and the cavity, the respective non-holonomic constraints of “perfect” rolling can be deduced after a longer manipulation:


where ρ=1r/R.

In order to substitute for accelerations uo into the Appell function (Eq. (54)), let us differentiate constraints Eqs. (57).

Several manipulations provide expressions for components of the ball center acceleration u¨G, which consist of acceleration in the moving origin A:u¨A representing the given external kinematic excitation and acceleration related to the point A being given by an expression: ρu¨C:


Because the kinematic excitation is supposed to be horizontal, u¨Az=0 into Eqs. (58) should be substituted.

Expressions Eqs. (58) are to be substituted into Eq. (54). Thereby, we obtain the Appell-Gibbs function G for the system investigated. The function G can be significantly simplified keeping only terms including second-time derivatives u¨G and ω̇, which represent second-time derivatives of respective rotations. This step provides the reduced Appell-Gibbs function Gr. Using Gr, one can write the Appell-Gibbs differential system:


where FG is the external force vector acting in ball center G. Vector FG is determined subsequently using the virtual displacements principle. Let us introduce the quasi-coordinates φx,φy,φz where ωx=φ̇x,ωy=φ̇y,ωz=φ̇z. The only external force acting in the ball center is the gravity. Therefore, the elementary work performed can be expressed as


Virtual displacement δuGz can be determined using the third non-holonomic constraint in Eqs. (57). It holds


and therefore


At the same time, the elementary work can be expressed in terms of quasi-coordinates:


Comparing coefficients at respective virtual components δφx,δφy,δφz, in Eqs. (61) and (63), one obtains


The damping will be introduced later in Section 6.3 in order to separate energy conservative approach and enable to discuss various stationary regimes with respect to parameter and excitation settings.

6.2. Governing system

Carrying out the differentiation outlined in Eqs. (59), and respecting Eqs. (64), it can be written after some adaptations:



Ω̇s=uCxω̇x+uCyω̇y+uCzRω̇z,Js=J+mρ2R2/mρ2.mass interia moment of the ballwithrespecttocenterofthecavityE66

It can be shown that Ω̇s=0 and therefore the second column on the left side of the system Eqs. (65) should be omitted. External excitations are specified by movement or acceleration in the point A. Hence, kinematic excitation in A is given as follows:


as it can be seen in Eqs. (65). Provided we need to investigate the response processes in a vertical plane, only one component remains non-zero and the second vanishes as well.

In order to obtain the system Eqs. (65) in the form with first-time derivatives concentrated on the left side, the first derivatives u̇C in its right sides should be expressed in displacements uC using non-holonomic constraints Eqs. (57):


Therefore, we obtained the system of six non-linear ODEs (Eqs. (65) and (68)) in a normal form with six unknown functions of time: uCx,uCy,uCz,ωx,ωy,ωz. Vector uC depicts displacements of the contact point and can be used to study the movement of the ball from a global point of view. Detailed behavior of the ball as a rotating body is given by angular velocities ω. If the time history of rotation should be traced, then a subsequent run is necessary to obtain rotations by means of Euler angles as solution of the system of three ODEs with an input of angular velocities ω.

6.3. Influence of the damping

Influence of the damping will be taken into account. Basically, two sources of the energy dissipation are ruling in the system: (1) dissipation due to air dynamic resistance and (2) energy loss in contact of the cavity and rolling ball. The former one can be neglected with respect to obvious geometric configuration of the device and relative velocity ball/cavity. Concerning the latter one, complicated energy dissipating processes are ruling in contact of the ball with cavity. Nevertheless, supposing that no slipping arises in the contact, the dissipation process can be approximated as proportional to relevant components of the angular velocity vector ω and the quality of the cavity/ball contact. Considering the obvious setting, the respective material coefficients characterizing the rolling movement of the ball can be considered constant regardless of the direction in the tangential plane to the cavity in the point C, see Figure 4. The coefficient determining the rotation resistance around the normal vector n in the contact point C is different as a rule. Therefore, the resistance moment vector D can be expressed in moving coordinates p,q,n, see Figure 4, as follows:


Components of the above vector can be written in a form as follows:


where κr,κs are coefficients of “viscous resistance” of rolling and spinning. Their meaning is: the moment for a unity rotation per second, that is (Nms/rad).

Turning of the vector DG=DGxDGyDGzT expressed in xyz coordinates into the vector D can be written as


The transformation matrix TC reads


where ν2=xC2+yC2. The matrix TC is orthogonal and, therefore, the inverse transformation goes using matrix TC1=TCT, in particular:


Components of the vector DG should be incorporated onto the right side of Eqs. (59), where right sides should be completed. It means that the elementary work δFG following Eq. (60) must be completed by a negative dissipating work due to DG.


Repeating the further derivation like in Section 6.2, one can revisit the system Eqs. (65) and (68), where the right sides are completed and instead (Eqs. (65)) they read:


Terms DGx/m,DGy/m,DGz/m which are linear functions of ωx,ωy,ωz determine the viscous type of the damping, although intensity in individual coordinates is variable depending on the position of the ball within the cavity.

6.4. Ball trajectories within the fixed cavity due to initial conditions

A large program of a ball trajectory investigation within a spherical cavity has been performed using the differential system (Eqs. (68) and (75)). Basically, it consists of two groups which are briefly illustrated in this and the next subsections. The first group concerns the fixed cavity (no excitation is applied). The only source of energy introduced is given by the initial deflection of the ball from equilibrium position in the point A (“southern pole”), or in other words by non-homogeneous initial conditions.

Differential system (Eqs. (68) and (75)) admits a number of singular solutions which can serve as separating limits of zones within which regular solutions exhibit certain character of trajectory shape. Some of them can be found analytically from the differential system taking into account their special properties concerning individual response component along the trajectory as a whole or in certain points of these curves. For details, special papers should be referred. Take a note that most of them emerge when no damping is considered. The reason is that the trajectory should be quasi-periodic (or cyclic-stationary), which is impossible when damping is respected and no external energy supply is considered. Trajectories start in a certain point on a meridian into which the ball is elevated. Then, it is thrown horizontally along the cavity parallel circle. Let us mention a few of the most important:

  1. circular trajectory in horizontal plane. No initial spin is considered ωn0=0. The impulse applied corresponds with the initial velocity ω=ωps,0,0, where it holds for ωps:


This case is the most important and can be called separating circle (SC).

  1. 2. circular trajectory in inclined plane, see Figure 5, ω0=100.0,0.0,0.0. This state is exactly valid for ωp0. The space spiral type trajectory changes from SC upwards successively into the upper hemisphere. Before the limit state for ωp0 is reached, the osculation plane of the trajectory can be recognized. It rotates around the vertical axis with a descending angular velocity as far as it vanishes and osculation and operating planes coincide.

  2. 3. trajectory of “kings crown form,” see Figure 5, ω0=5.817,0.0,5.0537. Cases, when the initial spin is considered. For a special value of ωn0=5.0537 takes a shape visible in the picture. The apexes of this curve correspond to ω=0.0,0.0ωn0 and u=0.0,0.0,0.0, which is a clue to find forms and parameters of this special case. This trajectory is reached from SC, increasing the initial spin velocity until the limit value. If it is lower, the trajectory has the spiral form. For a higher value, it became a curly form, see Figure 5, ω0=5.817,0.0,10.0. The limit state for infinite initial spin represents the ball apparently fixed in the initial point and not moving neither horizontally nor vertically.

Figure 5.

Illustration of the ball trajectories; cavity is not excited; energy only supply is due to non-homogeneous initial condition; in every triplet: movement time history of the contact point C: uCx,uCy,uCz; vertical view of trajectories uCx,uCy components; axonometric view of trajectories; parameters above triplets: initial values of ω0=ωp0ωq0ωn0 and damping parameters: κr,κs; line (a): no spin, no damping, line (b): spin considered, line (c): spin and damping considered.

Let us have a look at the bottom two pictures in Figure 5. They respect the influence of the damping. Coefficients κr,κs are different as it corresponds to conditions in the real system. The left demonstrates trajectory for positive initial spin and the right for negative initial spin. The transition through limit cases mentioned earlier is visible. The trajectory obviously finishes in the bottom “southern pole” of the cavity.

6.5. Ball trajectories within kinematically excited cavity

The second group of tests deals with the cavity which is undergone to kinematic excitation in a horizontal plane (only one-direction excitation is reported here).

Two extensive series of tests demonstrate the auto-parametric character of the system. In the first series, the response has been evaluated separately for every excitation frequency ω starting from homogeneous initial conditions. Figure 6 shows some selected results of numerical simulations which follow from the differential system (Eqs. (68) and (75)). We briefly point out a couple of features visible in Figure 6. In the picture (a), we can see the maximal horizontal amplitude of the ball trajectory, when the cavity is kinematically excited in the horizontal plane in x direction. The solid curve represents max uCx and the dashed curve is max uCy as functions of the exciting frequency ω. We can see that in the interval ω02.84, the semi-trivial solution is stable and so uy=0. The point ω=2.84 is a beginning of the resonance zone, which spans in ω2.842.99, where auto-parametric resonance occurs and amplitudes of both response components are commeasurable. For ω>2.99, the semi-trivial solution is regained. Samples of the trajectory shape are plotted in picture (b) for four frequencies ω=2.84,2.88,2.92,2.96. Their vertical views demonstrate the character of the semi-trivial and the auto-parametric resonance states. Take a note that the trajectory since ω=2.94 is a simple ellipse-like curve, which does not exhibit any symptom of a chaotic process. Compare this finding with analysis concerning the sweeping up and down excitation frequency for ω around and above B2 bifurcation point (BP) (see Figure 7 and explanation later in this subsection).

Figure 6.

Response of the ball in the resonance and adjacent zones due to harmonic horizontal excitation of the cavity: (a) amplitude of the displacement as a function of the excitation frequency; (b) vertical views of the ball trajectory for frequencies ω=2.84,2.88,2.92,2.96.

Figure 7.

Amplitudes of the ball displacement under cavity harmonic excitation, when the frequency is swept up and down: (a) amplitudes overview in the interval ω1.08.0, (b) zooming in the interval ω2.43.1; curves: solid red—max uCx, dashed red—max uCy, solid black—absolute displacement amplitude, blue dashed—attraction boundary between bu1 and bl1.

The second series has been controlled by sweeping the excitation frequency up and down in a large interval and in several detailed regimes in the area of the auto-parametric resonance zone. A few of the results are visible in Figure 7. Picture (a) demonstrates amplitudes max uCx (solid curves) and max uCy (dashed curve) and the total amplitude uCr in the interval ω1.08.5. Picture (b) is the magnified detail of picture (a) within the interval ω2.803.05 in order to make visible the resonance zone.

Let us pay attention to bifurcation points (BPs). There are obviously concentrated in the resonance zone. In principle, they can be classified into two categories. The most important reveal B1 and B2. In the latter one, two branches start. The lower one bl2 approaches zero for ω which indicates the non-moving ball in the vertical view. This branch takes place in the vertical plane and basically has a form of semi-trivial solution. Its stability increases with rising ω>ωB2 as it follows from decreasing negative values of the Lyapunov exponent and of inspection of the relevant stability basins. The upper branch bu2 is spatial. It follows from the resonance zone where the spatial response type has a chaotic character. The relevant attractor reveals as an annular concentric area with diminishing width with increasing ω. The trajectory very quickly approaches a circular form in the horizontal plane. Its level with respect to the vertical axis rises and approaches “equatorial” position. However, the stability of this trajectory decreases, and we can see in Figure 7 that around ω=8.0 even numerical perturbations of the integration process can overcome the stability limit (despite very small integration step) and the response trajectory falls down to the lower branch in the point D2. Its position is not fixed. If hypothetically zero perturbation occurs, it could shift to infinity and approach together with the branch bu2 the asymptote at the level R=1. Observing black max uCx and red max uCy parts of bu2, we can see that they are getting coincide with increasing ω. It means that trajectory approaches the circle with radius R=1.

Let us briefly discuss the shape of the response amplitudes for ω below BP B1 and B4. The BP B1 is reached sweeping up along the branch bl1, when it loses planar character passing through B4. In such a case, the spatial response type emerges, exhibiting a chaotic response since B1. This fact is obvious also looking at the dashed red curve representing uCy, which is trivial as far as B4 and can bring the system from the semi-trivial solution into the auto-parametric resonance starting B4. Take a note that passing BP B4, planar response can remain in force, if any perturbation is avoided. It meets in B3 the branch bl2 following also a planar path for swept up ω. Branch bu1 starts in B1. Its stability rapidly decreases with descending ω. The point D1 illustrates its limited extent in sub-resonance zone. This feature is visible observing the curve bs, which represents a limit separating the area of attraction to bu1 and to bl1. Take a note that bs starts in B4 and approaches D1, despite hypothetically it goes together with D1 as far as the vertical axis in the point R=1. The bs can be earned from stability basins for ω in the adequate interval for initial value ωp=0. It corresponds to amplitude of uCx as a testing value for decision about affiliation to bu1 or bl1 attractiveness.

The interval between B1 and B2 includes the spatial response, see non-trivial amplitude max uCy. The spatial response has a chaotic character, as it has been already outlined in the previous paragraph, when commenting the branch bu2.


7. Conclusion

The common physical origin of Lagrangian and Appell-Gibbs approaches has been shown. It originates from the equilibrium of energy-level evolution in time on one side and power supply together with energy dissipation on the other side. Various formulations of this principle lead finally to different variational principles, although they follow from the same minimization of the energy spent to system response portrait. Comparing individual sections of the chapter, we can see that each one of commonly used procedures based on particular energy formulations is preferable for a certain type of problems. It can be concluded that there does not exist a single universal approach which should be recommended.

Some detailed properties of both approaches have been demonstrated in Sections 5 and 6. Both of them discuss non-holonomic problem of the ball movement within the spherical cavity under external excitation. The former one deals with a simple planar problem and shows that the Lagrangian approach is easily applicable to obtain reasonable results as far as a wide parametric discussion, which enable to earn a detailed insight into the system dynamic properties. The latter alternative represents the full space problem with six DOFs and three non-holonomic constraints. Some earlier studies tried to formulate this problem also in Lagrangian style using Lagrange’s multipliers. Finally, it proved that the relevant governing differential system is too complex and does not enable appropriate detailed analysis of dynamic properties of the system. Therefore, the space problem outlined in Section 6 has been formulated using Appell-Gibbs approach. Transparent results have been obtained as needed for practical purposes in a device design and in further study of multi-body system dynamics. Take a note regarding the classification of singular solutions and their applicability for detailed analysis, stability of various regimes of the system under kinematic excitation, transitions among semi-trivial, auto-parametric, chaotic, and other states typical for nonlinear system. Let us add that both 2D and 3D problems have been investigated respecting the full nonlinearity without any simplifications of transcendent functions and thus enabling to study all effects without any limitations in amplitudes.

A certain shortcomings which apparently follow from the knowledge of rotation velocities only (no rotations themselves are calculated) can be disregarded, when displacements have been obtained. The rotation velocities represent mostly satisfactory information. Nevertheless, if rotations are still needed, there exist several variants of a simple differential system (following rotation vector definition) relating velocity and rotation vector components. This system can be subsequently easily solved, when necessary. A hidden complexity of the Lagrangian approach follows from an implicit connection of both parts, which are independent when Appell-Gibbs procedure is applied.



The kind support of the Czech Science Foundation project No. 17-26353 J and of the RVO 68378297 institutional support are gratefully acknowledged.


  1. 1. Lurie A. Analytical Mechanics (in Russian). Moscow: GIFML; 1961
  2. 2. Nejmark JI, Fufaev NA. Dynamics of nonholonomic systems (in Russian). Moscow: Nauka; 1967. English translation—A.M.S. Translations of Mathematical Monographs (Book) Vol. 33. Providence: AMS; 1972
  3. 3. Pars LA. A treatise on analytical dynamics. 2nd ed. Woodbridge: Ox Bow Press; 1972
  4. 4. Hamel G. Theoretische Mechanik. Berlin: Springer; 1978. DOI: 10.1007/978-3-642-88463-4
  5. 5. Bloch AM. Nonholonomic Mechanics and Control. New York: Springer; 2003. DOI: 10.1007/b97376
  6. 6. Rubio Hervas J, Reyhanoglu M. Controllability and stabilizability of a class of systems with higher-order nonholonomic constraints. IFAC Journal Automatica. 2015;54:229-234. DOI: 10.1016/j.automatica.2015.02.006
  7. 7. Bloch AM, Reyhanoglu M, McClamroch NH. Control and stabilization of nonholonomic dynamic systems. IEEE Transactions on Automatic Control. 1992;37(11):1746-1757. DOI: 10.1109/9.173144
  8. 8. Reyhanoglu M, van der Schaft AJ, McClamroch NH, Kolmanovsky I. Dynamics and control of a class of underactuated mechanical systems. IEEE Transactions on Automatic Control. 1999;44(9):1663-1671. DOI: 10.1109/9.788533
  9. 9. Cendra H, Ferraro S. A nonholonomic approach to isoparallel problems and some applications. Dynamic Systems. 2006;21(4):409-437. DOI: 10.1080/14689360600734112
  10. 10. Cendra H, Grillo S. Generalized nonholonomic mechanics, servomechanisms and related brackets. AIP, Journal of Mathematical Physics. 2006;47:2209-2238. DOI: 10.1063/1.2165797
  11. 11. Cendra H, Ferraro S, Grillo S. Lagrangian reduction of generalized nonholonomic systems. Journal of Geometry and Physics. 2008;58:1271-1290. DOI: 10.1016/j.geomphys.2008.05.002
  12. 12. Bates L, Sniatycki J. Nonholonomic reduction. Reports on Mathematical Physics. 1993;32:99-115. DOI: 10.1016/0034-4877(93)90073-N
  13. 13. Pirner M. Dissipation of kinetic energy of large-span bridges. Acta Technica CSAV. 1994;39:407-418
  14. 14. Lewis AD. The geometry of the Gibbs-Appell equations and gauss’ principle of least constraints. Reports on Mathematical Physics. 1996;38:11-28. DOI: 10.1016/0034-4877(96)87675-0
  15. 15. Udwadia F, Kalaba R. The explicit Gibbs-Appell equation and generalized inverse forms. Quarterly of Applied Mathematics. 1998;56:277-288
  16. 16. Desloge E. The Gibbs-Appell equation of motion. American Journal of Physics. 1988;56:841-846. DOI: 10.1119/1.15463
  17. 17. Emami M, Zohoor H, Sohrabpour S. Solving high order nonholonomic systems using Gibbs-Appell method. In: Balan V, editor. Proceedings of the International Conference Differential Geometry and Dynamical Systems. Mangalia: Geometry Balkan Press; 2009. pp. 70-79
  18. 18. Yong-Fen Q. Gibbs-Appells equations of variable mass nonlinear nonholonomic mechanical systems. Applied Mathematics and Mechanics. 1990;11:973-983. DOI: 10.1007/BF02115681
  19. 19. Pirner M, Fischer O. The development of a ball vibration absorber for the use on towers. Journal of the International Association for Shell and Spatial Structures. 2000;41(2):91-99
  20. 20. Náprstek J, Pirner M. Non-linear behaviour and dynamic stability of a vibration spherical absorber. In: Smyth A et al., editors. Proceedings of the 15th ASCE Engineering Mechanics Division Conference. New York: Columbia University; 2002. CD #150. 10 p
  21. 21. Náprstek J, Fischer C, Pirner M, Fischer O. Non-linear dynamic behaviour of a ball vibration absorber. In: Papadrakakis M, Fragiadakis M, Plevris V, editors. Proceedings of the 3rd International Conference on Computational Methods in Structural Dynamics and Earthquake Engineering. Corfu: COMPDYN; 2011. ID #180. 14 p
  22. 22. Xu Z, Cheung YK. Averaging method using generalized harmonic functions for strongly non-linear oscillators. Journal of Sound and Vibration. 1994;174(4):563-576. DOI: 10.1006/jsvi.1994.1294
  23. 23. Ren Y, Beards CF. A new receptance-based perturbative multi-harmonic balance method for the calculation of the steady state response of non-linear systems. Journal of Sound and Vibration. 1994;172(5):593-604. DOI: 10.1006/jsvi.1994.1201
  24. 24. Haxton RS, Barr A. The autoparametric vibration absorber. Journal of Engineering for Industry. 1972;94:119-125. DOI: 10.1115/1.3428100
  25. 25. Tondl A. Quenching of self-excited vibrations. Prague: Academia; 1991
  26. 26. Náprstek J, Fischer C. Dynamic behavior and stability of a ball rolling inside a spherical surface under external excitation. In: Zingoni A, editor. Proceedings of the 6th International Conference on Structural Engineering, Mechanics and Computation, SEMC. London:Taylor & Francis; 2016. pp. 214-219
  27. 27. Náprstek J, Fischer C. Non-holonomic dynamics of a ball moving inside a spherical cavity. In: Vestroni F et al., editors. Proceedings of the EURODYN 2017—10th International Conference on Structural Dynamics. Procedia Engineering. 2017;199:613-618. DOI: 10.1016/j.proeng.2017.09.105
  28. 28. Náprstek J, Fischer C. Dynamic response of a heavy ball rolling inside a spherical dish under external excitation. In: Zolotarev I, editor. Proceedings of the Engineering Mechanics. Svratka. Prague: IT ASCR; 2013. CD #46. pp. 96-106

Written By

Jiří Náprstek and Cyril Fischer

Submitted: 15 November 2017 Reviewed: 06 March 2018 Published: 18 July 2018