Experience of development and usage of packages of symbolic computations

intended for investigation of mechanical systems.

Banshchikov Andrej V., Bourlakova Larissa A., Ivanova Galina N.,

Irtegov Valentin D., Novichov Mikhal A., Titorenko Tatyana N.

Institute of Systems Dynamics and Control Theory SB RAS

134, Lermontov str., Irkutsk, 664033, Russia

E-mail: irteg@icc.ru

Full paper in compressed Postscript *.ps.gz

In the report long-term experience of the authors on development of algorithms and automation of a research of complicated mechanical and controlled systems is considered. We have created specialized systems of computer algebra and software packages "Dynamics", "Mechanic", "Normalization" possibility and the algorithms of which are described in [1-5 and other]. The base of the algorithms, which are realized in these packages, was formed by classical methods of analytical mechanics and stability theory.

The indicated software allows after geometric description of a system of rigid and deformable bodies to make their finite-dimensional models as the Lagrange equations of the second kind (in generalized coordinate) and Euler - Lagrange (in quasi-coordinates) ones and also will obtain the linear equations in a neighbourhood of steady-state motion. The application packages differ from another software by that they contain large units, permitting to carry out in a symbolic form qualitative researches of the differential equations : to build the first integrals, linear and quadratic in generalized velocities (quasi-velocities) ; for linear (ignorable) and quadratic integrals to determine invariant manifolds of steady motions and to investigate their stability ; to fulfill normalization in Poincaré sense in a neighbourhood of singular points of the equations.

While solving of these problems the computer algebra system (CAS) fulfills : operations of matrix - vectorial algebra, partial differentiation, expansion of functions in to Taylor series, selection of main diagonal minors, calculation of determinant, evaluation of conditions of a property of having fixed sign of the quadratic forms, replacement and substitution, removal of brackets, collecting terms, trigonometric transformations, symbolic-numerical interface and other. The algorithms of analytical transformations are simplified considering the specificity of the expressions of data domain.

The developed software allows to automatize, and consequently, essentially to speed up processes of modelling and dynamic analysis of complicated systems, to avoid errors at all stages of researches.

At the moment we develop the software "Analysis" for IBM PC/AT on basis of system "Mathematica" [6,7]. In the report we consider some problems and difficulties, connected to application of modern systems of computer algebra for solving mentioned above and other problems of common mechanics and stability theory.

For the package "Analysis" new algorithms for a research in symbolic-numerical form of stability and stabilization of mechanical systems under an operation of gyroscopic, dissipative, conservative and nonconservative positionary forces are developed [8]. Experience, accumulated at creation and maintenance of systems for the analysis of mechanical systems, was transferred by authors to electrical and electromechanical systems. The following operations are fulfilled for these systems :

- the graph of a circuit is analyzed ; - Routh or Hamilton function is created ; - the equations of motion are written in the Lagrange, Hamilton or Routh form ; - the first integrals are created ; - linear in velocities integrals are reduced to ignorable ones.

Structure of the dynamic analysis of systems

Figure 1.

Our experience with symbolic computation packages allows us to conclude that CAS are perspective tool for researches in the field of theoretical mechanics. Here we will not concentrate on specific software implementations. Instead we will give a brief description of the functionality of mentioned software packages ( Fig. 1 ).

1. Description of the models.

1.1. Model of the mechanical system. The mechanical system is a system of bodies $\ S_{j}\ \ (j = 1,\ldots,N)\ ,$connected by one-two-three-degree joints, i.e. for every $\ S_{i}\$there exists the point $\ O_i\ \ (\ O_{i}\in S_{j}\$ ¨ $\ O_{i}\in S_{i}\ \ (\ i,\ j \in \{1,\ldots ,N\},\ \ i \neq j\ )\$, or joints allowing translational displacements $\ S_{i}\$ relative to $\ S_{j}\$ ( Fig. 2 ).

Figure 2.

The body $\ S_{1}\$ is a carrier, $\ S_{i}\ \ (i=2,\ldots,N)\$ is carried. Let us introduce the systems of coordinates : a system $\ O_{i}\Sigma^{i}\$ is connected with $\ S_{i}\$ (with hardening body which will be received if all deformations would be equal to zero); a system $\ O_{i}\tilde{\Sigma }^{j}\$ is connected with $\ S_{i}\$ in a small vicinity of point $\ O_{i}\$ and such that if $\ S_{j}\$ is an absolute rigid body then $\ \Sigma^{j}\$ and $\ \tilde{\Sigma }^{j}\$ are parallel. The vocation of $\ S_{1}\$ in the inertial space $\ \Sigma^{in}\$(or in some system of coordinates $\ \Sigma^{0}\$ with given motion) assigns by coordinates of point $\ O_{1}\$ and by matrix of rotation $\ \alpha^{\ 01}\$ of axes $\ \Sigma^{1}\$. The location of body $\ S_{i}\$ assigns by the matrix of rotation $\ \alpha^{\ ji}\ (\ \tilde{\alpha}^{\ ji}\ )\$ with reference to $\ \Sigma^{j}\ (\ \tilde{\Sigma}^{j}\ )\$ or by coordinates of pole $\ O_{i}\$. It is assumed that for non-rigid body $\ S_{i}\$ in every point $\ M^{\prime}\$ of body $\ (\ \mbox {\boldmath$\rho$ }(x_{1},x_{2},x_{3}) = \overline{OM^{\prime}}\ )\$the deviation reflected due to deformation is represented by a vector

$${\bf u} = \sum_{\alpha=1}^{n^{i}}\ q_{\alpha}(t)\ {\bf U}^{\alpha }(\ x_{1},\ x_{2},\ x_{3}\ ) \eqno (1.1.1)$$

where $\ q_{\alpha }(t)\$ are generalized coordinates of deformations ; $\ {\bf U}^{\alpha }(x_{1},\ x_{2},\ x_{3})\$ are eigenforms of oscillations or their functions. One defines the choice of $\ q_{\alpha }(t)\$and number $\ n^{i}\$.

1.1.1. The kinetic energy of a system. The kinetic energy of a system is calculated as a sum of kinetic energies $\ T_{i}\$ for every body $\ S_{i}\$, connected with $\ S_{j}\$:

$$T_{i}= {1\over 2}\ M_{i}\ {\bf v}^{2}_{o_i} + {1\over 2}\ \mb... ...times \mbox{\boldmath $\omega$}_{i}\ ]\cdot {\bf r}^{i}_{c_i}+$$

$$+ k_{i}\ M_{i}\ {\bf v}_{o_i}\cdot \sum_{\alpha=1}^{n^{i}} {\... ... }\ \dot q^{(i)}_{\alpha }\ \dot q^{(i)}_{\beta } \eqno (1.1.2)$$

Here $\ M_{i}\$ is the mass of the body ; $\ k_{i}= 0\ ,$ if body is rigid ; $\ k_{i}= 1\ ,$ if body is non-rigid; $\ {\bf v}_{o_i}\$ is absolute velocity of point $\ O_{i}\$:

$${\bf v}_{o_i} = {\bf v}_{o_j} + [\ \mbox{\boldmath $\omega$}_... ..._i}\ \dot q^{(j)}_{\alpha }+ {\bf v}^{j}_{o_i}\ ; \eqno (1.1.3)$$

$\mbox{\boldmath$\omega$ }_{i}\$ is absolute angular velocity of body :

$$\mbox{\boldmath $\omega$}_{i} = \mbox{\boldmath $\omega$}_{j}... ...tilde{k}\ \mbox{\boldmath $\omega$}^{i}_{j\ def} \eqno (1.1.4)$$

$\mbox{\boldmath$\omega$ }^{i}_{j\ def}\$ is angular velocity of rotation $\ \tilde{\Sigma }^{j}\$ relative to $\ \Sigma^{j}\$;

$$\mbox{\boldmath $\omega$}^{i}_{j\ def}={1\over 2} \sum_{\alph... ...f U}^{\alpha j})_{o_i}\ \dot q^{(j)}_{\alpha }\ ; \eqno (1.1.5)$$

$\mbox{\boldmath$\omega$ }^{j}_{i}\$ is angular velocity of motion $\ \Sigma^{i}$ relative to $\ \Sigma^{j} \ \ (\tilde{k} = 0 )\$ or $\ \tilde{\Sigma }^{j}\ \ (\tilde{k} = 1)\$;

