Using Mathematica System in Mechanics Problems.

Irtegov V.D.,Titorenko T.N.

Full paper in compressed Postscript *.ps.gz

Mechanics as a science used computer-aided symbolic calculations since their appearance and indirectly helped to develop computer algebra in its modern form. It is possible to think that mutual interest of these disciplines is partially explained by the fact that Newton's dynamics for long period of its development offered a huge diversity of problems, methods and algorithms. The fields of authors interest presented in this paper are qualitative theory of ordinary differential equations, stability of motion and dynamics of solid body. Using the methods of qualitative analysis requires constructing and possible transformation of motion equations of complex mechanical systems in the Lagrange, Hamilton, Euler-Lagrange, Routh and other forms. It is labour-intensive but well algorithmizable job.

In many problems of gyroscopia, space dynamics and robotics the representation of motion equations in symbolic form allows to avoid unreasonably great numerical experiments, stimulates effective usage and development of qualitative methods. In these areas computer algebra systems usually offer the significant help both on modelling stage (construction of kinetic energy and force function of mechanical system, construction and transformation of motion equations) and during qualitative analysis of obtained equations. This help is significant even for objects of moderate dimension.

Authors of the report together with the colleagues have developed and use software package based on computer algebra system Mathematica. This package implements modelling of motion of solid bodies systems with holonomic links in Lagrange's formalism; construction of kinetic energy and force functions for certain force fields, construction of differential equations in the Lagrange, Hamilton and Routh forms; qualitative analysis of the constructed equations; construction and analysis of stationary solutions stability by second Lyapunov's method (including the first approximation methods that use the Routh-Hurwitz method and its generalizations); finding of first integrals and investigation of system spectrum, bifurcation sets of stationary manifolds, etc. with the help of these integrals. The algorithms of mechanical systems modelling and qualitative analysis of their motion equations are published in [1,2]

Let us discuss several applications of the above mentioned software package.

1. Modelling of solid body (satellite) with a point of mass m₂attached to main body with visco-elastic spring, moving in Newton's central force field (Fig.1).

Motion of holonomical mechanical system with ${\it n}$ degrees of freedom in potential force field is described by Lagrange's function

L=T+U.

$T=T(q,\dot q)$ is kinetic energy of the system; U(q) is force function of potential force; $q_{1},\ldots,q_{n}$ are generalized coordinates that describe a position of mechanical system in configuration space; ${\dot q_{1},\ldots,\dot q_{n}}$ are derivatives of the coordinates by ${\it t}$.

To calculate Lagrange's function of the mechanical system in Fig.1, we must supply the description of the system geometry to program, as it shown below. We use three coordinate systems for this description: inertial $\Sigma^0$, orbital $\Sigma^1$ and $\Sigma^2$ connected to the point m₂.

Program input.

The description of motion of the coordinate system $\Sigma^1$ relative to $\Sigma^0$:

mass = 0,

r_O1⁰ = $\{R \cos \lambda, R \sin \lambda, 0\}$,

r_C1¹ = $\{0, 0, 0\}$,

v_O1⁰ = $\{-R w_0 \sin \lambda, R w_0 \cos \lambda, 0\}$,

$J^{O_1} = \{\{ 0, 0, 0\}, \{0, 0, 0\}, \{0, 0, 0\}\}$,

the numbers of rotation axes = $\{3, 0, 0\}$,

the rotation angles = $\{\lambda, 0, 0\}$,

the list of coordinates = $\{\lambda\}$.

The description of the satellite motion in the coordinate system $\Sigma^1$:

mass = m₁,

r_O1¹ = $\{0, 0, 0\}$,

r_C1¹ = $\{0, 0, 0\}$,

v_O1¹ = $\{0, 0, 0\}$,

$J^{O_1} = \{\{ A, 0, 0\}, \{0, B, 0\}, \{0, 0, C\}\}$,

the numbers of rotation axes = $\{1, 2, 1\}$,

the rotation angles = $\{\psi, \theta, \varphi\}$,

the list of coordinates = $\{\psi, \theta, \varphi\}$.

The description of the coordinate system $\Sigma^2$ motion relative to $\Sigma^1$:

mass = 0,

