Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 UDK - UDC 534.2:004.41 Kratki znanstveni prispevek - Short scientific paper (1.03) Simbolno-številčno analiziranje nihanj sistemov z velikim številom stopenj prostosti A Symbolic-Numeric Vibrations Analysis of Systems with Many Degrees of Freedom Regina Kulvietiene - Genadijus Kulvietis - Inga Tumasoniene (Vilnius Gediminas Technical University, Vilnius) Za analizo nihanja sistemov z velikim številom prostostnih stopenj smo uporabili metode računalniške algebre. Z vidika računalniške algebre smo primerjali dve metodi reševanja in izbrali metodo harmonskega ravnovesja. Sistem smo razdelili v linearni in nelinearni del. Linearni del lahko popišemo na običajen način. Za rešitev nelinearnega dela v zaključeni obliki pa smo uporabili simbolični izračun. Izbrani simbolno-števični pristop ima veliko prednosti, še posebej pri sistemih z velikim številom stopenj prostosti: vodi v poenostavitev teoretičnih oblik modelov, opazno zmanjšanje obsežnosti dobljenih enačb in zato tudi skrajšanje časa računanja ter povečanje možnosti uporabe modelov z več objekti v drugih posebnih okoljih. © 2006 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: analize nihanj, metode računalniške algebre, računanje simbolično, stanja ustaljena) Computer algebra techniques were applied to analyze the vibrations of systems with many degrees of freedom. For this purpose, two solution methods were compared from the computer algebra point of view, and the harmonic balance method was chosen. The system is divided into linear and nonlinear parts. The linear part of the system can be formalized as usual, and symbolic computations were applied to perform a closed-form solution of the nonlinear part. The symbolic-numeric approach chosen, specially dedicated to systems with many degrees of freedom, affords various advantages: it leads to a simplification of the theoretical formulation of the models, a considerable reduction in the size of the generated equations, and hence in the resulting computing time, and also enhanced portability of the multibody models to other specific environments. © 2006 Journal of Mechanical Engineering. All rights reserved. (Keywords: vibration analysis, computer algebra, symbolic-numeric computations, steady state) 0 INTRODUCTION It is difficult to find a branch of modern engineering or science in which vibration problems can be neglected. The theory of oscillations even creates the foundation for many fields in technical sciences, for example, acoustics, radio engineering, automatic control, vibro-isolation. The development and design of new machines belongs to many domains of modern engineering, where the investigation of the level of vibrations is one of the main goals of basic and applied research [14]. The well-elaborated mathematical theory of linear differential equations allows us to solve, with sufficient accuracy, almost all the problems resulting from the fast development of various types of machines, mechanisms, machine tools, robots, traffic means, etc. The increased requirements on the speed, power, long life and reliability of mechanical systems change the situation and prove that the linear models of mechanical systems are only a first approximation of the real process. Using the linear theory, we cannot explain all the phenomena arising during operation. This unexpected behavior is often connected with catastrophic accidents [15]. In view of the accelerating trends of innovation in machine design, in new technologies, new electronic systems and in other branches of science, more complex approaches and more exact solutions are required. This cannot be done without consideration of the very detailed mathematical models respecting the nonlinear characteristics of 309 Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 real physical systems. For this reason it was necessary to use nonlinear models in the form of ordinary differential equations, differential equations with time-varying parameters, with random properties, with delayed arguments, [13] etc. The nonlinear oscillations (or more precisely, the oscillations of nonlinear systems) have been intensively studied for half a century ([7] and [15]). A lot of new methods for solutions were elaborated and much new knowledge about the properties, behavior and applications were presented ([13] and [14]). These include the properties of free and forced vibrations, sub- and ultra-harmonic resonances, self-excited vibrations, the interaction of self-excitation and external excitation, transients, bifurcation and jump phenomena, stability problems, identification procedures, etc. The natural vibration of nonlinear systems is of primary concern in studying resonance phenomena because the backbone curves (the amplitude-frequency relations) and the modes of vibrations, i.e., dynamic characteristics of systems, are determined. Analytical expressions for the backbone curves are very complex and numerical methods are not a convenient way to analyze nonlinear oscillations. In some cases, such as the one in the range where the internal resonance exists, the corresponding backbone curves have a very complex shape owing to the presence of sharp peaks, looping characteristics and rapidly changing slopes. It is difficult to determine these types of backbone curves using the developed numerical methods ([1] and [2]). Simulations by means of numerical methods are powerful tools for the investigation in mechanics; however, they have serious drawbacks, e.g., finite precision, difficulties in determining transient states and steady states, and the investigation of stability is error-prone and complex. The analytical steady-state solution by hand requires a lot of routine work, is error-prone and available only for very simple systems ([3] and [11]). Here, the computerized symbolic manipulation systems – so-called computer algebra – are indispensable tools. Symbolic manipulations provided by computer algebra systems in combination with the high-power number-crunching abilities of traditional hardware and software really opens a new way to the large-scale computations needed in steady-state solutions and stability analyses ([3], [5] and [9]). Subsequently, Gilsinn [12] demonstrated that with the help of a symbolic manipulator package approximations to the invariant tori can be developed by using Galerkin’s variational method. However, the amount of computation is so excessively large that it soon becomes impracticable with the increasing strength of the non-linearity. The aim of this paper is to describe the theoretical background of systematic computer algebra methods for analyzing the free and steady-state periodic vibrations of nonlinear structures. Many analytical steady-state solution methods are developed, but each of them has different capabilities, e.g., small parameter methods give a solution in closed form, and the harmonic balance method only converts nonlinear differential equations to algebraic. On the other hand, it is very important to assess the efficiency of analytical methods in terms of the view of a computer algebra system. For this reason, the VIBRAN computer algebra system [4] was used. The VIBRAN system has been developed by authors and used for engineering problems since 1979. The VIBRAN computer algebra system is a FORTRAN pre-processor for analytical computation with polynomials, rational functions, and trigonometric series. The main distinction between VIBRAN and other computer algebra systems is the representation form of the analytical expressions. The analytical expressions are stored in matrix form and the analytical perturbations are replaced by matrix operations. Analytical operations can be performed by calling VIBRAN procedures directly from FORTRAN code. Figure 1 illustrates the VIBRAN procedure fragment for performing a Fourier transformation. A special VIBRAN procedure can generate an optimized FORTRAN code from the obtained analytical expressions, which can be directly used in the programs for numerical analysis. In this paper the efficiency of two analytical methods is assessed from the point of view of a VIBRAN computer algebra system and a symbolic-numeric approach for vibrations analysis of the multibody systems is presented. 1 REALIZATION OF THE SMALL-PARAMETER METHOD The small-parameter method (Poincare’s method) [7] was developed to solve systems of nonlinear differential equations. Consider the 310 Kulvietiene R. - Kulvietis G. - Tumasoniene I. Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 POLINOM C,D RATIONAL F READ 1, N PRINT 17, N DO 101 I=1,NN... BINT(C,D,F,IAR,I) ... END PROCEDURE BINT(A,B,C,IAR,K) POLINOM A,B,C... FKCI(A,B,0) SCA1(A,2.) ... MXN(A,PI,1)... END Fig. 1. A fragment of the VIBRAN procedures algorithm of the small-parameter method realized by the VIBRAN computer algebra system for systems of nonlinear differential equations. For the sake of clarity, an algorithm of the small-parameter method is presented for one equation below: &x& + k2 x + f (t) = mF(t, x, x& ,m) (1), where x, x&, x&& - displacement, velocity and acceleration; m - small parameter; f(t) - continuous periodical time function with a period 2p; k - constant, not integer in the non-resonant case and integer in the resonant case; F(t, x, x&,m) - integer polynomial function and periodic with respect to t. The solution of Equation (1) according to the small-parameter method [7] can be found in the form: x=x0+mx1+m x2+m x3 (2), where x0 is the solution of Equation (1) without nonlinearities, xi, i = 1,2,3... are unknown functions of t with a period of 2p. To find these functions, the series (2) is substituted into Equation (1) and the coefficients with the corresponding m power are equalized. With respect to the m power, the linear differential equations are obtained: k2xi = Fi,i = 1,2,3,.. (3), where Fi are the integer rational functions of x0 , x1,..., xi-1; x&0 , x&1,..., x&i -1 and a continuous periodical function with respect to t. The solution of Equations (3) can be found in the form: 2k +z aij cos( jt) + bij sin( jt) (4), j=1 k -j where aij and bij are Fourier series coefficients of the function Fi. In the resonant case, where k is an integer and nearly equal or equal to n, this means: n2 -k2 = em where e has a finite value, the zero power solution must be expressed in the form: x0=j0(t) + M0cos(nt) + N0sin(nt) x1=j1(t) + M1cos(nt) + N1sin(nt) and the i-th power solution can be expressed in the same manner. In the general case, M0 and N0 can be found from the equations: 2p jF(t,j0+M0cos(nt) + N0sin(nt), j0 - M0n sin(nt) + N0n cos(nt), 0) sin(nt)dt = 0 jF(t,j0+M0cos(nt) + N0sin(nt ), j0 -M0nsin(nt) + N0ncos(nt),0)cos(nt)dt = 0 where j0 is the zero power solution of Equation (1) in the non-resonant case, excluding resonant harmonics. The other coefficients Mi and Ni can be found from the system of linear algebraic equations: M i -1 2pjfF]cos(nt)-nfF]sin(nt)lcos(nt)dt 8x ) dx N i -1 2p{Fjsin(nt)+n — \cos(nt)\-cos(nt)dt + Simbolno-številčno analiziranje nihanj sistemov - A Symbolic-Numeric Vibrations Analysis of Systems 311 Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 dF 2p I/ 2prdF JFi*cos(nt)dt = 0 cos( nt) Ni-J 0 dx sin(nt)+n dx dF) dx) J cos(nt)\sin(nt)dt + where: dF\ & j +axjj 1 sfA are known periodic functions. The small-parameter method was realized in VIBRAN for both resonant and non-resonant cases in systems of nonlinear differential equations [6]. 2 REALIZATION OF THE HARMONIC BALANCE METHOD The harmonic balance method is probably the oldest analytical method in the theory of nonlinear vibration ([6] and [15]). Consider the following nonlinear differential equation: &x& + f (x,x& ) = F(t) (5), where x,x,x denote displacement, velocity and acceleration; f(x,x) is a nonlinear function, expandable in a Fourier series; F(t) is assumed to be a periodic function: 00 F(t) = A0 + X (Ai cos(iwt) + Bi sin(iwt)) The solution of the above-mentioned Equation (5) can be expressed in a Fourier series in time: OO x(t) = a0 + ~YJ(ai cos(iwt) + bi sin(iwt)) i=1 (6). A nonlinear function is also expanded in the Fourier series in time: OO f(x, x) = a0 + ^ (ai cos(iwt) + bi sin(iwt)) (7), i=1 where the Fourier-series coefficients are calculated using the following formulae: a =— [ f(x,x)dt w 2p 2 w ai= wp\f(x,x)cos(iwt)dt w f bi = f(x, x) sin(iwt )dt p 0 A substitution of the formule (6) and (7) into Equation (5) gives an infinite number of algebraic equations to determine the unknown coefficients of Solution (6) : (8) or i2w2bi = bi - Bi,i = 1,2,3. a0 = a0 (a0 , a1,b1,...) ai = ai (a0,a1,b1,...) bi = bi (a0 , a1,b1,...) The above-mentioned version of the harmonic balance method was realized in VIBRAN for the system of nonlinear equations [5]. 3 EFFICIENCY MEASUREMENTS FOR BOTH METHODS Two VIBRANs programs that realize the above-mentioned methods were tested for the following equation: &x& +c1x& +c2x+c3x2 +c4x3 = d0 + d1 sin(t) + d2 cos(t) This equation describes the dynamics of an aerodynamically supported magnetic head in a recording device ([8] and [10]). We present below the result of the solution for the first harmonics: the first index for the solution coefficients A and B is the equation number and the second one is the harmonics number. The first result corresponds the first equation in Formula (8) and afterwards we present the result for the cosine and sine coefficients of the first term of the Fourier series. The result for the first equation of Formula (8) is: 0=A10*C2-.5*B11**2*C3-.5*A11**2*C3-C3*A10**2+.75*C4*B11**2*A12-1.5*C4*B11**2*A10-1.5*C4*B11*B12*A11-1.5*C4*B12**2*A10-.75*C4*A11**2*A12-1.5*C4*A11**2*A10-1.5*C4*A12**2*A10-C4*A10**3-D0 The result for the cosine coefficient of the first term of the Fourier series is: 312 Kulvietiene R. - Kulvietis G. - Tumasoniene I. Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 1400 1200 1000 800 M 600 400 200 1234 n Fig. 2. Number of terms in the solution’s expression 0=-A11-B11*C1+A11*C3-B11*B12*C3-2*A11*A10*C3-A11*A12*C3-.75*C4*B11**2*A11-3*C4*B11*B12*A10-1.5*C4*B12**2*A11-.75*C4*A11**3-1.5*C4*A11*A12**2-3*C4*A11*A10**2-3*C4*A11*A12*A10-D2 The result for the sine coefficient of the first term of the Fourier series is as follows: 0=-B11-A11*C1+B11*C3-B11*A12*C3-2*B11*A10*C3-A11*B12*C3-.75*C4*B11**3+3*C4*B11*A12*A10-1.5*C4*B12**2*B11-.75*C4*A11**2*B11-1.5*C4*B11*A12**2-3*C4*B11*A10**2-3*C4*A11*B12*A10-D1 In Figure 2 the number of terms in the solution expression (M) for the small-parameter (upper column) and for the harmonic balance method (lower column) with respect to the number of harmonics (n) in the solution (number of power in the case of small parameter method) is presented. In Figure 3 the presented graphic illustrates the convergence of the above-mentioned analytical methods for the four coefficients. The upper curve corresponds to the harmonic balance method and the lower one corresponds to the small-parameter method. In this case, the magnetic head construction parameters were: c1 = 0.1,c2 = 0.228*106,c3 = 0.167*104, c4 = -0.587*103,d0 = 0,d1 = 0.12,d2 = 0.12 The comparison of the above-mentioned analytical methods illustrates the similarities and differences of their application. The similarities illustrate the tolerance curves that are of the same shape. There are 12045 terms in the solution’s expression for the small-parameter method (5 harmonics), whereas only 1524 terms are observed for the harmonic balance method. This means that the harmonic balance method is much more convenient for computer algebra realization, especially for the multibody system, but needs a special stability analysis procedure for the steady-state solution. 5 SYMBOLIC-NUMERIC REALIZATION OF THE HARMONIC BALANCE METHOD Computerized symbolic manipulation is a very attractive means to reliably perform analytical 15 10 5 0 C ks ^ ^ H ^ ----< >~ ,—-- ¦---------------< > 0,2 0,4 0,6 0,8 w Fig. 3. Tolerance for harmonic balance and small-parameter analytical methods 0 Simbolno-številčno analiziranje nihanj sistemov - A Symbolic-Numeric Vibrations Analysis of Systems 313 Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 calculations, even with complex formulae and expressions. But often, a semi-analytical approach, combining the features of analytical and numerical computations, is the most desirable synthesis. This allows the analytical work to be pushed further, before the numerical computations start. For symbolic-numeric computation of the nonlinear oscillation multibody systems the VIBRAN computer algebra system was used. VIBRAN’s special procedure can generate optimized FORTRAN code from obtained analytical expressions, which can be directly used in programs for numerical analysis. Figure 4 illustrates the scheme of the proposed approach. Let us discuss the multibody system with s degrees of freedom. After good known perturbations the equations of motion can be rewritten in the matrix form: [A]{x} + [B]{x} + [C]{x} {H(x,x,t)} + {f(t)} (9), where: { f (t)} = { f0} +{ f1}sin(wt) +{ f2}cos(wt) +... {f(t)} - continuous periodical time function or expandable into a Fourier series, and vector H (x, x&,t) is the result from block 1. The solution of the above-mentioned system can be expressed using the harmonic balance method in the form: {x} ={D0}+{D1}sin(wt) +{D2}cos(wt) +... where {Di} are the unknown vectors that can be found from the nonlinear algebraic equations. Following the harmonic balance method ([6] and [15]) these equations for the first three vectors’ coefficients are: [U]{D} + {f(D)}={Q} (10), where fi are the coefficients of the function f(t) Fourier expansion. Analogously, the equations for other harmonics can be found using VIBRAN’s program, shown in block 4. The expressions of Hi and their required derivatives are expressed in closed form using computer algebra techniques with a FORTRAN code generation procedure. The programs’ structure is shown in Figure 4. Obviously, the matrix of coefficients [U] consists of independent sub-matrix blocks located at the main diagonal. Therefore, the linear part of the matrix equation decomposes into m separate systems (m is the number of harmonics in the solution vector) of size 2s and one system (corresponding to the zero harmonic) of size s. The obtained solution of these systems serves as an initial approximation for further computation. The Newton’s iteration formula claims: {E(Di)} + [E'(Di)]{Di+1 -Di} = { 0 } (11), where {E(D)} is the system of equations and [E’(D)] is the Jacobi matrix of the whole system (including the linear part as well). Substituting these two formulae the Newton iteration becomes ([U] + [f'(Di)]){Di +1} = [f '(Di)]{Di {f(Di)}+{Q} The iterations start with: [U]{A0} = {f} (12). (13). This linear matrix equation must be solved for every iteration step. This equation can be rewritten in the simple form ([U] + [f']){D} = {F} (14). It is not expedient to sum up the matrices [U] and [f] in advance, because they have quite different structures. Matrix [U], as was already mentioned, is block-banded, i.e., consists of separate sub-matrix blocks located at the main diagonal. The structure of matrix [f] depends on the type and location of nonlinear terms in the initial differential equation system, but it is always sparse and can possess nonzero elements far from the main diagonal. Thus, matrix [D] is not stored in the computer memory at all, and its every element is computed by a reference to the special subroutine. Matrix [f’] is stored in the main memory as a sparse matrix. Of course, such storage demands corresponding modifications to the solution algorithm itself. In many applications the solution of the differential equation system must be obtained in a specified domain of some varied parameters (frequency, stiffness, mass, etc.). Therefore, the program is designed in such a way that any parameter of the initial system can be varied with a regular or logarithmic step. Note that the analytical computation is performed only once, while the numerical calculations are repeated every time when the value of any parameter of the system is being changed. 314 Kulvietiene R. - Kulvietis G. - Tumasoniene I. Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 The proposed symbolic-numeric method was tested for the dynamics of a double-head recording device with aerodynamic support [8]. The following equations describe the system behavior: & x&1 +k1x&1 +k2x&2 +k3x1 +u0 +u1x1 +u2x12 +u3x13 = f1(t ) & x&2 +k6x&1 +k7x&2 +k8x1+k9x2 = f2(t ) where: f1(t ) = k4 sint + k5 cost , f2 (t ) = k10 sint + k11 cost Analytical expressions obtained according to the VIBRAN program conclude the part of the analytical calculations. The analytical result for Hi in Equation (5) is: H0=k3*a10+u1*a10+.5*u2*a11**2+.5*u2*b 11**2+u2*a10**2+1.5*u3*a10*a11**2+1.5 *u3*a10*b11**2+u3*a10**3+u0 H1=a20*k9+k8*a10+u1*a10+.5*u2*a11**2 +.5*u2*b11**2+u2*a10**2+1.5*u3*a10*a1 1**2+1.5*u3*a10*b11**2+u3*a10**3+u0 H2=b21*w*k2+w*b11*k1+k3*a11+a11*u1+ 2*u2*a11*a10+.75*u3*a11**3+.75*u3*a11* b11**2+3*u3*a11*a10**2 H3=-a21*w*k2-w*a11*k1+k3*b11+b11*u1+2*u2*b11*a10+ .75*u3*a11**2*b11+.75*u3*b11**3+3*u3* b11*a10**2 H4=b21*w*k7+a21*k9+w*b11*k6+k8*a11+ a11*u1+2*u2*a11*a10+.75*u3*a11**3+.75 *u3*a11*b11**2+3*u3*a11*a10**2 H5=-a21*w*k7+b21*k9-w*a11*k6+k8*b11+b11*u1+2*u2*b11*a10+ .75*u3*a11**2*b11+.75*u3*b11**3+3*u3* b11*a10**2 The corresponding derivatives are very simple and there are only 25 nonzero terms: H23 = k1*w+1.5*u3*b11*a11 H26 = k2*w H31 = 2*u2*b11+6*u3*a10*b11 H32 = -k1*w+1.5*u3*b11*a11 H41 = 2*u2*a11+6*u3*a10*a11 H43 = k6*w+1.5*u3*b11*a11 H45 = k9 H46 = k7*w H51 = 2*u2*b11+6*u3*a10*b11 H52 = -k6*w+1.5*u3*b11*a11 H55 = -k7*w H56 = k9 Figure 4 illustrates the generated program code that could be directly used for numerical investigations. This code is very simple and contains only 89 floating-point product operations. B(2)=+O(10)+O(11)+Y5 ... END___________________________________ Fig. 4. A fragment of the code generated by VIBRAN 5 CONCLUSION On the basis of non-linear differential equations solved by a harmonic balance method and the synthesis of the VIBRAN analytical calculation system an investigation method for non-linear systems was created. This method combines the advantages of analytical calculation methods and computer algebra. They are compared on the principle of a symbolic-numerical calculation, where the analytical rearrangements are applied only to the nonlinear part of the system, and at the same time the linear part of the system could be easily solved in a numerical way. The proposed method provides smaller expressions for the analytical computation and allows the analysis of systems with greater order. SUBROUTINE ISRA01(A,O) IMPLICIT REAL(A-Z) DIMENSION A(20),O(60) u0 =A( 1) u3 =A( 2) u2 =A( 3) k3 =A( 4) ... O(17)=a11*a10*u3 O(19)=b11*a10*u3 O(20)=O(31)*b21 O(21)=O(29)*b1 ... END SUBROUTINE ISRA(A,B,O) DIMENSION A(20),B(42),O(60) Y1=+.5*O(3) Y2=+.5*O(4) Y3=+1. 5*O(6) ... Y26=+.75*O(13) Y27=+2.25*O(1) B(1)=+O(1)+O(2)+O(5)+O(8)+O(9)+Y1+Y2+Y3+Y4 Simbolno-številčno analiziranje nihanj sistemov - A Symbolic-Numeric Vibrations Analysis of Systems 315 Strojniški vestnik - Journal of Mechanical Engineering 52(2006)5, 309-316 The problem considered in this paper clearly de- proach is that it provides both quantitative and quali- monstrates the power of the symbolic-numeric per- tative results regarding the dynamic behavior of turbation method. An important feature of this ap- multibody systems. 6 REFERENCES [I] R. Lewandowski (1997) Computational formulation for periodic vibration of geometrically nonlinear structures, Part 1: Theoretical background, Int. J. Solids Structures, 34 (15) (1997) 1925-1947. [2] E.Riks (1984) Some computational aspects of the stability analysis of nonlinear structures, Computational Methods in Applied Mechanical Engineering, 47 (1984) 219-259. [3] D.M. Klymov, V.M. Rudenko (1989) Metody kompjuternoj algebry v zadačah mehaniki, Nauka, Moscow. [4] R. Kulvietiene, G. Kulvietis (1989) Analytical computation using microcomputers, LUSTI, Vilnius. [5] R. Kulvietiene, G. Kulvietis, J. Galkauskaite (1997) Computer algebra application to large nonlinear oscillation systems. Mechanine technologija, XXV, Kaunas, 126-130. [6] R. Kulvietiene, G. Kulvietis (1996) Symbolic solution of nonlinear oscillation system using harmonic balance method, Proc. of the 2nd European Nonlinear Oscillations Conference, Vol. 2, 1996, 109-112. [7] I.G. Malkin (1956) Some problems of the theory of nonlinear oscillations, Gostehizdat, Moscow. [8] R. Kulvietiene (1982) Dynamics of aerodynamically supported magnetic head in the recording device, Ph.D. Thesis, KTU, Kaunas. [9] F. San-Juan, A.Abad (2001) Algebraic and symbolic manipulation of poisson series. Journal of Symbolic Computation, 32 (5), (2001) 565-572. [10] A. Cepulkauskas, R. Kulvietiene, G. Kulvietis (2003) Computer algebra for analyzing the vibrations of nonlinear structures, Proc. CASA’2003, Lecture Notes in Computer Science, Vol. 2657, Springer-Verlag, Berlin Heidelberg New York, 747-753. [II] J.-C. Samin, P. Fisette (2003) Symbolic modeling of multibody systems, Kluwer Academic Publishers. [12] Gilsinn, David E. (1995) Constructing Galerkin’s approximations of invariant tori using MACSYMA. Nonlinear Dynamics 8 (2) (1995) 269-305. [13] N.N. Bogolyubov, JA. Mitropolskij (1995) Asymptotical methods in the theory of nonlinear oscillations. TTL, Moskva. [14] L. Pust, F Peterka, G. Stepan, G.R. Tomlinson, A. Tondl (1999) Nonlinear oscillations in machines and mechanisms theory. Mechanism and Machine Theory 34 (1999) 1237-1253. [15] C. Hayashi (1964) Nonlinear oscillations in physical systems. McGraw Hill, New York. Authors’ Address: Doc. Dr. Regina Kulvietiene Prof. Hab. Dr. Genadijus Kulvietis Inga Tumasoniene Vilnius Gediminas Technical University Department of Information Technologies rku@fm.vtu.lt gk@fm.vtu.lt tinga@gama.vtu.lt Sauletekio 11 Vilnius 2040, Lithuania Sprejeto: Odprto za diskusijo: 1 leto 23.2.2006 Accepted: Open for discussion: 1 year Kulvietiene R. - Kulvietis G. - Tumasoniene I. Prejeto: 7.11.2005 Received: 316