${\bf v}^{j}_{o_i}\$ is the relative velocity of point $\ O_{i}\$, if $\ S_{i}\$ moves translationally (or freely) relative to $\ S_{j}\$;

$${\bf r}^{j}_{o_i}=\mbox {\boldmath $\rho$}^{\ j}_{o_i} + k_{j... ...}} {\bf U}^{\alpha i}_{o_i}\ q^{(j)}_{\alpha }\ ; \eqno (1.1.6)$$

${\bf r}^{i}_{c_i}\$ is the radius-vector of centre of mass of the body :

$${\bf r}^{i}_{c_i}=\mbox {\boldmath $\rho$}^{\ j}_{c_i} + k_{i... ...i)}\ \int_{M_{i}}{\bf U}^{\alpha i}dm\ )/M_{i}\ ; \eqno (1.1.7)$$

$\Theta^{o_i}\$ is the inertia tensor of the body relative to $\ O_{i}\$:

$$\Theta^{o_i}= \Theta^{o_{i}}_{*} + 2k_{i}\ \sum_{\alpha=1}^{n... ...a \beta }\ q^{(i)}_{\alpha }\ q^{(i)}_{\beta }\ ; \eqno (1.1.8)$$

$\Theta^{o_{i}}_{*}\$ is the inertia tensor of rigid (solidified) body,

$$\Lambda ^{o_{i}}_{\alpha} = \left( \begin{array}{ccc} \int\li... ...}_{2}^{\alpha}+x {\bf U}_{1}^{\alpha} ) dm \end{array} \right)$$

$$\eqno (1.1.9)$$

$$Q^{o_i}_{\alpha \beta } = \left( \begin{array}{lll} \int\limi... ...f U}_{1}^{\alpha} {\bf U}_{1}^{\beta} ) dm \end{array} \right)$$

$$\eqno (1.1.10)$$

$${\bf B}^{i}_{\alpha } = \int\limits_{M_{i}}\ [\ \mbox{\boldma... ...f U}^{\beta } \times {\bf U}^{\alpha }\ ]\ dm\ ; \eqno (1.1.11)$$

$$A^{i}_{\alpha \beta } = \int\limits_{M_{i}}\ (\ {\bf U}^{\alpha i} \cdot {\bf U}^{\beta i}\ )\ dm\ ; \eqno (1.1.12)$$

1.1.2. The force function of a system. Let the system is under Newton's gravity ( approximate ) to the motionless centre $\ O^{\prime}\$. The force function of a system equals a sum of force functions $\ {\cal U}_{i}\$ of each body $\ S_{i}\ \ (i=1,\ldots,N)\$:

$${\cal U}_{i}= {\nu M_{i}\over r_{o_1}}- {\nu M_{i}\over r^{3}... ...{2}\over 2} + ({\bf r}^{1}_{o_i}\cdot {\bf r}^{i}_{c_i})\ ]\ +$$

$$+\ {3\nu M_{i}\over r^{5}_{o_1}}\ [\ {({\bf r}^{0}_{o_1}\cdo... ...}_{o_i})\cdot ({\bf r}^{0}_{o_1}\cdot {\bf r}^{i}_{c_i})\ ]\ +$$

$$+ {\nu \over 2r^{3}_{o_1}} (\ \Theta^{o_i}_{11} + \Theta^{o_i... ...3}_{o_1}}(\beta^{i})^{T}\Theta^{o_i}\beta^{i}\ . \eqno (1.1.13)$$

Here $\ \nu \$ is the constant of gravitation ; $\ r_{o_1}\$ is the distance $\ O_{1}\$ from centre $\ O^{\prime}\$; $\ \beta^{i} = \alpha^{1i}\beta^{1}\ ,\ \ \beta^{1}= \alpha^{01}\gamma\ , \ \ \beta^{i}\$ is matrix $\ 3\times1\$ of cosines of angles between $\ {\bf r}^{0}_{o_1}\$ and axes $\ \Sigma^{i}$; $\ \gamma\$ is matrix $\ 3\times1\$ of cosines of angles between $\ {\bf r}^{0}_{o_1}\$ and axes $\ \Sigma^{0}\ ;\ \ \Theta^{o_i}_{jj}\$ are elements of inertia tensor $\ \Theta^{o_i}\ (1.8)\ ;\ \ \Sigma^{i}\$ are the principal axes of deformation.

The force function of a body $\ S_{i}\$ in field of constant gravity is calculated under the formula :

$${\cal U}_{i} = M_{i}\left({\bf g}\cdot {\bf r}^{o}_{c_i}\right)+ const \ . \eqno (1.1.14)$$

Here $\ {\bf {g}}\$ is vector of acceleration of constant force of gravity ; $\ \ {\bf r}^{0}_{c_i}={\bf r}^{0}_{o_1}+{\bf r}^{1}_{o_i}+{\bf r}^{i}_{c_i}\ .$

If body $\ S_{i}\$ is non-rigid, then the force function of deformation is assumed as the quadratic form of the generalized coordinates of deformation :

$${\cal U}^{def}_{i} = -{1 \over 2}\sum^{n^{i}}_{\alpha ,\ \bet... ...pha \beta }\ q^{(i)}_{\alpha }\ q^{(i)}_{\beta } \eqno (1.1.15)$$

1.2. Description of an electric circuit. Electromechanical analogies [10] allow to use the same apparatus - methods of analytical mechanics - for describing and investigation of both mechanical systems and electric circuits.

Consider a linear electric circuit, in which resistors (R), inductors (L), capacitors (C), power sources of current (I) and voltage (E) are interconnected arbitrarily. For the purpose of describing the circuit let us choose a set of independent variables characterizing its state at any time moment. Such variables may be represented by either currents in the loops or voltages in the nodes or else currents and nodal voltages. Selection of a set of variables defines the structure of the equations.

Consider an algorithm of constructing the Lagrangian for the electric circuit. A set of currents will be used as the state variables. The computations are conducted in the two stages. On the 1st stage, the set of state variables of the electrical circuit is determined by the method of loop currents; on the 2d - the Lagrangian is computed.

The general idea of the method of loop currents consists in decomposing the electrical circuit into independent geometric loops and assigning the current in each of the loops. This problem may be solved via finding a set of fundamental cycles in a graph.

Let us use the graph of the electrical circuit ( Fig. 3 ) for finding the fundamental cycles.

Let us first arbitrarily number graph's vertices and then assign the order of its path-tracing: we intend to trace graph's vertices in the direction of growth of their numbers (or vice versa). The computer implementation of the graph has the form a list of pairs of vertices corresponding to its edges.

Let $\ {\it C_e}\$ be a list of constructed fundamental cycles, S contains a list of contiguous vertices for each of the graph's vertex, D a list "free" graph's edges: $\ \{n_i,\ n_j\} \in D\$, but $\ \{n_i,\ n_j\} \not \in C_e\$ ( $\ n_i,\ n_j\ -$numbers of the graph's vertices, $\ \{n_i,\ n_j\}\ -$ a graph's edge ). The list elements D, S are ordered in accordance with the technique of path-tracing the graph. Let $\ D \neq$ Ø. Let us describe the constructing of a new fundamental cycle.

Choose the 1st element $\ \{n_0, n_1 \}\$ from D and remove it from the list. Using S find the vertex $\ n_2\$ contiguous with $\ n_1$. The vertex $\ n_2\$ is called used and is removed from S. In accordance with the same rule find the vertices $\ n_3,\ n_4$, etc. Each new chosen vertex shall not coincide with those already chosen. The path traced is stored in the form of a list of pairs $\ \{n_i, n_j\}\$.

If on the ${\it k}-$ th step it appears that there is not any unused vertex contiguous with $\ n_k$, we return to the vertex $\ n_{k-1}\$ and continue the process. When $\ n_k=n_0\$ the construction of the cycle is over. The latter is always possible since the graph corresponds to the closed electrical circuit. The constructed cycle is stored in the list $\ {\it C_e}\$. All the elements included in the new cycle are removed from D . If $\ D \neq$ Øthen we select a new element from it and the process of constructing the cycle is repeated.