r_O2¹ = $\{0, 0, r\}$,

r_C2² = $\{0, 0, 0\}$,

v_O2¹ = $\{0, 0, 0\}$,

$J^{O_2} = \{\{ 0, 0, 0\}, \{0, 0, 0\}, \{0, 0, 0\}\}$,

the numbers of rotation axes = $\{0, 0, 0\}$,

the rotation angles = $\{0, 0, 0\}$,

the list of coordinates = $\{\}$.

The description of motion of the point m₂ in the coordinate system $\Sigma^2$:

mass = m₂,

r_O2² = $\{0, 0, z\}$,

r_C2² = $\{0, 0, 0\}$,

v_O2² = $\{0, 0, \dot z\}$,

$J^{O_2} = \{\{ 0, 0, 0\}, \{0, 0, 0\}, \{0, 0, 0\}\}$,

the numbers of rotation axes = $\{0, 0, 0\}$,

the rotation angles = $\{0, 0, 0\}$,

the list of coordinates = $\{z\}$,

Rayleigh's function = $\frac{1}{2} k_1 \dot z^2$,

the potential energy of the spring = $-\frac{1}{2} k_2 z^2$.

Here r_Oi^j is the radius-vector of attachment point O_i of i-th and j-th bodies in coordinate system $\Sigma^j$; r_Ciⁱ is the radius-vector of mass center of i-th body in coordinate system $\Sigma^i$; v_Oi^j is the relative speed of point O_i in projections on axes of $\Sigma^j$; J^O_i is the inertia tenzor of solid body; O₁ is the inertia center of the satellite; $w_0 = \dot \lambda$ is the angle velocity of mass center O₁; m₁ is the mass of the satellite; A, B, C are the main central inertia moments of the satellite; $R = \vert O_0 \ O_1 \vert$; $r = \vert O_1 \ O_2 \vert$; z is the extension of the spring; k₁, k₂ are the viscosity and the elasticity coefficients of the spring.

Program output.

Lagrange's function:

$$(-(m_2 R w_0 (r + z) \cos \varphi \cos \psi) + m_2 R w_0 (r + z) \cos \theta \sin \varphi \sin \psi - w_0 (A + m_2 r^2 +$$

$$+ 2 m_2 r z + m_2 z^2 ) \cos \psi \sin \theta) \dot \varphi + \frac{ A + m_2 r^2 + 2 m_2 r z + m_2 z^2 }{2} \dot \varphi^2 +$$

$$+ (-(m_2 R w_0 (r + z) \cos \varphi \cos \psi \cos \theta) + m_2 R w_0 (r + z) \sin \varphi \sin \psi - w_0 (A + m_2 r^2 +$$

$$+ 2 m_2 r z + m_2 z^2) \cos \psi \cos \theta \sin \theta + C w_0 \cos^2 \varphi \cos \psi \cos \theta \sin \theta -$$

$$- C w_0 \cos \varphi \sin \varphi \sin \psi \sin \theta + w_... ... m_2 z^2 ) \sin \varphi (\cos \psi \cos \theta \sin \varphi +$$

$$+ \cos \varphi \sin \psi) \sin \theta) \dot \psi + (A + m_2 r^2 + 2 m_2 r z + m_2 z^2 ) \cos \theta \dot \varphi \dot \psi +$$

$$+ \frac{ A \cos^2 \theta + C \cos^2 \varphi \sin^2 \theta + B... ...cos^2 \theta + \sin^2 \varphi \sin^2 \theta)}{2} \dot \psi^2 +$$

