=Paper=
{{Paper
|id=Vol-2638/paper3
|storemode=property
|title=Application of software tools for symbolic description and modeling of mechanical systems
|pdfUrl=https://ceur-ws.org/Vol-2638/paper3.pdf
|volume=Vol-2638
|authors=Andrei Banshchikov,Alexander Vetrov
|dblpUrl=https://dblp.org/rec/conf/iccs-de/BanshchikovV20
}}
==Application of software tools for symbolic description and modeling of mechanical systems==
https://ceur-ws.org/Vol-2638/paper3.pdf
Application of software tools for symbolic description
and modeling of mechanical systems
Andrei Banshchikov, Alexander Vetrov
Matrosov Institute for System Dynamics and Control Theory of Siberian Branch of Russian
Academy of Sciences, Irkutsk, Russia
E-mail: bav@icc.ru
Abstract. The paper presents two software tools (graphical editor and software package).
The editor is designed for the formation of a symbolic description of a mechanical system
using the Lagrange formalism. A system of the absolutely rigid bodies connected by joints
is considered as a mechanical system. The editor is a user interface by which one sets the
structure of the interconnection of bodies (system configuration) as well as the geometric and
kinematic characteristics for each body of the system. The created structure and the entered
data are automatically presented in the form of a text file, which is used as an input file for the
software package for modeling mechanical systems in a symbolic form with a computer. The
use of these software tools is shown in detail in the example of the analysis of the dynamics of
a satellite with a gravitational stabilizer in a circular orbit. For this system, the kinetic energy
and force function of an approximate Newtonian gravitational field were obtained, nonlinear
and linearized equations of motion were constructed, and the question of the stability of the
relative equilibrium position was considered.
1. Introduction
In recent years, the importance of computer algebra systems (CAS) [1, 2] and the methods
of their use [3] for scientific computations in various fields of natural sciences and engineering
knowledge has grown significantly. Many years of experience of the authors [4] show that CASs
turned out to be quite effective at the stages of constructing mathematical models and analyzing
their properties. The combination of methods of analytical mechanics and the theory of stability
of motion with the algorithms of CAS allows one to create quite powerful tools for the study of
mechanical systems.
The characteristic feature of our investigations is that practically all the computations on a
computer are conducted in symbolic (analytical) form. When it is difficult to conduct research
exclusively in symbolic form, one can proceed to “digitizing” some (or all) of the parameters of
the problem under consideration and, thereby, close the question of rounding errors. This allows
the researcher to penetrate into the qualitative nature of the processes as deeply as possible and
move on to a quantitative (numerical) assessment at a later stage.
As can be seen from the publications, more often there is an approach to using the CAS
as a calculator for solving a particular problem. There is another approach when, based on
Copyright c 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution
4.0 International (CC BY 4.0).
Figure 1. The structure of interconnection between the bodies.
the internal programming language of CAS (in our case – “Mathematica” [5]), a software is
developed to solve a specific class of problems.
Some algorithms for analyzing the dynamics and stability of motion of mechanical systems
using symbolic computations are presented in [6, 7]. The application of computer algebra
methods to the problems of celestial mechanics has rich history and till today attracts academic
attention (see, for example, [8, 9]).
2. Graphic editor for the description of mechanical systems
As a mechanical system, a system of the absolutely rigid bodies Bi (i = 1, N ), connected by
one-two-three-degree joints is considered. The body B1 is a carrier. Its position is determined
relative to the inertial coordinate system (CS) Σ0 . Let us associate the coordinate system Σ1
to the rigid body B1 . The starting point of that system will be at the point O1 . The bodies
are connected so that for each Bk there is a common point Ok ( Ok ∈ Bk , Ok ∈ Bj , j = 1, N ,
k = 2, N , j < k ) (see Figure 1), or joints allow translational displacements Bk relative to Bj .
With the body Bk in a point Ok ∈ Bk we will associate a CS Σk . The angular position of Σk
relative to Σj is determined by the matrix ℜjk , which elements are the functions of the selected
generalized coordinates.
It is assumed that the system is under the influence of forces of potential and non-potential
nature.
In accordance with the configuration of the mechanical system, the number of bodies N in
the system is defined, and for each body Bk (k = 1, N ) are specified:
• j is the number of body Bj , to which body Bk is connected; or the coordinate system
number Σj with respect to which the position of the body Bk is determined ( j < k for k > 1
and j = 0 for the body B1 ),
• the numbers of rotation axes and its corresponding identifiers of rotation angles,
• rjOk is the radius vector of point Ok of connection of bodies Bj and Bk in the axes Σj ,
• rkCk is the radius vector of the mass center of the body Bk in the axes Σk ,
j
• VO k
is the vector of relative linear velocity of point Ok in the projections onto axes Σj
(for the body B1 is the vector of absolute linear velocity),
• the body mass mk and its tensor of inertia ΘOk relative to the point Ok ,
• the list of generalized coordinates.
In the integrated development environment Embarcadero Delphi using the programming
language Object Pascal, a graphical editor was created. The software was created within the
classical paradigm of object-oriented programming using an incremental development model.
The editor is intended for the formation of a symbolic description of the studied mechanical
system using the Lagrange formalism. As an example, we will consider the description of a
satellite with a stabilizer in the paragraph 3.
The user interface of the editor allows
1) Calling the dialog box for a general description of the mechanical system and filling in the
form, indicating, for example, the type of potential forces influencing the system,
2) Defining the structure of the interconnection of bodies (system configuration) using editor
tools,
3) Assigning semantic data (geometric and kinematic characteristics) to the bodies. In order
to set the semantic data, the context menu is used. These characteristics include, for example,
the radius vectors of connection points of the bodies and the vector of relative linear velocities
of these points; the radius vectors of the centers of mass; the numbers of the rotation axes and
their corresponding rotation angle identifiers; and a list of generalized coordinates,
4) Generating automatically a symbolic description of the investigated mechanical system.
The final description is saved as a text file and served as an input file for another software
package for modeling and qualitative analysis of mechanical systems in a symbolic form.
3. Description of a satellite with a gravitational stabilizer
The mass center (point O1 ) of the system moves along the circular Keplerian orbit of radius
R with constant angular velocity. The stabilizer is a rigid rod with point mass at its free end.
The rod is connected to the satellite at point O3 with a 2-degree-of-freedom suspension (see
Figure 2). Here l > 0 is the rod length; r ≥ 0 is the distance from the system mass center to
the point of attachment of the rod. The system is influenced by a gravitation moment.
So, for the system under scrutiny (Figure 2), input data assigning the geometric and kinematic
characteristics have the following form:
The number of bodies is 4 .
Body 1 (orbital coordinate system O1 x1 y1 z1 ) is connected to the body numbered as 0 (inertial
coordinate system Σ0 or Ōx0 y0 z0 ). The mass and the inertia tensor are zero.
The numbers of the axes of rotation are { 3 , 0, 0 } and angles of rotations {−χ , 0 , 0 }, i.e.
the orientation of the orbital coordinate system with respect to the inertial coordinate system
is determined by one rotation at the angle −χ about the third axis Ōz0 . Zeros among the
numbers of the rotation axes indicate that there are no rotations about other axes.
The radius vector of the point O1 in the coordinate system Σ0 is r0O1 = ( R sin χ , R cos χ , 0 ) .
The radius vector of the body mass center is r1C1 = ( 0 , 0 , 0 ) . The vector of relative linear
velocity of point O1 is VO 0 = ( R ω cos χ , −R ω sin χ , 0 ) , where ω = χ̇ = const is the
1 0 0 0
module of orbital angular velocity.
Body 2 (satellite) is described with respect to Body 1. The body mass is M2 .
The numbers of rotation axes are { 2 , 3 , 1 } and the respective angles of rotations are { ψ , θ , ϕ },
i.e. the orientation of the satellite with respect to the orbital coordinate system is determined
by the sequence of three rotations: about the second axis O1 y1 at an angle ψ ; about the third
axis O1 z1 at the angle θ ; about the first axis O1 x1 at the angle ϕ . The radius vector of the
point O2 with respect to the orbital coordinate system is r1O2 = ( 0, 0, 0 ), i.e. points O1 and O2
coincide. The radius vector of the mass center of satellite is r2C2 = ( 0 , 0 , 0 ), i.e. the satellite
mass center is at point O1 . The vector of relative linear velocity of point O2 is VO 1 = (0,0,0).
2
Figure 2. The satellite with a stabilizer.
Jx 0 0
The body’s tensor of inertia at point O2 is written as ΘO2 = 0 Jy 0 , i.e. Jx , Jy , Jz
0 0 Jz
are satellite principal inertia moments.
Body 3 (rod) is connected to Body 2. The numbers of rotation axes are { 3 , 1 , 0 } and the
angles of rotations are { δ , σ , 0 }, i.e. orientation of the body with respect to the coordinate
system O1 x2 y2 z2 associated with the satellite is determined by the sequence of two rotations:
about the third axis at an angle δ and about the first axis at an angle σ .
The radius vector of the point O3 of the connection of the bodies is r2O3 = ( 0 , r , 0 ) . The
radius vector of the mass center is r3C3 = ( 0 , l/2 , 0 ). The vector of relative linear velocity of
the point O3 is VO 2 = (0,0,0).
3
ml2 /3 0 0
The body inertia tensor at point O3 is written as ΘO3 = 0 0 0 , where m is
0 0 ml2 /3
the mass of rod.
Body 4 (the point of the mass m0 on rod) is connected to Body 3. The point inertia tensor
is zero. The radius vector of the point O4 of the joint of the bodies, its relative linear velocity,
and the radius vector of the center of mass, respectively, have the following form:
r3O4 = ( 0 , l , 0 ) ; 3
VO 4
= r4C4 = ( 0 , 0 , 0 ) .
Thus, the described conservative mechanical system has five degrees of freedom (aircraft angles
ψ , θ , ϕ as well as rod rotation angles δ , σ ). It is located in the Newtonian field of gravitation
to the center of the circle Ō.
We have to emphasize that the description of the mechanical system is formed in interactive
mode, when the user inputs the above mentioned data into corresponding windows of the editor
(see, for example, Figure 3). This important process is conducted just once. Later, it is always
possible to correct input data stored in a file (see Figure 4).
Figure 3. The geometric and kinematic
characteristics of the second body.
4. Construction of a symbolical model for a satellite with a stabilizer
While constructing the mathematical model for complex mechanical objects, searching for its
solutions and conducting the qualitative analysis of their properties, it is necessary to operate
with bulky analytical expressions. Thus, it is relevant to apply computer algebra systems and
develop specialized software on the basis of these systems. Another software tool used in this
paper is designed to model the systems of interconnected absolutely rigid bodies, and also to
examine the questions of stability and stabilization of solutions of linearized models on the basis
of classical theorems of motion stability. The software package is a set of interactive programs
executed in the interpretation mode in the CAS “Mathematica” environment.
By the construction of a symbolical model, one implies the obtaining of nonlinear and
linearized differential equations of motion of mechanical systems in symbolic form in computer
memory.
Let us present a list of problems solved using the software package for the mechanical system
described in paragraph 3.
Problem 1. The rotation axis numbers and the corresponding rotation angle identifiers for
each body Bk are used for automatic construction of the
⊲ Matrix of directional cosines ℜjk which determines the angular position of CS Σk relative
to Σj ,
⊲ Vector of the relative angular velocity ωkj of the k-th body in CS Σj and vector of the
absolute angular velocity ωk = ωj + ωkj of the body ,
⊲ Absolute linear velocity vector VOk = VOj + [ ωk × rjOk ] + VO j
k
of point Ok .
For example, we write out the matrix of directional cosines that determines the angular
Figure 4. The final symbolic description of a satellite with a stabilizer.
position of the coordinate system O1 x2 y2 z2 with respect to the coordinate system O1 x1 y1 z1 :
cos θ cos ψ sin θ − cos θ sin ψ
12
ℜ = sin ϕ sin ψ − cos ϕ cos ψ sin θ cos θ cos ϕ cos ψ sin ϕ + cos ϕ sin θ sin ψ
cos ψ sin θ sin ϕ + cos ϕ sin ψ − cos θ sin ϕ cos ϕ cos ψ − sin θ sin ϕ sin ψ
and the projections of the vector of satellite absolute angular velocity ω2 (of the second body)
onto the axes O1 x2 y2 z2 connected to the body
ωx = ϕ̇ + ψ̇ sin θ + ω0 cos θ sin ψ
ωy = θ̇ sin ϕ + ψ̇ cos θ cos ϕ − ω0 (sin ϕ cos ψ + sin θ cos ϕ sin ψ)
ωz = θ̇ cos ϕ − ψ̇ cos θ sin ϕ + ω0 (sin θ sin ϕ sin ψ − cos ϕ cos ψ) .
Before executing Problem 1, the generalized coordinates are replaced with the corresponding
Greek letters, as well as other necessary substitutions are made.
Problem 2. The kinetic energy for the system of bodies as a quadratic form with respect
to generalized velocities
5 X 5 5
1X X
T (q, q̇) = Aij (q) q̇i q̇j + Ai (q) q̇i + A0 (q) (1)
2 i=1 j=1 i=1
and the force function U (q) of the approximate Newtonian field of gravitation have been
obtained. Here q = (ψ , θ , ϕ , δ , σ) is the vector of generalized coordinates. It is established
that the rotation angle χ is a cyclic generalized coordinate, i.e. the system’s Lagrangian
L = T (q, q̇) + U (q) does not depend on this coordinate. The formulas used in solving Problems
2 and 3 are in [4].
Problem 3. Five nonlinear equations of motion in Lagrange form of the 2nd kind have
been constructed. Each of equations takes up several screens of the monitor. Therefore, it is
impossible to present them in the paper. As a confirmation, let’s demonstrate, for example, the
terms with respect to the accelerations (the notation q ′′ , see Figure 5) in the nonlinear equation
with respect to the θ coordinate.
1
l H- l Cos@∆D Sin@2 ΣD Sin@jD Hm + 3 m0 L +
6
Cos@ΣD Cos@jD H3 r Cos@∆D Hm + 2 m0 L + 2 l Cos@ΣD Hm + 3 m0 LLL ∆¢¢ +
I2 m II3 r2 + l Cos@ΣD H3 r Cos@∆D + l Cos@ΣDLM Cos@jD2 - l H3 r + 2 l Cos@∆D Cos@ΣDL
Cos@jD Sin@ΣD Sin@jD + l 2 ISin@∆D2 + Cos@∆D2 Sin@ΣD2 M Sin@jD2 M +
3 Cos@jD2 I2 Jz + Il 2 + 2 r2 + 4 l r Cos@∆D Cos@ΣD + l 2 Cos@2 ΣDM m0 M +
6 Sin@jD ISin@jD Jy + l I- 2 Hr + l Cos@∆D Cos@ΣDL Cos@jD Sin@ΣD +
l ISin@∆D2 + Cos@∆D2 Sin@ΣD2 M Sin@jDM m0 MM Θ¢¢ +
l Sin@∆D HH- 3 r Cos@jD Sin@ΣD Hm + 2 m0 L + 2 l Sin@jD Hm + 3 m0 LL Σ¢¢ +
Hl m Cos@jD Sin@2 ΣD + m Cos@ΣD H3 r + 2 l Cos@∆D Cos@ΣDL Sin@jD +
3 Hl Cos@jD Sin@2 ΣD + 2 Cos@ΣD Hr + l Cos@∆D Cos@ΣDL Sin@jDL m0 L j¢¢ L -
1
I4 l m I3 r Cos@∆D2 Cos@ΘD Cos@2 jD Sin@ΣD + 3 r Cos@ΘD Cos@2 jD Sin@∆D2 Sin@ΣD +
4
l Cos@∆D Cos@ΘD Cos@2 jD Sin@2 ΣD - l Cos@jD Sin@∆D Sin@ΘD Sin@2 ΣD -
Cos@ΣD H3 r + 2 l Cos@∆D Cos@ΣDL Sin@∆D Sin@ΘD Sin@jDM +
m Cos@ΘD I- l 2 + 12 r2 + 12 l r Cos@∆D Cos@ΣD + l 2 HCos@2 ∆D + H3 + Cos@2 ∆DL Cos@2 ΣDLM
Sin@2 jD + 12 Cos@ΘD Sin@2 jD I- Jy + Jz M -
6 I- 2 l 2 Cos@∆D ICos@ΘD Cos@2 jD Sin@2 ΣD - 2 Cos@ΣD2 Sin@∆D Sin@ΘD Sin@jDM + 4 Cos@∆D2
Cos@ΘD Hl Cos@jD Sin@ΣD + r Sin@jDL H- r Cos@jD + l Sin@ΣD Sin@jDL + 4 l Sin@∆D
H- r Cos@ΘD Cos@2 jD Sin@∆D Sin@ΣD + Cos@ΣD Sin@ΘD Hl Cos@jD Sin@ΣD + r Sin@jDLL +
Cos@ΘD II- l 2 + r2 M Cos@2 ∆D - r Hr + 4 l Cos@∆D Cos@ΣDL - l 2 Cos@2 ΣDM
Sin@2 jDM m0 M Ψ¢¢ + ...
Figure 5. Part of the equation of motion with respect to the θ coordinate.
Problem 4. Let us assign unperturbed motion in the form
ψ = 0 , θ = 0 , ϕ = 0 , δ = 0 , σ = 0 , ψ̇ = 0 , θ̇ = 0 , ϕ̇ = 0 , δ̇ = 0 , σ̇ = 0 . (2)
When moving unperturbed, the system principal central axes of inertia coincide with the
axes of orbital coordinate system, and rod is oriented along the radius of the orbit. This is the
equilibrium position of a satellite with a stabilizer in regard to orbital coordinate system.
The given unperturbed motion is automatically substituted into nonlinear equations of
motion. If motion (2) is the solution, then each equation of motion turns into an identity.
Otherwise, the program issues the conditions of existence of a motion, which are to be taken
into account by the user in any further computations or which fulfillment has to be ensured.
Problem 5. Constructed are equations of perturbed motion in the first approximation.
Linearized in the vicinity of the equilibrium position (2), equations of motion for a satellite with
a stabilizer are decomposed into two subsystems. Respectively, a “pitch” subsystem ( θ ) and a
“yaw-and-roll” subsystem ( ψ , ϕ ) are
(
M1 q̈1 + K1 q1 = 0
(3)
M2 q̈2 + G q̇2 + K2 q2 = 0 ,
where all derivatives are calculated on dimensionless time τ = ω0 t ;
ψ a 0 0
θ c f b−a f
q1 = ; q2 = ϕ ; M1 = ; K1 = 3 ; M2 = 0 b f ;
δ f d f f
σ 0 f d
c−b 0 0 0 c−b−a 0
K2 = 0 4(c − a) 4f ; G= a+b−c 0 0 .
0 4f 3f + d 0 0 0
Here, we introduce the following notations:
1
a = Jy ; b = Jx + m r (l + r) + m l2 + m0 (l + r)2 ; c = b + Jz − Jx ;
3
m m
d= + m0 l2 ; f= + m0 rl + d ; c − b − a = Jz − Jx − Jy .
3 2
We can also derive linear equations of motion from the Lagrangian of the system which is
expanded in a Taylor series in the vicinity of the solution up to the 2nd order of smallness with
respect to deviations of coordinates ∆q = q − q o from their values in the unperturbed motion.
Namely, taking into account the Lagrangian structure (1), the following coefficients are expanded
in a series and calculated in the vicinity of equilibrium (2): those with the second degrees of
generalized velocities Aij (q) are up to a constant; those with the first degrees of velocities
Ai (q) are up to linear values with respect to deviations; those with zero degrees of velocities
A0 (q) + U (q) are up to quadratic components with respect to deviations. This procedure has
been called “quadratization” of the Lagrangian.
In the process of simplification of the expression of “quadratized” Lagrangian we have taken
into account the condition of circular orbit ν/R3 = ω0 2 , where, we have to recall, R is the orbit’s
radius, ω0 is a constant value of orbital angular velocity and ν is the gravitation constant.
5. Necessary conditions of stability of equilibrium position
Equations (3) may be interpreted as equations of oscillations of a mechanical system influenced
by potential (with the matrices K1 , K2 ) and gyroscopic (with the matrix G) forces. These forces
are determined by gravitation forces as well as by orbital motion. The matrices M1 , M2 play
the role of diagonal blocks of a positive definite matrix of kinetic energy.
We introduce the following four dimensionless parameters:
c−b Jz − Jx b−a d f
α= = ; γ= ; p1 = ; p2 = . (4)
a Jy c f c
The physically obtainable values of the parameters are located within the intervals:
−1 < α < 1 , 0 < γ < 1 , 0 < p1 ≤ 1 , 0 < p2 < 1 .
It is obvious that the characteristic equation of system (3) is factorized as Λ(λ) ≡ Λ(1) Λ(2) = 0 .
After performing elementary transformations with the characteristic matrices (multiplying their
rows by positive factors), we obtain the characteristic determinants in notation (4), respectively,
in the “pitch” subsystem and in the “yaw-and-roll” subsystem
λ2 + 3γ p2 (λ2 + 3)
Λ(1) = = λ4 (p1 − p2 ) + 3 λ2 (1 + γp1 − 2p2 ) + 9 (γ − p2 ) ;
λ2 + 3 λ2 p1 + 3
λ2 + α λ(α − 1) 0
(2)
Λ = λ(α − 1)(γ − 1) λ2 (1 + γα) + 4(α + γ) (λ2 + 4)p2 (α + 1) =
0 λ2 + 4 λ2 p1 + (3 + p1 )
= v 6 λ6 + v 4 λ4 + v 2 λ2 + v 0 , where v6 ≡ det M2 = (1 + γ α) p1 − (α + 1)p2 ;
v4 = 3 (1 + αγ) + (α + 1) (p1 (α + 3γ + 2) − (α + 8)p2 ) ;
v2 = (1 + 3γ + α (α + 2γ + 3))(3 + p1 ) + 4 α p1 (γ + α) − 8 p2 (α + 1)(α + 2) ;
v0 ≡ det K2 = 4 α((α + γ)(3 + p1 ) − 4 (α + 1) p2 ) .
The characteristic equation of system (3) contains λ only in even degrees. The stability of a
trivial solution takes place when all roots with respect to λ2 , being simple, will be real negative
numbers. The algebraic conditions providing specified properties of roots (necessary conditions
of stability), represent the system of inequalities
p1 − p2 > 0 , 1 + γ p1 − 2 p2 > 0 , γ − p2 > 0 ,
(5)
p2 γ 2 + ( 4 p (1 − p ) − 2 p ) γ + ( 1 − 4 p (1 − p ) ) > 0 ;
1 2 1 1 2 1
v6 > 0 , v4 > 0 , v2 > 0 , v0 > 0 ,
(6)
Dis ≡ v 2 v 2 − 4 v 3 v − 4 v 3 v + 18 v v v v − 27 v 2 v 2 > 0 .
4 2 2 6 4 0 6 4 2 0 0 6
It is worth noting that the first conditions in (5), (6) is always satisfied by virtue of the positive
definiteness of the kinetic energy matrix. The latter inequality in (6) has also been obtained in
analytical form and represented in the form of an explicit dependence on the parameters (4).
This polynomial is not presented due to being immense.
Emphasize that the calculation of the coefficients vi (i = 0, 3) and the discriminant Dis from
(6) was also performed using the software package.
A parametric analysis of the stability conditions (5), (6) has been conducted in [10]. This
article also considers the possibility of gyroscopic stabilization of an unstable equilibrium
position. In contrast to the presented passive stabilization system, the possibility of active
control of a gravitational stabilizer is investigated in [11]. Assuming an instability of a potential
system, this article considers the problem of stabilization of the system up to asymptotic
Lyapunov stability.
6. Conclusion
The software allows to automatize, and, consequently, essentially speed up the processes of
modeling and dynamic analysis of complicated systems and to avoid errors at all stages of
research. The effectiveness of the presented software tools is confirmed, in particular, by their
use in the study in symbolic form of (i) the dynamics of a satellite with gyrodines, (ii) the
stability of the relative equilibriums of the gyrostat in a circular orbit, (iii) obtaining invariant
manifolds and analyzing their stability in the Clebsch-Tisseran-Brun problem [12].
Acknowledgments
The work was supported by the Russian Foundation for Basic Research (grant No. 19-01-00301).
References
[1] Buchberger B, Collins G E, Loos R and Albrecht R 1983 Computer Algebra: Symbolic and Algebraic
Computation (New York: Springer-Verlag)
[2] Davenport J H, Siret Y and Tournier E 1988 Computer Algebra: Systems and Algorithms for Algebraic
Computations (New York: Academic Press)
[3] Cohen J S 2003 Computer Algebra and Symbolic Computation: mathematical methods (A K Peters, Ltd.)
[4] Banshchikov A V, Bourlakova L A, Ivanova G N, Irtegov V D, Novickov M A and Titorenko T N 1998 Experi-
ence of development and usage of packages of symbolic computations intended for investigation of mechan-
ical systems Elect. Proc. of the 4th Int. IMACS Conf. on Applications of Computer Algebra (Prague: Czech
Technical University) URL https://math.unm.edu/ACA/1998/sessions/history/bourlakova/index.html
[5] Wolfram S 1999 The Mathematica Book. Fourth Edition (New York: Cambridge University Press)
[6] Banshchikov A V and Bourlakova L A 1997 Algorithms of Symbolic Computation Used in Stability Analysis
Programming and Computer Software 23(3) 173–179
[7] Irtegov V D and Titorenko T N 2001 Using the system “Mathematica” in problems of mechanics Mathematics
and Computers in Simulation 57(3) 227–237
[8] Gutnik S A, Guerman A and Sarychev V A 2015 Application of computer algebra methods to investigation
of influence of constant torque on stationary motions of satellite Lecture Notes in Computer Science 9301
198–209 Springer
[9] Prokopenya A N, Minglibayev M Z and Mayemerova G M 2014 Symbolic calculations in studying the problem
of three bodies with variable masses Programming and Computer Software 40(2) 79–85
[10] Banshchikov A V 2002 Parametric analysis of stability conditions for a satellite with a gravitation stabilizer
Proc. of the 5th Int. Workshop on Computer Algebra in Scientific Computing ed V Ganzha et al. (Munich:
Technische Universität München) pp 1–6
[11] Banshchikov A V 2017 On the Asymptotic Stability of a Satellite with a Gravitational Stabilizer Lecture
Notes in Computer Science 10490 16–26 Springer
[12] Irtegov V D and Titorenko T N 2012 Invariant manifolds in the Clebsch-Tisserand-Brun problem Journal of
Applied Mathematics and Mechanics 76(3) 268–274