It may happen that the number of loops constructed by such algorithm is less than that required by the method of loop currents. In this case the list $\ C_e$ is complemented up to the necessary number under the condition that the new constructed loop does not coincide with any of existing loops.

When denoting by $\ \dot q_{i}\$ the current in the i- th loop and considering it positive in the edge $\ \{n_k,\ n_m\}\$, when $\ n_k < n_m\$ and otherwise negative, we shall find the magnitude of the current in each branch as the sum or the difference of the corresponding loop currents.

According to [10] the kinetic energy $\ T$, the force function $\ { \cal U}$, the Rayleigh function $\ {\tilde R}$ for the linear electrical circuit has the form :

$$T = \frac{1}{2} \sum^{k}_{i,j=1} L_{ij} \dot q_{i} \dot q_{j}... ...= \frac{1}{2} \sum^{k}_{i,j=1} R_{ij} \dot q_{i} \dot q_{j}\ ,$$

where $\ L_{ij}\ -$ inductance, $\ C_{ij}\ -$ capacitance, $\ R_{ij}\ -$ resistance, $\ q_i\ -$ charge.

For constructing the function characterizing the state of the nonlinear electric circuit we use the algorithm [11,12].

2. Equations of motion.

After an evaluation of kinetic energy and force functions we discover characteristic function of the system and create the equations of motion.

2.1. Lagrange's equations. Let characteristic function be a Lagrangian $\ L=T + \ {\cal U} \$, which is shown in the form:

$$L =T_{2}+T_{1}+T_{0}+{\cal U}= {1\over 2} \sum^{k}_{i=1}\sum^... ...k}_{i=1}\ A_{i}(q,t)\ \dot q_{i} + A_{0}(q,t)\ , \eqno (2.1.1)$$

where $\ k\$ is number of degrees of freedom of a system ; $\ A_ {0}(q, t) = T_ {0} (q, t) + \ {\cal U} (q) \ ;\ \ T_ {0} (q, t) \$is part of a kinetic energy, not depending on velocities ; $\ A_ {ij} = A_ {ji} \$. Then the equations of motion are created as the Lagrange equations of the second kind :

$$\sum^{k}_{r=1}\ A_{ri}\ \ddot q_{r} + {1\over 2} \sum^{k}_{r=... ...j}} - {\partial A_{j}\over \partial q_{i}}\right)\ \dot q_{j}+$$

$$+ \sum^{k}_{r=1}\ {\partial A_{ri}\over \partial t}\ \dot q_{... ...l q_{i}} = Q_{i}\ ,\ \ \ \ (\ i=1,\ldots,k\ )\ , \eqno (2.1.2)$$

$Q_ {i} \$ are generalized nonpotential forces.

2.2. Euler-Lagrange's equations. Let the configuration of a system be set by independent parameters $\ q_ {1}\ ,\ldots,\ q_ {k} \$and quasi-velocities - linear forms on generalized velocities $\ \dot q_{i}\$ (in common case inhomogeneous):

$$\Omega_{s} = \sum^{k + 1}_{i=1} \ \tilde {a}_{si} (q) \ \dot q_{i} \ , \ \ \ \ (\ s=1,\ldots, k + 1\ ) \ , \eqno (2.2.1)$$

here $\ \Omega_{k + 1} = \dot q_ {k + 1} = 1\ \ \ \ (\ \tilde {a}_{k + 1, i} (q) = \delta_ {k + 1, i} \ {\rm\ is\ Kronecker's\ symbol\ ;} \ \ q_{k + 1} = t\ ) \$.

Jt is supposed, that $\ \tilde {a} =\ \parallel \tilde {a}_{si} (q) \parallel\$ is nondegenerate. The linear forms of differentials of generalized coordinates

$$d\pi _{s}= \sum^{k+1}_{i=1}\ \tilde{a}_{si}(q)\ dq_{i}\ ,\ \ \ \ (\ s = 1,\ldots,k+1\ )\ ,$$

determine differentials of quasi-coordinates.

The kinetic energy of system $\ T$ in quasi-velocities in common case looks like :

$$T={1\over 2} \sum^{k}_{s=1}\sum^{k}_{i=1}\ A^{*}_{si}(q)\ \Om... ...q)\ \Omega _{s} + A^{*}_{0}(q)\ =\ T_2+T_1+T_0\ . \eqno (2.2.2)$$

We shall write out an explicit aspect of the Euler-Lagrange's equations for the case, when the kinetic energy does not contain time $\ t\$:

$$\sum^{k}_{i=1}\ A^{*}_{si}\ \dot \Omega_{i} + \sum^{k+1}_{m=1... ..._{m=1}\sum^{k}_{r=1}\ \gamma^{r}_{ms}\ A^{*}_{r}\ \Omega_{m} +$$

$$+ \sum^{k}_{m=1}\sum^{k}_{r=1}\ ( {\partial A^{*}_{sm}\over \... ...s}} - {\partial {\cal U}\over \partial \pi_{s}} = Q^{*}_{s}\ ,$$

$$\dot q_{s}= \sum^{k+1}_{r=1}\ b_{sr}\ \Omega_{r}\ ,\ \ \ \ (\ s=1,\ldots, k\ )\ ,\eqno (2.2.3)$$

where $\ Q^{*}_{s} = \sum^{k}_{r=1} \ b_{rs}\ Q_ {r} \$is a generalized nonpotential force, referred to a quasi-coordinate $\ \pi_{s}\$;

$$\gamma^{s}_{nm}= \sum^{k}_{r=1} \sum^{k}_{i=1}\ (\ {\partial ... ...rn}\ b_{im}\ ,\ \ \ (\ m,s=1,\ldots,k\ ;\ \ n= 1,\ldots,k+1\ )$$

are three-index Bolzman-Gamel's symbols, $(\ b_{ij}$ are the elements of matrix $\ b=\tilde {a}^{-1} \ ).$Reducing the Euler-Lagrange's equations is automated also at presence of nonholonomic connections and surplus generalized coordinates.

2.3. Hamilton's equations. Let characteristic function be Hamilton's function

$$\ H (t, q, p) = \sum_{i=1}^{k} \ p_ {i} \ \dot q_ {i} \ -\ L\ ,$$

where $\ \dot q_{i}\$ is the function of coordinates and impulses $\ p_ {j} =\ \partial L (q, \dot q) / \partial \dot q_ {i} \$. Then differential equations of motion look like :

$${\dot q_{i}} = {\partial H(q,p) \over \partial p_{i}}\ \ ,\ \... ...artial q_{i}} + Q_{i}\ , \ \ \ \ \ (i=1,\ldots,k) \eqno(2.3.1)$$

where $Q_ {i} \$ are generalized nonpotential forces.

In electric circuits it is often the case that

$$\det \Vert {\partial^2 L(q, \dot q) \over \partial \dot q_{i}... ...rt = 0\ , \hskip 2cm (\ i=1,\ldots,k\ ; \ \ j=1,\ldots,k\ )\ .$$

Therefore it is not always possible to execute standard transformation of generalized velocities to Hamilton variable. Partial transformation is in this case carried out. The singularities of a conclusion of the Hamilton's equations for such situations are considered in [13].

3. Qualitative investigations. Principal formulas and algorithms.

3.1. Investigations in the first approximation. The differential equations obtained after their expanding in series in the neighbourhood of the given or found steady-state solution may be represented in the form

$$A \ddot x\ + \ B\dot x \ + \ C x = Q(x,\dot x)\ , \eqno(3.1.1)$$

where $\ x\ ,\ \dot x\ ,\ \ddot x\ -$ column matrix of generalized coordinates, velocities, accelerations, respectively, $\ Q(x,\dot x)\ -$ summands nonlinear with respect to $\ x, \ \dot x\$, $\ Q(0,0)=0\$. When $\ Q(x,\dot x)=0\$, we have the equations of the first approximation.

Now construct the characteristic equation (CE):

$$det(A \Lambda^{2} + B \Lambda + C )= h_{0}\Lambda^{2k} + h_{1}\Lambda^{2k-1} + \ldots + h_{2k}= 0\ . \eqno (3.1.2)$$

Using classical Lyapunov's theorems, with the help of CE (2) it is possible to conduct the investigation of stability in the first approximation. For this purpose we have algorithmically implemented both the Routh-Hurwitz criterion of absence of roots with positive real parts and the criterion, in accordance with which all the roots of the polynomial (2) are negative real. By the complete eq. (1) we investigate the problem of effect of the structure of forces on the stability of the trivial solution of the system (1) [6]. If, using the potentialities of the system of analytical computations, one performs the substitution $\ \Lambda = \mu+\xi\$, where $\ \xi\$ is an assigned shift from the given imaginary axis, then it is possible to solve the problems of estimation of a given degree of stability.

3.2. Construction of first integrals for differential equations of motion. Algorithm to construct first integrals for the systems described is based on invariant relations (differential consequences) of equations of motion [5]. For equations (2.2.3) we have invariant correlation :

$${d\over dt}\left(\ T_{2}- T_{0}- {\cal U} - \sum^{k}_{r=1} \t... ...Q^{*}_{r}\ \left(\ \tilde{a}_{r,k+1} - \Omega_{r}\ \right) \ -$$

$$-\ \sum^{k}_{r,m,s=1}\ {\partial T\over \partial \Omega_{r}}\... ...} {\partial b_{mr} \over \partial t}\ \right) \ . \eqno(3.2.1)$$

In certain conditions from (1) we can obtain quadratic first integral with respect to quasi-velocities. For example,

Statement 1. If lagrangian (2.1.1) and kinematic correlations (2.2.1) do not depend explicitly on time and nonpotential forces $\ Q^{*}_{r}\$ are absent, there exists generalized Jacobi integral:

$$T_{2}(\Omega ,q) - T_{0}(q) - {\cal U}(q) - \sum^{k}_{r=1} \t... ...} \over \partial \Omega_{r}}\ \right)\ = const\ . \eqno(3.2.2)$$