$$+ (-(C w_0 \cos \varphi \cos \psi \cos \theta \sin \varphi) +... ...w_0 \sin^2 \varphi \sin \psi + w_0 (B + m_2 r^2 + 2 m_2 r z +$$

$$+ m_2 z^2 ) \cos \varphi (\cos \psi \cos \theta \sin \varphi ... ... w_0 (r + z) \cos \varphi \sin \psi \sin \theta) \dot \theta +$$

$$+ (B - C + m_2 r^2 + 2 m_2 r z + m_2 z^2 ) \cos \varphi \sin \varphi \sin \theta \dot \psi \dot \theta +$$

$$+ \frac{(B + m_2 r^2 + 2 m_2 r z + m_2 z^2 ) \cos^2 \varphi + C \sin^2 \varphi}{2} \dot \theta^2 -$$

$$- m_2 R w_0 (\cos \psi \sin \varphi + \cos \varphi \cos \theta \sin \psi) \dot z + \frac{m2}{2} \dot z^2 +$$

$$\frac{A \nu + B \nu + C \nu + 2 m_1 \nu R^2 + 2 m_2 \nu R^2 +... ...- k_2 R^3 z^2 }{ 2 R^3} - \frac{3}{2} A w_0^2 \cos^2 \theta -$$

$$- \frac{ m_2 w_0^2 (r + z)^2 ( \cos^2 \varphi \cos^2 \theta +... ...^2 \cos \varphi \cos \psi \cos \theta \sin \varphi \sin \psi +$$

$$+ \frac{B w_0^2 (cos^2 \psi \cos^2 \theta \sin^2 \varphi + co... ... \cos^2 \psi \cos^2 \theta + \sin^2 \varphi \sin^2 \psi)}{2} +$$

$$+ \frac{m_2 (-\nu + R^3 w_0^2 ) (r + z) \cos \varphi \sin \th... ...i \sin^2 \theta + \frac{A w_0^2 \cos^2 \psi \sin^2 \theta}{2}-$$

$$- \frac{3}{2} B w_0^2 \sin^2 \varphi \sin^2 \theta + m_2 w_0... ...\theta \sin \varphi \sin \psi + \cos \varphi \sin^2 \theta) +$$

$$+ \frac{ m_2 w_0^2 (r + z)^2 (\cos^2 \psi \cos^2 \theta \sin^... ... \cos^2 \varphi \sin^2 \psi + \cos^2 \psi \sin^2 \theta )}{2}.$$

Rayleigh's function: $\frac{1}{2} k_1 \dot z^2$

The list of generalized coordinates: $\{\psi, \theta, \varphi, z\}$

Elapsed time: 0 h. 0 min. 43:1 sec.

2. Modelling of RLC-circuits.

Authors tried to apply the methods of modelling and analysis of mechanical systems to electrical ones. The results of computer-aided modelling and analysis of linear circuits are presented in [3]. Here nonlinear circuits are discussed.

According to [5,6], the state of nonlinear RLC-circuit can be described by the equations system:

$$L_{k}(i_{k}) \frac{di_{k}}{dt} = \frac{\partial P}{\partial i_{k}}, \ \ \ (k=1,\ldots,r)$$

$$C_{l}(v_{l}) \frac{dv_{l}}{dt} = -\frac{\partial P}{\partial v_{l}}, \ \ \ (l=r+1,\ldots,s)$$

where L_k(i_k) is the inductance; C_l(v_l) is the capacitance; $i=(i_1,\ldots,i_r)$ denotes the vector of currents; $v=(v_{r+1},\ldots,v_{r+s})$ denotes the vector of voltages; P(i,v) is the function called the mixed potential that is defined like

$$P(i,v) = \int \limits_{\Gamma} \sum_{\mu>r+s} v_{\mu} di_{\mu} + \sum_{l=r+1}^{r+s} i_{l} v_{l} \vert_{\Gamma}.$$

$\Gamma$ is a curve in the space of currents and voltages of circuit where Kirchhoff's laws are correct.

The function P(i,v) can be used not just for the construction of circuit equations, but also for qualitative analysis of them. In latter case it could be treated as an analogue for Lyapunov's function.

In [5] it is discussed the construction of the mixed potential for the class of complete circuits. The complete means that a circuit is described by the complete set of variables. The complete set of variables is chosen independently without leading to a violation of Kirchhoff's laws and determines current or voltage in any branch of circuit.

The set of currents in inductors and voltages across capacitors is complete. If a circuit is not complete then it is transformed to the complete one by adding inductors in series and capacitors in parallel. The original circuit is considered as a limiting case of the new one.

Let us describe the procedure of constructing the mixed potential for complete circuits. It consists of the following steps.

Choose in the graph a maximal tree $\tau$ with the set of links $\alpha$.

Choose in $\tau$ a subtree $\grave \tau$ with the set of links $\grave \alpha$. $\grave \alpha$ is all branches which connect two nodes of $\grave \tau$ and form a loop with branches of $\grave \tau$ only.

According to [5], the currents $i_1,\ldots,i_r$ in $\alpha - \grave \alpha$ and the voltages $v_{r+1},\ldots,v_{r+s}$in $\grave \tau$ form the complete set of variables. This is correct also for r=0,s>0; r>0,s=0; rs>0.

Calculate the mixed potential by the formula:

$$P(i^*, v^*) = F(i^*) - G(v^*) + (i^*, \gamma v^*), \eqno(1)$$

$$F(i^*) = \sum_{\mu > r+s, \Omega_i} \int \limits_{\Gamma} v_{... ...sum_{\Lambda_{\rho} \cap \grave \tau} \pm v_{\sigma}, \eqno(2)$$

where $i^* = (i_1,\ldots,i_r)$; $v^* = (v_{r+1},\ldots,v_{r+s})$; $\gamma = (\gamma_{\rho \sigma})$ is matrix $r \times s$; $\gamma_{\rho \sigma}$ = +1,-1,0;

$\Omega_i$ is the set of branches where the currents can be determined from $i_1,\ldots,i_r$ by Kirchhoff's node law; $\Omega_v$ is the set of branches where the voltages can be determined from $v_{r+1},\ldots,v_{r+s}$ by Kirchhoff's loop law; $\Lambda_{\rho} \ (\rho \in \alpha - \grave \alpha)$ is a loop; $i_{\rho}$ is a loop current; $v_{\sigma}$ is the voltage of branches of $\grave \tau$ in $\Lambda_{\rho}$.

To illustrate the procedure of constructing the mixed potential we will consider the circuit in Fig.2 and its graph in Fig.3.

The graph of the circuit of Figure 2

Here the nonlinear elements are resistors $f_1(v_5), \ f_2(v_6), \ f_3(v_7), \ f_4(v_8)$.

The maximal tree $\tau$ consists of branches $\{5,6,7,8,9,10,11,12,17\}$, $\alpha$ = $\{1,2,3,4,13,$ $14,15,16\}$, $\grave \tau$ = $\{5,6,7,8\}$, $\grave \alpha$ = $\{13,14,15,16\}$, $\alpha - \grave \alpha$ = $\{1,2,3,4\}$. So the currents i₁,i₂,i₃,i₄ in the inductors L₁,L₂,L₃,L₄and the voltages v₅,v₆,v₇,v₈ across the capacitors C₁,C₂,C₃,C₄are the complete set of variables.

Now we can construct the mixed potential by the formula (1).

$\Omega_i$ = $\{1,2,3,4,9,10,11,12,17\}$. According to (2), F(i₁,i₂,i₃,i₄) is calculated only for the resistive branches $\{9,10,11,12,17\}$.

$$F(i_1,i_2,i_3,i_4) = E_1 i_1 - \frac{1}{2}R_1 i_{1}^2 - \fra... ... i_{2}^2 - \frac{1}{2} R_3 i_{3}^2 - \frac{1}{2} R_4 i_{4}^2.$$

$\Omega_v$ = $\{5,6,7,8,13,14,15,16\}$. We choose the branches $\{13,14,15,16\}$ with the resistors.

$$G(v_5,v_6,v_7,v_8) = - \int \limits _0^{v_5} f_1(v) dv - \in... ...\limits _0^{v_7} f_3(v) dv - \int \limits _0^{v_8} f_4(v) dv.$$

To calculate $(i^*, \gamma v^*)$ we construct the loops $\Lambda_1, \Lambda_2, \Lambda_3, \Lambda_4$ and consider their intersections with $\grave \tau$.

$\Lambda_1$ = $\{1,5,8,17,9\}$, $\Lambda_2$ = $\{2,6,5,10\}$, $\Lambda_3$ = $\{3,7,6,11\}$, $\Lambda_4$ = $\{4,8,7,12\}$. $\Lambda_1 \cap \grave \tau$ = $\{5,8\}$, $\Lambda_2 \cap \grave \tau$ = $\{5,6\}$, $\Lambda_3 \cap \grave \tau$ = $\{6,7\}$, $\Lambda_4 \cap \grave \tau$ = $\{7,8\}$.

According to (2), $(i^*, \gamma v^*)$ = i₁ (v₅ - v₈) + i₂(v₆ - v₅) + i₃( -v₆ - v₇) + i₄ (v₇ + v₈) and the mixed potential

$$P(i^*,v^*) = E_1 i_1 - \frac{1}{2}R_1 i_{1}^2 - \frac{1}{2}R... ...\limits _0^{v_5} f_1(v) dv + \int \limits _0^{v_6} f_2(v) dv +$$

$$+ \int \limits _0^{v_7} f_3(v) dv + \int \limits _0^{v_8} f_4... ... + v_8) + i_2(v_5 - v_6) + i_3( v_6 + v_7) + i_4 (-v_7 - v_8).$$

Ability to apply the analytical methods of mechanical systems to electrical ones allows to formulate a problem of modelling of mechanical systems with circuits. This problem can be applied, for example, to analysis of transient states of mechanical systems after losing stability near the boundaries of stationary motions family. The number of these problems could be increased. It is important that parameters of circuits can be measured easier than parameters of mechanical ones.

3. Finally, we present the results of investigation of one generic mechanical system with first integrals. Let us discuss the system of absolutely solid bodies, connected by spherical and cylindrical hinges to single carrying body, that has a fixed point. The positions of carried bodies we will define as generalized coordinates q₁,q₂,...,q_n. Let us suppose that forces acting on the system are potential and could be determined by the force function U(q₁,q₂,...,q_n). The full mechanical energy of such system looks like

$$2H=\sum_{\alpha=1}^{3} \sum_{\beta=1}^{3} J_{\alpha \beta}(q)... ...1}^{n} c_{ij}(q) \dot q_{i} \dot q_{j} - 2U(q_{1},...q_{n}).$$

Here $\omega_{\alpha} \ ( \alpha =1,2,3 )$ are the projections of angle velocity of carrying body to coordinate axes connected to this body; $q_{i}, \dot q_{i} \ (i=1,\ldots,n)$ are the coordinates and the velocities of carried bodies relative to carrying one and to each other. There is another possible interpretation of given characteristic function [4].

The differential equations of described system motion have the following form:

$$\sum_{j=1}^{n} c_{ij} \ddot q_{j}+ \sum_{\alpha=1}^{3} e_{i \... ...c{ \partial c_{kj}}{ \partial q_{k}} ) \dot q_{k} \dot q_{j} -$$

$$\hskip 2cm - \frac{1}{2} \sum_{ \alpha=1}^{3} \sum_{ \beta=1}... ...\partial U}{ \partial q_{i}} = 0 \hskip 1.5cm (i=1,\ldots,n),$$

$$\sum_{\beta=1}^{3} J_{1 \beta} \dot \omega_{\beta} + \sum_{k... ... J_{1 \beta}} {\partial q_{k}} \omega_{\beta} ) \dot q_{k} +$$