If kinematic correlations are also homogeneous $\ (\ \tilde{a}_{s,k+1} = 0\ ,\ \ s=1,\ldots,k\ )\$, Jacobi integral takes the form of :

$$T_{2}(\Omega,q) - T_{0}(q) - {\cal U}(q) = const\ . \eqno (3.2.3)$$

Assuming that quasi-velocities are equal to generalized velocities, from (1) - (3) we get appropriate statements for systems in Lagrange variables.

Statement 2. If for $\ q_{i}\ \ (\ i=1,\ldots, l\ ; \ \ l < k\ )\$ $\ \partial {\cal U}/ \partial q_{i}= 0\$ and

$$\sum^{k}_{r=1} \left(\ {\partial \tilde{a}_{r,k+1} \over \par... ...a_{r}}\ + {\partial T\over \partial q_{i}}= 0\ , \eqno (3.2.4)$$

then system (2.2.3) has ignorable integral

$$\sum^{k}_{s=1}\ {\partial T\over \partial \Omega_{s}}\ \tilde{a}_{si}\ =const\ .\ \ \ \eqno (3.2.5)$$

Statement 3. If in stationary mechanical system following correlation

$$\sum^{k}_{r=1}\ \left(\ \left(\ {\partial T\over \partial q_{... ...al T\over \partial \Omega_{r}}\ \right) \ = 0\ , \eqno (3.2.6)$$

holds, then equations system allows the integral

$$\sum^{k}_{s=1}\ {\partial T\over \partial \Omega_{s}}\ \tilde{a}_{s,k+1}\ = const\ . \eqno (3.2.7)$$

Statement 4. If body system is described by Lagrange equations in generalized coordinates $\ q\$ and velocities $\ \dot q\$, and lagrangian does not contain coordinates $\ q_{i}\ \ \ (\ i= 1,\ldots,l\ )\$, there exist ignorable integrals

$${\partial L\over \partial \dot q_{i}}\ =\ const\ . \eqno(3.2.8)$$

Let for carrying body of the system generalized coordinates $\ q_{1}\ ,\ q_{2}\ ,\ q_{3}\$ describe rotational motion, $\ q_{4}\ ,\ q_{5}\ ,\ q_{6}\$ - motion of point $\ O_{1}\$, and coordinates $\ q_{j}\ \ (j=7,\ldots, k)\$ describe relative motions of carried bodies. We choose quasi-velocities for carrying body equal to projections of absolute angular velocity $\ {\mbox{\boldmath$\omega$ }}_{1}\$and linear velocity $\ {\bf v}\$ of the point $\ O_{1}\$ to axes associated with body. Let quasi-velocities of carried bodies be equal to projections of relative velocities or to generalized velocities. In these conditions for Euler-Lagrange equations following correlations take place :

$${d\over dt}\left(\left( {\partial T\over \partial{\bf v}}\rig... ...left( {\partial T\over \partial {\bf v}} \right)^{T}{\bf F}\ ;$$

$${d\over dt}\left(\left( {\partial T\over \partial {\bf\Omega}... ...}} = 2 {\bf M}^{T} {\partial T\over \partial {\bf\Omega }}\ ;$$

$${d\over dt}\left(\left({\partial T\over \partial {\bf\Omega} ... ...al T\over \partial {\bf v}}\right)^{T}{\bf M}\ ,\ \eqno(3.2.9)$$

where $\ {\bf M} = col\ (P_{1},\ P_{2},\ P_{3})\ ,\ \ {\bf F}= col\ (P_{4},\ P_{5},\ P_{6})\$,

$$\breve{ v}=\left( \begin{array}{ccc} 0 & -v_{3} & v_{2} \\ v... ...rtial \Omega_{2}},\ {\partial T\over \partial \Omega_{3}})\ .$$

Let us choose two constant directions $\ \zeta \$ and $\ {\tilde \zeta}\ ,$with matrices of direction cosine in coordinates of carrying body $\ \Gamma \$ and $\ \tilde {\Gamma}\$ of size $\ 3\times 1\ .$If we complement Euler-Lagrange equations system with Poisson equations

$$\dot \Gamma = - \breve{\Omega}\Gamma \ \ \ {\rm and} \ \ \ \d... ...ega_{1} \\ -\Omega_{2} & \Omega_{1} & 0 \end{array}\right)\ ,$$

we easyly can get following correlation :

$${d\over dt}\left(\ \Gamma^{T} {\partial T\over \partial {\bf\... ...\Gamma^{T}{\bf M} + \tilde{\Gamma}^{T}{\bf F}\ . \eqno(3.2.10)$$

From correlations (9), (10) in appropriate conditions we can get first integrals of differential equations of motion.

Statement 5. If $\ ( \partial T/\partial {\bf v} )^{T}\ {\bf F} = 0\$, then

$$\left( {\partial T\over \partial{\bf v}}\right)^{T} {\partial T\over \partial {\bf v}}= const.$$

Let point $\ O_{1}\$ be fixed.

Statement 6. If $\ {\bf M}^{T}\ (\partial T/ \partial {\bf\Omega})\ = 0\$, then

$$\left( {\partial T\over \partial {\bf\Omega }}\right)^{T} {\partial T\over \partial {\bf\Omega }} =const\ .$$

In these conditions also takes place generalized integral of areas:

$$(\ {\bf e\ \cdot }\ {\partial T\over \partial \bf\Omega }\ )\ =\ const\ ,$$

where $\ {\bf e}\$ is unit vector (movable), directed parallel to the vector of kinetic momentum of the body system.

Statement 7. If $\ {\bf M} = c^{T}\ \breve{\Omega}\ ,$where $\ c\ -$ is a constant matrix $\ 3\times1\$, then

$$(\ c + {\partial T\over \partial {\bf\Omega }}\ )^{T} (\ c + {\partial T\over \partial {\bf\Omega }}\ ) = const\ .$$

Statement 8. If $\ \Gamma^{T}\ {\bf M} = 0\$ , then $\ \Gamma^{T}\ \left(\partial T/ \partial {\bf\Omega }\right) = const\ .$

Statement 9. If $\ {\bf M} = c^{T}\ \breve{\Omega}\$, then

$$\Gamma^{T}\ \left(\ {\partial T\over \partial {\bf\Omega }} + c\ \right) = const\ .$$

Using above differential consequences and statements it's possible to automatically construct quadratic and linear first integrals with respect to quasi-velocities, if conditions of their existence are met. Verification of these conditions and explicit form of integrals is possible only with the help of CAS. The CAS should be able to do complex trigonometrical transformations.

For equations of the form (3.1.1) algorithms of obtaining differential consequences and invariant correlations are described in [6, 8]. For linear systems (3.1.1) obtaining quadratic first integrals can be reduced to systems of linear algebraic equations.

3.3. Investigation of stability by Lyapunov's second method. Invariant relations, first integrals obtained from it and other invariants may be used for finding steady motions or invariant manifolds of steady motions as well as for constructing Lyapunov functions in investigation of stability of the latter [14]. To this end, an algorithm of Routh-Lyapunov's theorem is used. In accordance with the latter theorem we come to the problem of conditional extremum.

Let for the system of ordinary differential equations

$$\dot x= X(x)\ ,\ \ \ \ x \in R^{k}\ . \eqno(3.3.1)$$

some number of first integrals

$$V_{0}(x) = c_{0}\ ,\ \ V_{1}(x) = c_{1}\ , \ \ldots\ ,\ \ V_{m}(x) = c_{m}\ .$$

be known. Having chosen these integrals as the basis, we define the linear space of system's integrals over $\ R\$ :

$$K_{J} = \sum^{n}_{j=1}\lambda _{i_{j}} V_{i_{j}}(x)\ ,\ \ \ J = \{i_{0},\ldots ,i_{n}\} \in (0,1,\ldots ,m)\ .$$

Let us find steady motions (invariant manifolds of steady motions - IMSM) as solutions of the system of equations

$${\partial K\over \partial x_{i}} = f_{i}(x_{1},\ldots x_{k}, ... ...i_{1}},...,\lambda _{i_{n}}) = 0\ ,\ \ \ \ (i=1,\ldots \ ,\ k)$$

$$V_{i_{1}}(x) = c_{i_{1}}\ ,\ \ldots \ ,\ V_{i_{n}}(x) = c_{i_{n}} \eqno (3.3.2))$$

In the regular case, for a system of interconnected bodies the steady motions are found as solutions of equations of stationarity of the bundle $\ K_{J}\$ of first integrals of motion differential equations. Sufficient conditions of stability for these motions may be found as the conditions of sign-definiteness of the second variation $\ \delta^{2}K_{J}\$ for constant values of the rest of the first integrals. This problem has a good algorithm of solving in the case when side by side with the energy integrals there are linear first integrals with respect to velocities or quasi-velocities (e.g. ignorable ones) [2, 4]. In the case, when there are (i) first integrals of another structure or (ii) singularity of the system (2) (this is characteristic of electromechanical systems), difficulties of the investigation increase. When there are no first integrals in the system, for the purpose of constructing the Lyapunov functions in the process of investigation of stability by Lyapunov's method it is possible to use differential corollaries of equations of disturbed motion [6].

3.4. Reduction to the normal form. The preliminary reduction of the system to one of the normal forms is an efficient aid in solving problems of qualitative investigation of ordinary differential equations. Poincare's form [15, 16] is one such normal forms for differential equations of the general form (3.3.1).

Normalization of systems of analytical ordinary differential equations

$$\dot x\ =\ A x\ +\ f(x)\ \ ,$$

where $\ A\ -$ constant $\ n \times n \ -$ matrix, $\ \ f(x)\ -$ nonlinear terms of the 2d and higher degrees of smallness, undergoes decomposition into the 2 parts : (i) the reduction of the matrix $\ A\$ to the Jordan canonical form and (ii) nonlinear nonsingular substitution of variables, which removes most of nonlinear terms in the right-hand side of the equation. The linear normalization determines the eigenvalues $\ \Lambda_1,\ldots,\Lambda_n\$ of the matrix $\ A\ ,\$ and, as soon as eigenvectors and adjunct vectors (if any) are found, allows to obtain the Jordan canonical matrix and determines the corresponding linear transformation of variables (as well as its converse).

As a result, the system of ordinary differential equations takes the form

$$\dot y\ =\ \bar {A} y\ +\ \varphi(y)\ . \eqno {(3.4.1)}$$

Furthermore, in Hamilton-type systems the linear transformation is assumed canonical.

In practical investigations the precision of computing of the eigenvalues was chosen $\ \epsilon \ = \ 10^{-11}\ .$Zero and multiple roots were determined with the same precision. Eigenvectors and adjunct vectors were found with the precision of $\ 10^{-8}\ .$

Linear normalization has been conducted only in the numerical form.

Nonlinear normalization has been conducted by Poincare's method [15,16], which retains only the resonance terms of the form $\ ({\bf s},\ {\bf\Lambda})\ -\ \Lambda_i\ =\ 0\ ,\$ where $\ {\bf s}\ -\$ integer numeric vector of nonnegative indices of powers of the variables. Nonlinear transformation is realized in the form of the formal series:

$$y\ =\ z\ +\ \psi(z)\ . \eqno {(3.4.2)}$$

The normalized system of ordinary differential equations for the variables $\ z$

$$\dot z\ =\ \bar A z\ +\ g(z)\ ,$$

contains in $\ g(z)\$ only resonance terms $\ \sum_ {\vert s \vert \geq 2}^{\infty} {b_s\ z^s}\ .$The coefficients $\ b_s\$ are sought simultaneously with the coefficients of the series (3.4.2) by recurrent formulas [16] obtained from the identity with respect to $\ z$:

$$g_{is}(z)\ +\ [(s,\Lambda)- \Lambda_i]\ \psi_{is}(z)\ = \ \varphi(z+\psi(z)) \ +\ \sigma_i \psi_{(i+1) s}\ -$$

$$-\ \sum_{k=1}^{n} {\sigma_{k+1}\ {{\partial \psi_i(z)}\over {... ...+1} {{\partial \psi_{iR}(z)}\over {\partial z_k}} g_{kQ} }\ \ ,$$

where $\ \sigma_i\ -\$ nondiagonal elements of Jordan cells.

Nonlinear normalization is executed in both numerical and symbolic forms.

4. Examples.

With the help of packages in a symbolic form the authors investigated rather large number of concrete mechanical systems. In particular:

- nonlinear equations of motion of a series of robotic systems are obtained (for example, robot - manipulator of 6 links, walking platform on four legs with two degrees of freedom each) ;

- problems of dynamics and stability of steady motions in space systems were considered (for example, satellite with a gravitational system of the stabilization on circular orbit are considered at various control laws by a gravitational stabilizer, gyroscopic frame in Newtonean central field of forces );

- linear equations of motion for a mechanical system with 32 degrees of freedom consisting of 20 absolutely of rigid bodies are constructed.

Application of the software system "Normalization" has allowed us:

- to investigate stability of a nonlinear autonomous mechanical system with two pairs of purely imaginary multiple roots (with two groups of solutions);

- to solve the problem of stabilization of one mechanical Hamilton-type system with the resonance of 1:3 and with nonlinear control.

Let us illustrate our reasoning with a few applied examples of usage of the above software.

Example 4.1. Let us consider gyroscopic frame, placed in Newtonean central field of forces. The mechanical system consists of a carrying frame (a body with a fixed point) and two identical connected two-degree gyroscopes, symmetrically installed in a frame (Fig. 4). Masses and the moments of inertia of housing of gyroscopes are neglected.

Figure 4.

Before exposition of input information about the system we note that the matrix of rotation $\ \alpha^ {ij} \$ and relative angular velocities of bodies $\mbox{\boldmath$\omega$ }^{j}_{i}\$ are not required to be introduce. They are calculated automatically on specific "sequence of rotations". Hereby we specify :

a) Number of axis of rotation (one of numbers 0, 1, 2, 3), i.e. 1 - the rotation is carried out around the axis $\ Ox_ {i}\ ;\$2 - around the axis $\ Oy_ {i} \ ;\$ 3 - around the axis $\ Oz_{i}\ ;\$0 - there is no rotation;

b) angle of rotation.

Input data :

Number of bodies in the system : $\ N = 5\$.

Body 1 is a frame. $\ M_1\$ is mass ;

$${\bf r^{0}_{o_1}} = \left( \begin{array}{c} 0 \\ 0 \\ R \end{... ...ay}{c} \Omega_1 \\ \Omega_2 \\ \Omega_3 \end{array} \right)\ ;$$

$${\rm the\ inertia\ tensor\ :}\hskip 1.5cm \Theta^{o_1}_{*} = ... ... & 0 & 0 \\ 0 & B_{1} & 0 \\ 0 & 0 & C_{1} \end{array} \right)$$

"sequence of rotations" : $\ \ 3\ ,\ \psi\ ;\ 1\ ,\ \theta\ ;\ 3\ ,\ \varphi\$.

Body 2 is the housing of the first gyroscope, it is connected to the frame. $\ M_2 = 0\$;

$${\bf r^{1}_{o_2}} = \left( \begin{array}{c} 0 \\ 0 \\ l \end{... ...}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right)$$

"sequence of rotations" : $\ \ 1\ ,\ \delta\$.

Body 3 is the rotor of the first gyroscope, it is connected to the second body. $\ m\$ is mass of a rotor;

$${\bf r^{2}_{o_3}} = {\bf r^{3}_{c_3}} = {\bf v^{2}_{o_3}} = \... ... A_3 & 0 & 0 \\ 0 & A_3 & 0 \\ 0 & 0 & C_3 \end{array} \right)$$

"sequence of rotations" : $\ \ 3\ ,\ \beta_1\$.

Body 4 is the housing of the second gyroscope, it is connected to a frame. $\ M_4 = 0\$;

$${\bf r^{1}_{o_4}} = \left( \begin{array}{c} 0 \\ 0 \\ -l \end... ...}{ccc} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right)$$

"sequence of rotations" : $\ \ 1\ ,\ -\delta\$.

Body 5 is the rotor of the second gyroscope, it is connected to the fourth body. $\ m\$ is mass of a rotor;

$${\bf r^{4}_{o_5}} = {\bf r^{5}_{c_5}} = {\bf v^{4}_{o_5}} = \... ... A_3 & 0 & 0 \\ 0 & A_3 & 0 \\ 0 & 0 & C_3 \end{array} \right)$$

"sequence of rotations" : $\ \ 3\ ,\ \beta_2\$.

Program output :

Matrices of rotation :

$$\alpha^{01} = \left( \begin{array}{ccc} (\cos\varphi\cos\psi-... ...\psi & -\cos\psi\sin\theta & \cos\theta \end{array} \right)\ ;$$

$$\alpha^{12} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos... ...\beta_1 & \cos\beta_1 & 0 \\ 0 & 0 & 1 \end{array} \right)\ ;$$

$$\alpha^{14} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos... ...\beta_2 & \cos\beta_2 & 0 \\ 0 & 0 & 1 \end{array} \right)\ .$$

Relative angular velocities :

$$\mbox{\boldmath $\omega$}^{1}_{2} = \left( \begin{array}{c} \... ...\begin{array}{c} 0 \\ 0 \\ \dot \beta_2 \end{array} \right)\ .$$

The kinetic energy of the system :

$$T = {1\over 2}[ A \Omega^{2}_{1}+ B \Omega^{2}_{2}+ C \Omega^... ...+ \Omega_{3}\cos \delta (\dot{\beta }_{1} + \dot{\beta }_{2})]$$

£¤¥ $\ \ A = A_{1}+ 2ml^{2}+ 2A_{3}\ ;\ \ \ \ C = C_{1}+ 2A_{3}\sin^{2}\delta + 2C_{3}\cos^{2}\delta\ ;$

$\hskip 1.2cm B = B_{1}+ 2ml^{2}+ 2A_{3}\cos ^{2}\delta + 2C_{3}\sin ^{2}\delta$

$\Omega_ {1},\ \Omega_ {2},\ \Omega_{3} \$ are the projections of absolute angular velocity of the frame on the axes $\ x_ {1},\ y_ {1},\ z_ {1} \$. They are selected as quasi-velocities.

$$\Omega_1 = \dot \psi\sin\theta\sin\varphi+\dot \theta\cos\var... ...rphi\ ; \ \ \ \Omega_3 = \dot \psi\cos\theta + \dot \varphi\ ,$$