$$+ \sum_{\beta=1}^{3} (J_{3 \beta} \omega_{2} - J_{2 \beta} \... ...{\partial e_{k1}}{\partial q_{j}} \dot q_{j} \dot q_{k} = 0,$$

$$\dot \gamma_{1}= \gamma_{2} \omega_{3}- \gamma_{3} \omega_{2}.$$

Another four equations (not shown here) can be obtained from last two by circular substitution of indices 1,2,3. For this differential equations system it is possible to obtain three first integrals in addition to the integral of the full mechanical energy H=h=const. They are:

the integral of kinetic moment projection to vertical axis:

$$V_{1}= \sum_{\beta=1}^{3}[ \sum_{\alpha_=1}^{3} J_{\beta \al... ...a}+ \sum_{k=1}^{n} e_{k \beta} \dot q_{k}] \gamma_{\beta}=m;$$

the integral of square of modulus of kinetic moment:

$$V_{2}= \sum_{\beta=1}^{3}[ \sum_{\alpha=1}^{3}J_{\beta \alph... ...pha}+ \sum_{k=1}^{n} e_{k \beta} \dot q_{k} ]^{2} = {k}^{2};$$

and the cosine integral:

$$V_{3}= \sum_{k=1}^{3}{ \gamma_{k}}^{2}=1.$$

For the entire class of problems that are described by the above differential equations we can make qualitative analysis of phase space by selecting a set of invariant manifolds that provide stationary values for different first integrals from algebra of first integrals with basis formed by the integrals H,V₁,V₂,V₃. Also we can study stability and bifurcations of selected invariant manifolds of stationary motions using Lyapunov and Poincaré methods. For example, let us discuss bundle of the first integrals

$$K=\frac{1}{2} V_{2}- \lambda_{1} V_{1}+\frac{1}{2} \lambda_{1}^{2} V_{3}, \ \lambda_1 = const.$$

The stationary conditions of function K can be written in the form:

$$\frac{\partial K}{\partial \omega_{i}}= \sum_{j=1}^{3} \varphi_{j} J_{ij}=0, \hskip 3cm$$

$$\frac{ \partial K}{\partial \gamma_{i}}=- \lambda_{1} \varphi_{i}= 0 \hskip 2cm ( i=1,2,3),$$

$$\frac{\partial K}{\partial \dot q_{m}}= \sum_{i=1}^{3} \varphi_{i} e_{mi}=0, \hskip 4cm \eqno(A)$$

$$\frac{\partial K}{\partial q_{m}}= \sum_{i=1}^{3} \varphi_{i}... ...i}}{\partial q_{m}} \dot q_{k})=0 \hskip 1cm (m=1,\ldots,n).$$