The remaining quasi-velocities coincide with generalized velocities $\ \dot{\delta }\ ,\ \dot{\beta }_{1}\ ,\ \dot{\beta }_{2}\$.

Approximate force function $\ { \cal U}$ (1.1.13) in Newtonean field of gravitation to the fixed center $\ O$

$${\cal U} = {3\nu \over R^{3}} \ [\ {1 \over 6} (A_{1}+B_{1}+C... ...C_{1}\gamma ^{2}_{3}) + {1 \over 3} (2A_{3} + C_{3}- ml^{2}) -$$

$$- A_{3}(\gamma^{2}_{1} + \gamma^{2}_{2} \cos^{2}\delta + \gam... ...\gamma ^{2}_{3}\cos ^{2}\delta ) + ml^{2} \gamma ^{2}_{3}\ ] -$$

$$- M_{1}{\nu\over R^{2}} (x_{c}\gamma _{1} + y_{c}\gamma _{2}+ z_{c}\gamma _{3}) + (M_{1}+ 2m) {\nu \over R} + const\ .$$

Here $\ \gamma_{1} = \sin \varphi \sin \theta \ , \ \ \gamma_{2} = \sin\theta \cos \varphi\ , \ \ \gamma_{3} = \cos \theta \$ are direction cosines of axis $\ Oz\$ with axes $\ x_ {1},\ y_ {1},\ z_ {1} \$;

The force function of elastic forces of the spring device :

$${\cal U}(\delta )\ =\ -2\ \chi\ \sin^{2} {\delta \over 2}\ \ , \hskip 1.2cm {\rm where } \hskip 0.8cm \chi\ = \ const\ .$$

Further four first integral of motion equations (integral of an energy and three ignorable on coordinates $\ \psi\ ,\ \beta_1\ ,\ \beta_2\$) are found, equations of steady motions and sufficient conditions of stability of steady motions are obtained by computer [17]. All the expressions obtained by computer have the symbolic form.

Example 4.2. The modelling of the electrical circuit consisting of two parts connected inductively (Fig. 5).

Figure 5.

1, 2, 3, 4 are the numbers of the circuit's nodes. The values in the interior parentheses are the branches $\{i,j\}$ and the elements of these ones.

Input data.

The number of parts of the electrical circuit connected inductively: 2

The first part of the electrical circuit:

$\hskip 0.2cm \{\{\{1, 2\},\ \{L_1, C_1, R_1, e_1\}\},\ \{\{1, 2\},\ \{C_5\}\},$ $\{\{1,\ 2\},\ \{C_2, L_2, R_2, M\}\}\}$

The second part of the electrical circuit:

$\hskip 0.2cm \{\{\{3, 4\},\ \{R_3, L_3, C_3, M\}\},$ $\{\{3, 4\},\ \{e_2\}\},\ \{\{3, 4\},\ \{C_4, L_4, R_4\}\}\}$

The order of path-tracing the graph:

"+" - to trace the graph's vertices in the direction of growth of their numbers ("-" - inverse order).

Program output.

The Lagrange's function:

$$\frac{1}{2} L_1 (\dot q_1 - \dot q_2)^2 + \frac{1}{2} L_2 \d... ...{2} L_3 (\dot q_3 - \dot q_4)^2 + \frac{1}{2} L_4 \dot q_4^2 -$$

$$-\frac{q_1^2}{2 C_5} + e_1 (q_1 - q_2) - \frac{(q_1 - q_2)^2}... ...- e_2 q_3 - \frac{(q_3 - q_4)^2}{2 C_3} - \frac{q_4^2}{2 C_4}$$

The Rayleigh's function:

$$\frac{1}{2} R_1 (\dot q_1 - \dot q_2)^2 + \frac{1}{2} R_2 \do... ...1}{2} R_3 (\dot q_3 - \dot q_4)^2 + \frac{1}{2} R_4 \dot q_4^2$$

The differential equations of the electrical circuit in Lagrange's form of 2-d kind:

$$L_1 \ddot q_1 - L_1 \ddot q_2 + R_1 \dot q_1 - R_1 \dot q_2 + (\frac{1}{C_1} + \frac{1}{C_5}) q_1 - \frac{q_2}{C_1} - e_1 = 0$$

$$-L_1 \ddot q_1 + (L_1 + L_2) \ddot q_2 + M \ddot q_3 - M \ddot q_4 - R_1 \dot q_1 + (R_1 + R_2) \dot q_2 - \frac{q_1}{C_1} +$$

$$+ (\frac{1}{C_1} + \frac{1}{C_2}) q_2 + e_1 = 0$$

$$L_3 \ddot q_3 - L_3 \ddot q_4 + M \ddot q_2 + R_3 \dot q_3 - R_3 \dot q_4 + \frac{q_3}{C_3} - \frac{q_4}{C_3} + e_2 = 0$$

$$- L_3 \ddot q_3 + (L_3 + L_4) \ddot q_4 - M \ddot q_2 - R_3 ... ...4 - \frac{q_3}{C_3} + (\frac{1}{C_3} + \frac{1}{C_4}) q_4 = 0$$

Example 4.3. In the process of investigation of the phase space of a mechanical or electric system there often appears the need to solve algebraic equations (or those reducible to them). In the contemporary universal computer algebra systems there exists a tool for solving such problems, which allows to simplify such systems. Gröbner's bases or their analogies may be considered as such tools.

In applied problems, direct usage of Gröbner's bases does not always entail in satisfactory result. For example, for the systems containing parameters the investigation requires greater efforts and does not guarantee a success even when a computer is used. In such situations the following technique of analysis of algebraic systems with parameters may be useful.

Let a family of solutions be known for an algebraic system. The existence of such a family in concrete cases may be grounded by the mechanical or physical character of the problem. We take an interest in such solutions of the system which (i) contain some elements of this family or (ii) adjunct to the family. Such problems arise in different directions of bifurcation theory. Now, let there be given a system of algebraic equations with parameters

$$f_{i}\ (\ x_{1},\ldots,x_{n},\ \lambda_{1},\ldots,\lambda_{k}\ )\ =\ 0 \ \ \ \ \ \ \ (i=1,\ldots,n) \eqno(4.3.1)$$

for which there is a known family of solutions $\ x_{1}^0,\ldots,x_{n-1}^0\$with parameters $\ x_{n}\$ for some relationships between $\ \lambda_{1},\ldots,\lambda_{k}$

$$f_{i}\ (\ x_{1}^0,\ldots,x_{n-1}^0,\ x_{n}, \ \lambda_{1},\ld... ...lambda_{k}\ )\ =\ 0 \ \ \ \ \ \ (i=1,\ldots,n)\ . \eqno(4.3.2)$$

Assume that the rank of the Jacoby matrix for the system (2) with respect to $\ \lambda_{1},\ldots,\lambda_{k}$ is $\ l< k\$ and consider the system of equations

$$f_{i}\ (\ x_{1}^0,\ldots,x_{n-1}^0,\ x_{n},\ \lambda_{1},\ldots,\lambda_{k}\ )\ =\ 0 \ \ \ \ \ \ \ (i=1,\ldots,n)$$

$$\Delta = \tilde \varphi\ (\ x_{1}^0,\ldots,x_{n-1}^0,\ x_{n},\lambda_{1},\ldots,\lambda_{k}\ )\ =\ 0\ , \ \ \ \eqno(4.3.3)$$

where $\ \Delta\ -$ Jacobian with respect to $\ x_{1},\ldots,x_{n}\$for the system (1), which has been computed for the solutions of our family. It is assumed that it is not identically zero.

Analysis of the system (3) allows us to find, generally speaking, those relationships between $\ \lambda_{1},\ldots,\lambda_{k}$, for which the Jacoby matrix's rank $\parallel \partial f_{i}/\partial x_{j} \parallel$on the chosen family $\ x_{1}^0,\ldots,x_{n-1}^0,\ x_{n}\$becomes smaller.

These bifurcation values of the parameters may be substituted into the initial system (1), what allows to simplify it and find new solutions or their families.

The proposed procedure may successfully be employed, for example, in analysis of equations of stationarity of mechanical systems' first integrals.

Consider a concrete example of application of above approach in analysis of bifurcations in the neighbourhood of the family of permanent rotations of S.V.Kovalevskaya's top about the vertical axis $\ x\$ [14].

To this end, let us compose a complete linear bundle of the problem's first integrals:

$$K = H - \lambda_{1}V_{1}-\lambda_{2}{V_{2} \over 2} - \lambda_{3}{V_{3} \over 3},$$

where $\ \lambda_{1},\ \lambda_{2},\ \lambda_{3}\ -$ constants,

$$H = 2p^2 + 2q^2 + r^2 +2x_{0}\gamma_{1}=2h\ , \ \ \ \ \ V_{1} = 2p\gamma_{1} + 2q\gamma_{2} + r\gamma_{3}=m\ ,$$

$$V_{3} = \gamma_{1}^2 + \gamma_{2}^2 + \gamma_{3}^2 = 1\ , \ \... ...\gamma_{1})^2 + (2pq - x_{0}\gamma_{2})^2=k^2\ . \eqno(4.3.4)$$

Now write the steady-state conditions $\ K\$ with respect to the phase variables:

$${\partial K \over \partial p} = 2[(\lambda_{2}x_{0}p-\lambda_... ... + \lambda_{2}x_{0}q\gamma_{2} - \lambda_{2}p(p^2 + q^2)]=0\ ;$$

$${\partial K \over \partial q} = 2[(\lambda_{2}x_{0}p-\lambda_... ... - \lambda_{2}x_{0}q\gamma_{1} - \lambda_{2}q(p^2 + q^2)]=0\ ;$$

$${\partial K \over \partial \gamma_{1}} = -(\lambda_{3}+x_{0}^... ...2\lambda_{1}p + \lambda_{2}x_{0}(p^2 - q^2)]=0\ ; \eqno(4.3.5)$$

$${\partial K \over \partial \gamma_{2}} = -(\lambda_{3}+x_{0}^2\lambda_{2})\gamma_{2}-2q(\lambda_{1}-\lambda_{2}x_{0}p)=0\ ;$$

$${\partial K \over \partial r} = r - \lambda_{1}\gamma_{3}=0\ ... ...rtial \gamma_{3}} = -\lambda_{1}r - \lambda_{3}\gamma_{3}=0\ .$$

It can be easily determined that the system of algebraic eqs. (4), (5) has the elements of the family of the body's permanent rotations about its vertical axis $\ x\$ in the capacity of its solutions:

$$\gamma_{1}=1\ ;\ \ \ p^0=\lambda_{1}^0=m\ ; \ \ \ q=r=\gamma... ...gamma_{3}=0\ ; \ \ \ \lambda_{3}^0=x_{0}-2m^2\ . \eqno(4.3.6)$$

Having substituted the solutions (6) into eqs. (5), we obtain the condition for the parameters

$$p - \lambda_{1} + \lambda_{2}p(x_{0} - p^2)=0\ ; \ \ \ x_{0}... ...\lambda_{3} - \lambda_{2}x_{0}(x_{0} - p^2)=0\ . \eqno(4.3.7)$$

The Jacobian for the system (5) with respect to $\ p,\ q,\ r,\ \gamma_{1},\ \gamma_{2},\ \gamma_{3}\$ for the solution (6) writes:

$$\Delta_{0}= -4(\lambda_{3}+\lambda_{1}^2 )\ \{\ [\ 2\lambda_{... ...lambda_{1}\lambda_{2}x_{0}p-2\lambda_{2}\lambda_{3}p^2\ ]^2\ -$$

$$-(\lambda_{3}+x_{0}^2 \lambda_{2})^2 \lambda_{2}^2 (p^2-x_{0})^2\ \} \ . \eqno(4.3.8)$$

Hence, for obtaining bifurcation values of the parameters $\ \lambda_{1},\ \lambda_{2},\ \lambda_{3}\$ it is necessary to find solutions of the algebraic equations (7), which are complemented with the condition $\ \Delta_{0} = 0\$. Furthermore, as obvious from (8), we managed to factorize the expression of the determinant.

Complete analysis of the system (7), (8) is rather cumbersome. So, we restrict our investigation, while reducing it to only the first obvious subcase. Consider the following system of algebraic equations with unknown variables $\ \lambda_{1},\ \lambda_{2},\ \lambda_{3},\ p\$ :

$$p-\lambda_{1}+\lambda_{2}p(x_{0}-p^2)=0\ ; \ \ \ x_{0}-2\lam... ...{2}x_{0}(x_{0}-p^2)=0\ ; \ \ \ \lambda_{3}+\lambda_{1}^2=0\ .$$

The standard procedure puts the following collection of equations in correspondence to the system in the capacity of its Gröbner basis for the corresponding ordering of unknown variables (computations have been conducted with the CAS on IBM PC/486):

$$(\lambda_{2}x_{0}+1)p + p^3\lambda_{2}+\lambda_{1}\ =\ 0\ ;$$

$$(8p^2-\lambda_{2}x_{0}^2)\lambda_{1}^2 - (3\lambda_{2}x_{0}+... ...+ 3\lambda_{2}^2x_{0}^4 + 3\lambda_{2}x_{0}^3 - 6x_{0}^2=0\ ;$$

$$(16p^4-20p^2x_{0}-5\lambda_{2}x_{0}^3+8x_{0}^2)\lambda_{1}^2+... ...2}x_{0}-4)p^2x_{0}^2+(\lambda_{2}x_{0}+5)\lambda_{1}^4x_{0}\ -$$

$$-\ \lambda_{1}^6-8\lambda_{2}x_{0}^4+8x_{0}^3=0\ ; \eqno(4.3.9)$$

$$(\lambda_{2}x_{0}-4)p^2-\lambda_{2}x_{0}^2+x_{0})\lambda_{1}^... ...a_{2}x_{0}^2+\lambda_{1}^4-2\lambda_{2}x_{0}^3+2x_{0}^2\ =0\ ;$$