Where

$$\varphi_{i}= \sum_{\alpha=1}^{3}J_{i \alpha} \omega_{\alpha} ... ...{ki} \dot q_{k}- \lambda_{1} \gamma_{i} \ \ ( i=1,\ldots,3).$$

According to the theorem of [3] it is possible to state that the elements of manifolds family

$$\varphi_{1}=0, \ \ \varphi_{2}=0, \ \ \varphi_{3}=0$$

with parameter $\lambda_{1}$ that provide extremal value to the first integrals from bundle K with any fixed $\lambda_{1}$form the family of invariant manifolds for the original differential equations system, because the rank of system (A) matrix is equal to 3. This is the consequence of submatrix elements properties:

$$\left( \begin{array}{ccc} J_{11}(q) & J_{12}(q) & J_{13}(q) ... ... \\ J_{31}(q) & J_{32}(q) & J_{33}(q) \end{array} \right),$$

that are interpreted as the inertia moments of mechanical system.

Using this approach we can obtain stability conditions for the elements of found family of invariant manifolds of stationary motions (IMSM) by second Lyapunov's method. Non-writing the equations of disturbed motion in the neighbourhood of our IMSM, let us consider the first integral Kas Lyapunov's function. Let us write it in the terms of deviation from the manifold in question, i.e (in our case) in variables $\varphi_{i}$. The direct calculation gives:

$$K= \frac{1}{2}({ \varphi_{1}}^{2}+{\varphi_{2}}^{2}+{\varphi_{3}}^{2}).$$

Because K is the sign-definite function of variables $\varphi_{1}, \varphi_{2}, \varphi_{3}$, using well known Zubov's theorem on stability of sets we conclude, that all elements of discussed IMSM family are stable with all values of $\lambda_{1} \neq 0$.

References

1. Banshchikov A.V., Bourlakova L.A., Ivanova G.N., Simonov S.A. MECHANIK package of symbolic computations. The problems and the structure //Packages of applied programs. Results and applications. Novosibirsk. Nauka, 1986, p.96-105. (in Russian)

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

3. Irtegov V.D., Titorenko T.N. On modelling and research of some problems with the use of computer algebra //Programming (Programirovanie), 1, 1997, p.68-74. (in Russian)

4. Litvin-Sedoi M.Z. Systems mechanics of connected solid bodies //Results of science and engineering. General mechanics. Moscow.: VINITI, vol.5, 1992, p.3-52. (in Russian)

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

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

7. Gaylord Richard J., Kamin Samuel N., Wellin Paul R. Introduction to Programming with Mathematica. Springer-Verlag, 1996.

8. Wolfram S. Mathematica: A System for Doing Mathematics by Computer. Addison - Wesley Publ. Co., 1988.