$$(\lambda_{2}x_{0}+2)p^2-2p^4\lambda_{2}-\lambda_{1}^{2}+\lambda_{2}x_{0}^{2}- x_{0}\ =\ 0\ ;$$

$$(4p^4+x_{0}^{2})\lambda_{1}^{2}-(\lambda_{1}^{4}+4x_{0}^{2})p^{2}\ =\ 0\ .$$

Analysis of these equations allows one to find the bifurcation values of parameters, which correspond to the branches from the family of permanent rotations (6) of invariant manifolds of the problem's steady-state motions.

With the aid of the CAS it is possible to find all the solutions of the above nonlinear algebraic system. Without suggesting the result in complete form, let us outline the plan of further actions. From the last equation of (9) we have the following solutions for $\ \lambda_{1}^2 (p)\$:

$$\lambda_{1}^2=4p^2\ ; \ \ \ \lambda_{1}^2 = x_{0}^2/p^2\ .$$

The penultimate equation of (9) for $\ \lambda_{1}^2=\ x_{0}^2/p^2\$gives the following value $\ \lambda_{2}(p)\$:

$$\lambda_{2}=1/p^2\ \ \ {\rm for}\ \ \ x_{0}p^{2}\ -\ 2p^{4}\ +\ x_{0}^{2}\ \neq\ 0\$$

Having removed $\ p\$, we obtain

$$x_{0}^{2}\lambda_{2}-\lambda_{1}^{2}\ =\ 0\ ; \ \ \ \lambda_{3}=-\lambda_{1}^{2} \ \ , \eqno (4.3.10)$$

which is one-parameter family of bifurcation values of the parameters $\ \lambda_{1}, \ \lambda_{2},\ \lambda_{3}\ .$Substitution of these relation in the rest of equations of (9) converts them into identities. Eqs. (5) for $\ \lambda_{1}\ ,\ \lambda_{2}\ ,\ \lambda_{3}\ ,$which satisfy the relations (10), define the family of IMSMs

$$q=0\ ;\ \ \ \ p=x_{0}/\lambda_{1}\ ;\ \ \ \ r=\lambda_{1} \gamma_{3}$$

with the vector field

$$\dot{r}\ =\ -x_{0}\ \gamma_{2}\ ; \ \ \ \dot{\gamma}_{1}\ =\... ...dot{\gamma}_{2}\ =\ r\ (\ (x_{0}/\lambda_{1})-\gamma_{1}\ )\ ,$$

on it, which branches from the permanent rotations (6) for

$$p\neq 0 \ ; \ \ \ p^2\neq \ x_{0}\ ; \ \ \ p^2\neq \ - x_{0}/2\ .$$

This result serves as the complement to the analysis of bifurcations in the neighbourhood of the family (6) given in [7].

REFERENCES

1. Bourlakova L.A., Irtegov V.D., Pochtarenko M.V. Use of the symbolic calculations on the computer in some problems of a mechanics. // Proceeding International Conference on Systems and Techniques of Analytical Computing and Their Application in Theoretical Physics. Dubna, JINR, $\ D-11-80-13,\ 1980,\ 137-142.$ (in Russian)

2. Banshchikov A.V., Bourlakova L.A., Irtegov V.D. Qualitative investigations of Systems with the Aid of Computer Algebra. // IV Intern. Conf. Computer Algebra in Physical Research. World Scientific Publishing Co. Pte. Ltd, Singapore, $\ 1991,\ 394-398.$

3. Banshchikov A.V., Bourlakova L.A., Irtegov V.D. The use of computer symbolic calculations in problems of motion stability. // Computer Algebra and Its Applications to Mechanics. Nova Science Publishers, Inc., USA, $\ 1992, \ 153-158.$

4. Banshchikov A.V., Bourlakova L.A., Ivanova G.N., Irtegov V.D., Usova Y.O. Package of symbolic computations for investigation of dynamics of system bodies. // Application packages. Software of mathematical modeling. Algorithms and algorithmic languages. Moscow: Nauka, $\ 1992,\ 63-71\ .$(in Russian)

5. Banshchikov A.V., Bourlakova L.A., Irtegov V.D., Novikov M.A. Algorithms of qualitative investigations of complex systems. // Cybernetics and system analysis. (Kibernetika i sistemnii analiz.) $\ 1992,\ N^ {0} \ 1,\ 129-139.$ (in Russian).

6. Banshchikov A.V., Bourlakova L.A. Information and Research System "Stability". //Informations RAS. The theory and control systems. (Izvestia RAN. Teoria i sistemi upravlenia.) $\ 1996.\ N^ {0} \ 2.\ 13 - 20.$ (in Russian).

7. Irtegov V.D., Titorenko T.N . On the modelling and the research of some problems with the use of computer algebra. // Programming. (Programirovanie.) $\ 1997.\ N^ {0} \ 1.\ 68-74.$ (in Russian).

8. Banshchikov A.V., Bourlakova L.A. On algorithms of symbolic computation at the research of the stability. // Programming. (Programirovanie.) $\ 1997.\ N^ {0} \ 3,\ 72-80.$ (in Russian).

9. Lurie A.I. Analytical Mechanics.- Moscow: Phizmathgiz, 1961.

10. Hantmaher F.P. The lectures on analytical mechanics. Moscow: Nauka, 1966. (in Russian)

11. Brayton R.K., Moser J.K. A theory of nonlinear networks - 1 //Quaterly of Applied Mathematics, 1964, vol.22, 1, p.1-33.

12. Brayton R.K., Moser J.K. A theory of nonlinear networks - 2 //Quaterly of Applied Mathematics, 1964, vol.22, 2, p.81-104.

13. Titorenko T.N. The Hamilton's equations in the case of the degenerated dynamic systems. //Thesises of the reports of the 4-th International conference of women-mathematicians. Volgograd, 1996, p.122-123. (in Russian)

14. Irtegov V.D. Invariant manifolds of stationary motions and their stability. Novosibirsk. Nauka. 1985. ( in Russian )

15. Poincare H. The new methods of celestial mechanics. Moscow: Nauka, vol.1, 1971. (in Russian) ( Les méthodes nouvelles de la Mécanique céleste, t. 1. Paris, Gauthier-Villars, 1892.)

16. Bruno A.D. The analytical form of differential equations. //Proceeding Moscow mathematical society, 1971, vol.25, p.119-262. (in Russian)

17. Banshchikov A.V. Research of a stability of gyroframe by tools of computer algebra. // Mathematical modelling. (Matematicheskoe modelirovanie.) $\ 1996.\ N^{0}\ 4,\ 67-78.$ (in Russian).