Solving differential equations of motion for Constrained mechanical systems Текст научной статьи по специальности «Математика»

Теоретическая механика

Differential Equations of Motion for Constrained Mechanical Systems

R. G. Mukharlyamov*, A. W. Beshaw^

* Department of Theoretical Mechanics Peoples' Friendship University of Russia Miklukho-Maklaya, 6, Moscow, Russia, 117198

^ Department of Mathematics Bahr Dar University Ethiopia, Bahr Dar

This paper presents an investigation of modeling and solving system of differential equations in the study of mechanical systems with holonomic constraints. A method is developed for constracting equation of motion for mechanical system with constraints. A technique is developed how to approximate the solution of the problem that is obtained from modeling of kinematic constraint equation which is stable. A perturbation analysis shows that velocity stabilization is the most efficient projection with regard to improvement of the numerical integration. How frequently the numerical solution of the ordinary differential equation should be stabilized is discussed. A procedure is indicated to get approximate solution when the systems of differential equations can't be solved analytically. A new approach is applied for constructing and stabilyzing Runge-Kutta numerical methods. The Runge-Kutta numerical methods are reformulated in a new approach. Not only the technique of formulation but also the test developed for its stability is new.Finally an example is presented not only to demonstrate how the stability of the solution depends on the variation of the factor but also how to find an approximate solution of the problem using numerical integration.

Key words and phrases: numerical integration, kinematic constraint, stable solution, Taylor expansion, row decomposition.

UDC 531.3 Solving

1. Introduction

The theory of mechanical systems with kinematic constraints goes back to the last century, with important contributions by Herzt (1894), Ferrers (1871), Neumann (1888), Vierkandt (1892), Chaplygin (1897). Several recent papers [1] show a strong renewal of interest in that theory, in relation with new developments in control theory.

In the construction of equations of motion, one can treat the constraints imposed on a mechanical system as servo-constraints, and the constraint reactions can be treated as the corresponding controls. Then the construction of equations of motion can be reduced to finding expressions for the constraint reactions on the right-hand sides guaranteeing that the solutions of the system satisfy the constraint equations with the desired accuracy [2]. This is, in a sense, the constraint stabilization problem, and to solve it, one should take account of not only the deviations from the constraint equations but also their derivatives [2]. It turns out that the constraint stabilization is possible only if the constraint equations are particular integrals of the equations of motion of the mechanical system and the differential equations for the deviations from the constraint equations have an asymptotically stable trivial solution [2]. Still this does not guarantee that the deviations from the constraint equations remain bounded in the numerical solution of the differential-algebraic equations comprised by the equations of motion and the constraint equations. In the present paper, to estimate the deviations from the constraint equations, we introduce the equations of constraints and the equations of constraint perturbations and define the notions of stable and asymptotically stable constraints.

2. Modeling of Kinematic Mechanical System

Kinematic state of mechanical system can be described by ordinary differential equations:

y = v(y,t), where y(to) = y0, y G Rn. (1)

Assume the constraints to be imposed:

f (y,t) = 0, where f G Rm, m < n, f (y0,to) = 0. (2)

The vector on the right hand side of equation (1) should be taken so that for all t > t0 and q = q(t) satisfying equation (2) it also satisfies the following equality:

fy V + ft = 0, (3)

where fy = (fij), ft = (fi t), fij = fy, i = 1,2,...,m, j = 1, 2,...,n. The solution of

equation (3) can be constructed directly if some unknowns v2 = (vm+1,..., vn) given arbitrary and determine the rest v1 = (v1,..., vm) using

fylV1 = -(fy2V2 + ft),

where, fyi = (fik), fy2 = (fi i), i,k = 1,...,m, I = m + 1,...,n. Then the system in equation (1) will have the form [2]:

y1 = -f-1(fy2V2 + ft), y2 = v2(y,t), (4)

To use the above system (4), it is important to require the validity of initial condition (2) and the relation det(/y1) = 0. In some cases, this condition is satisfied. Unfortunately, the requirement det(/y1) = 0 is not the only requirement to use the method mentioned above for the construction of differential equations. In the application of this method one should have in mind the possibility of integration error accumulation for equation (4), which in the course of time, leads to the destruction of the constraint equation (2). Consequently, equation (1) should be constructed so as to ensure that the solution y = y(t) satisfying the initial conditions y(t0) = y° with

||/(yVo)|| < e (5)

deviates from constraint equation (2) in the course of numerical integration by a quantity of the order of e:

11/(^)ll < (6)

To define the right-hand side v(y,t) of equation (1) instead of equation (3), one uses the relation [2]

fy v + ft = a, (7)

where a = a(f,y,t), a(0,y,t) = 0. It follows from (7) that the vector v should be determined as a solution of m linear equations with n unknowns:

Av = b, (8)

where A = fy, b = a — /t.The structure of the general solution of equation (8) is described by the following theorem [2]:

Theorem 1. The set of all solutions of the linear systems (6), in which the matrix A has rank r is determined by the relation

v = c [AC] + A+b = cvT + vu, (9)

where c is an arbitrary scalar quantity, [AC] = [A1,..., AmCm+1,..., Cn-\] is vector product,Ai = (Aij), and arbitrary CT = (CTj), t = m + 1,... ,n — 1, j = 1,... ,n, A+ = A? (AAT )-1 and the component [AC] • of a determinant is

5j1 6j2 . . .

An A12 ... A1n

[AC], = Am1 Am2 ... A , where

Cm+1,1 Cm+1,2 . . Cm+1, n

-Cn-1,1 Cn-1,2 . . C n — 1 ,n -

(i

0, if j = k; 1, if j = k.

Example 1. A slider crank mechanism is given in the figure below.

Figure 1. Slider crank mechanism

The translational displacement x of a unit mass and the angular displacement d of crank are joined by the rod I. Then the constraint equation will be

f (y) = r sin d1 — I sin d2 = 0

and Av = b, where A = fy, b = a = kf, as ft = 0, the solution of this can be put as

v = cvT + Vv, vT = [fy] , V" = f+a

where

v\ = —I cos d2, v2 = — r cos 0\,

kr cos 6\(r sin 6\ — I sin 62) v —klcos d2(r sin 6\ — I sin 62)

V2 =

1 (r cos (9i)2 + (l cos d2)2 ' where k > 0. Then

V1 = cv\ + Vi = —cl cos d2 +

(r cos 0\)2 + (I cos d2)2

kr cos 0\(r sin 6\ — I sin d2) (r cos 0\)2 + (I cos d2)2 '

_ kl cos d2 (r sin 6\ — I sin d2)

V2 = CV2 + v2 = —cr cos 01------—— .

(r cos 0\)2 + (/ cos 02)2

Taking c =1, r = 1,1 = 2 and the value of k to be choosen depending on the stability of the solution.

_ k cos 01 (sin d1 — 2sin d2) 2k cos 02(sin d1 — 2sin d2)

v1 = —2cos02 +—--^ i -^ , v2 = — cos01 — -

(cos Ô1)2 + (2 cos Ô2)2

(cos 01)2 + (2 cos 02)2

3. An Approximate Solution of Differential Equations

Many systems involving differential equations are so complex, or the systems that they describe are so large, that a purely mathematical analysis is not possible. It is in these complex systems where computer simulations and numerical approximations are useful [3].

Most often, systems of differential equations can not be solved analytically. Algorithms based on numerical methods are therefore needed. By numerical integration, we mean to compute, from y0 (the initial condition), each successive result y1, y2, y3, ... that satisfy equation

§ = «y,«. d0)

An algorithm is thus a program that computes as precisely as possible yn+1 from yn,where, yn = y(tn), tn+i = tn + r, and y = c [fyC] + f+(Kf - ft).

Of course, y can be a vector and the equations can be non-linear differential equations [3].

3.1. Euler's method

Let us assume that the initial conditions t0 and y(t0) be known for the solution of the equation (10). Setting a = Kf and using the right side of (10) construct the equation

yn+1 =yn + ryn. (11)

For a differential equation the following assertion holds [2]:

Theorem 2. There are constants a, r1te and matrix K(y, t), which if

1) || f(y0, io)|| < £,

2) t < T]_, and for all y = yn,t = tn,n = 0,1,... the inequality is fulfilled

3) ||J + tK(y, i)|| < a < 1,

4) T2 ||f(2)H < (1 — a)£,

where f(2) = vTfyTyv + 2fyt + ftt, then the solution of the difference equation (11) will satisfy the condition || f(yn, in)|| < e, for any n = 1, 2,....

The above theorem helps us to estimate the value of the constant matrix K in our example1, so that the solution is stable.

Example 2. Consider the problem what we have constructed in Ex. 1, that is,

01 = V1(61, 62), O2 = V2(O1, 02),

kr cos d1 (r sin d1 — I sin d2)

V\ = cvl + v" = — cI cos d2 +

(r cos di)2 + ( I cos 02)2

T , V n COS ^(rsin^i- ¿sin 02)

V2 = cv2 + v2 = —cr COS 61--------— (12)

(r cos 0\)2 + ( I cos 02)a

System (12) has a particular solution

f(01, #2) = r sind 1- ¿sin62 = 0 (13)

dand be calculated as 6= 6?+rv?, = 92^+tv2{, 60 = 0, and = 0. From (13) it follows that condition 1. holds clearly: |/(0, 0)| = 0 < e. Condition 3 imposes a limit and k: |1 + rfc| < a < 1, that is,

1 — a , a + 1 iN - < k <--(14)

To check condition 4 for this problem, consider the right hand side of (12) and rewrite it as |r sin Q1 — ¿sin 02| < e, (r cos 01)2 + ( I cos 62)2 > I2 cos2 02 > I2 — r2 sin2 Q1 — 2r |sin 01| — s2 > I2 — r2 — 2r£1 — e 1 = I2 — ( r + £1 )2. Then it follows:

| V11 < cl +

kr£

I2 — (r + £1)2

= W1, |V2| <cr +

kl£

So that f(2) = 0v22 + 0

I2 — (r + ei)2

= W2

f(2)

d2f 2, 92f

de 1

12 +

de 2

<

dH de2

+

du

de 2 d 2f

d 12

12 +

d2

d 22

< rW1 + IW2 = F.

Conditions 2 and 4 take the form:r ^ T1, t2F ^ 2(1 — a)e. Taking r = 1,1 = 2, c = 1 in equation (12) and e = £1 = 10—4,setting (approximate values):W1 = 3, W2 = 2.Then F = 17, correspondingly condition 4 takes the form 17 x 104rf ^ 2(1 — a) and again let us assume n = 10—3, a < 0.9. Taking a = 0.9, we obtain a restriction for k: 100 €k < 1900.

2

2

2

2

1

2

2

2

3.2. Runge-Kutta Method

Runge-Kutta methods introduce values between tn and tn+1, and evaluate v at these intermediate points [4]. The general Runge-Kutta method is defined by

yn+1 =yn + rh(tn,yn, t), (15)

where h(tn,yn, t) = ^2^ crkr, with

^¿n + rar, yn + r ^ brsks^j

k1 = v(t n,yn), kr = v [tn + rar ,yn + r }^brsks) , (16)

and ar = Y1 r—1 brs, for r = 2, 3,... ,R. These constants, cr, ar and brs, are determined to ensure the highest order accuracy for the method. To establish the order of accuracy, consider the Taylor series expansion of y(tn+1)

yn+1 = yn + ryn + yyn + ... = yn + r (yn + 2yn + ...)

From this we may recall that the Taylor's series method of order p for h in (15) can be written with

p—1

h(tn, yn, r) = hT(tn, yn, t) = v(tn, yn) + 2v(1)(tn, yn) + ...+v(p—1)(tn, yn), (17)

where v(l) = jrv(tn, yn), i = 1, 2,..., (p — 1). We further expand k2, k3 and k4 upto order 3 and represnt the summary as follows:

k2 = v( tn + T(l2,yn + T(l2k1), T2 T3

k2 = vn + Ta2Fn + —a2Gn + —as2Hn + r4Rkv4, (18)

where vn = v(tn, yn), Fn = f™ + vnfy and

Gn = fn + 2vnfn + (vn)2 , Hn = (vn)3 f^yy + 3(vn)2 f^yt + 3vnfntt + Ku, h = v(tn + таз, yn + r(b3iki + Ьз2к2)), аз = 631 + 632,

k3 = vn + ra3Fn + у (a2Gn + 2a2&32 fnFn)+

3

+ — (a3#n + 6a2b32a3lnFn + 3a3&32 fnGn) + т4Д£4, (19) 6 y

where In = vnfZy + fynv

k4 = v(i n + ™4 ,yn + r( b4lkl + 642 k2 + b43k3)), a4 = 641 + b42 + &43,

2 3

k4 = vn + Td4Fn + —a4Gn + T2(a2&42 + a3&43) fnFn + у(a2&42 + a3&43) /yGn+

2

+ T2a2b32( fy)2Fn + T3(a2b42a4 + a3&43a4)/nFn + у a2#n + T4#J4 . (20) using equation (15) and Regrouping similar terms we get h(tn, yn, t) = (ci + C2 + C3 + C4)wn + t(C2a2 + C3a3 + C4a4)Fn+

T2

+ у (a2c2 + а2сз + a4c4)G" + r2( сза2Ъз2 + C4 (a2b42 + a3b32))( fyF )n+

t3

+ -TT(( C20-2 + C3a3 + C4al)Hn + 2( C3a2b32tt3 + 0402642^4 + C4a,3b43a,4)(F I)n+ 6

+ 3( C3a2b32 + C4a2&42 + C4a3&43)( fyG)n + 6 C402b32( f?)2Fn) + T4^4. (21) Now we have to match equation (17) with (21) to find unkown parameters.

3.2.1. When R = 2

As c3 = 0, equation (21) reduces to

T2

h(tn, yn, t) = (C! + C2)vn + Ta2C2Fn + y C2a2Gn + r3Rk3, (22)

and comparing this with (17), we have

Cl + C2 = 1; a2C2 = 2. (23)

This gives a set of two equations in three unknowns and there exists a one parameter family of solutions.

Combining (22) with (15), we can rewrite A yn as

yn+1 -yn = Th(tn, yn, t) = t(d + C2)vn + T2a2C2Fn + T3Rk3

2

n n 3 k3

A yn = TVn + y Fn + r3Rk3. (24)

Consider the row decomposition of vector fn+1 = f(yn+1, tn+1) [2]: fn+1 = fn + fn^yn + Tfn + 2((Ayn)T4%Ayn + 2rfnyAyn + ^№ + r3Rk/, (25) where R}3 = (R^ ,---,Rk/m), R1}^ = (EP,g,r fTpqr Ay^y^

+ T Hp,qfnpqtA 2/pA Vq + 3r X] p fi,pt tA Vp + T fi^tt t), and /i^pqr, fi,pqt, fi,pt t, fi^ttt

are the values of third partial derivative, to be determined for y = yn + BAyn,t = tn + 6t. The matrix B and the scalar d admit the corresponding intermediate values of derivatives.

Theorem 3. If a solution of (25) was obtained using a difference scheme of the second order of accuracy (24) and for all values of the variables y = yn,t = tn, n = 0,1,...,N, the values T,q > 0,c2,a2, the vector f°, the matrix K(y,t) and the remainder in Taylor series expansion T3Rk3 satisfy the inequalities: || f0\\ < e, 2a2c2 = 1, t3 \\Rk31| < (1 - q)e, I + tKn + ^(Kn + (Kn)2) < q< 1, then \\fn\\ for alln = 1,2,...,N.

Proof. Assume that the inequality || f°\\ ^ e holds for the initial condition y (t0) = y0. The equality A yn = r(1 — c2)vn + tc2v(tn + ra2, yn + ra2vn), which is obtained from (24) can be written in the form

A yn = r(vn + t C2a2Fn) + r3R},

(26)

where R}3 = (R}3,...,R}3m), R} = ^(Ep.q ZiVW +2Y,pv^tptv^ + v^t). After substituting (26) in equation (25) and rearranging terms we get

n+1 _ rn

fn + rFn + r2C2a2( fyF )n + -(v+fyyv + 2 ftyV + ftt )n + r3Rk/, (27)

further using the equalities f = K(y, t)f = fyV + ft = F, f = Kf + K2f = fyF + G, equation (27) can be represented as

2

fn+1 = (I + rKn + T2C2a2(Kn + (Kn)2))fn + ^-(1 — 2a2C2)Gn + t3R}3 ,

2 J

But from equation (23): a2c2 = 1/2,

2

/n+1 = (I + rKn + — (K + (Kn)2 ))/n + T3Rkf3, 2 ^

Now, taking the hypothesis of the theorem in to account one gets

n+1 \

<

I + TK + y(Kn + (Kn)2)

\\/n\\ + r3

R k

< qe + (1 — q)e = e.

2

2

3.2.2. When R=3

We can match (21) with (17) and the following set of equations are satisfied:

C1 + C2 + C3 = 1, a2C2 + a3C3 = 1/2, a^C2 + a?3C3 = 1/3, a2b32C3 = 1/6. (28)

These are four equations in six unknowns and there exists a two-parameter family of solutions. Combining (21) with (17), and substitute k1,k2 and expanding by Taylor

we can rewrite A n as

A yn = 7~( Cl + C2 + C3)vn + T2(a2C2 + a3C?)Fn + T3a2b32C3Fn +

-3

+ y(a2c2 + a2C3)Gn + r4Rk4 (29)

Equation (29) can further be simplified as

A yn = tMivn + t2M2F? + t3M4F2 + ym?G? + T4Rk4,

(30)

where Mi = ci + C2 + C3, M2 = a2C2 + a3C3, M3 = C2a2 + C3a2, M4 = a2CSbS2. Now, recall the decomposition of vector fn+i = f(yn+i, tn+i);

fn+i = r + f:Ayn + rfn + ^((Ayny fnTyAyn + 2^Ayn + r2/£)+

1

n\T rn A „.n 1 o>_ i*n A „.n 1 2 rnx

1

/y—y 1 ' Jt T- y > JyTy y ^ Jty^y "T / Jit;

fn \n,n \n,n \n,n ^ 1rfn \n,n \n,n ^ 1 r2fn rr-3 fn \ _1__4

+ -( fnyy AynAynAyn + 3rfnyt AynAyn + 3r2 /yntt Ayn + r3 f^t) + r4R^. (31)

Theorem 4. If a solution of (31) used a difference scheme of the third order of accuracy (30) and for all values of the variables y = yn,t = tn, n = 0,1,... ,N, the values t,5 > 0, f0, the matrix K(y, t) and the remainder in Taylor series expansion T4Rki satisfy the inequalities:

||f0H < e, t4 ||Rk41| < (1 — S)e,

2 T3

i + rKn + -—(Kn + (Kn)2) + ^-(Kn + 3KnKn + (Kn)3) < 1, 2 6

then llfnll < e, for all n = 1, 2,...,N.

Proof. Substituting (30) in (31) and further simplifying, one gets

T2 T3

n+i = fn + tFn + t2M2 fy!Fn + —Gn + —M3 fnGn+

2

2

+ t3M4( CYF2 + T3M2(FI)n + —Hn + T4Rk4, (32) y 6

where Hn = (vn)3f^ + 3(vn)2f njt + 3vnf^ + /ttt. And recall the equalities

f = K(y, t) f, f = Kf + K2 f, '/ = (K + 3Kk + K3) f, f = fyv + ft = F, f = fyF + G, f= flF + fyG + 3FI + H.

fn+i = + tKn + y (kn + (Kn)2) + y (Kn + 3KnKn + (Kn)3^ fn+

+ T4Rk4. (33)

Taking into account the given conditions and assuming that K is constant [2] the conclusion follows

n+i \

<

T2 T3

i + tkn + ^-(Kn)2 + —(Kn)3 26

\\fn\\ + t4 \\Rk4 II <fe + (1 - S)e = e.

3

3.2.3. When R = 4

Matching equation (21) with (17) we get the following equations:

C1 + C2 + C3 + C4 = 1,

C2<i2 + C3d3 + C4tt4 = 1/2, a2,C2 + a?3C3 + a4c4 = 1/3, c^aO^ + C3a33 + 04a33 = 1/4,

0312^32 + C4a2b42 + C4a3b43 = 1/6, C3a2b32a3 + 04(12^)42(14 + C4a3b43a4 = 1/8,

C3a2>b32 + C4a2&42 + C4a2&43 = 1/12, C4a2b32 = 1/24.

But A yn = rh(yn, tn, t) = t(c 1k1 + c2fc2 + c3k3 + c4k4). After replacing corresponding expressions for k1,k2 ,k3, k4 and simplifying the result we get:

2 3 4

Ayn = TVn+2pn(Gn+/yn^n)+24(^n+(FJ)n+fyn°n+(/yn)2^n)+T5Rk5 (34)

and recalling the decomposion of fn+1 = f(yn+1, tn+1);

r+1=r+/ynAyn+r/n+2f(n2)+3!/(n3)+4!f(n4)+r5Rks, (35)

where

/( n2) = ((A yn)TryTy Ayn + 2r fty Ayn + r2 fli), f( n3) = (/ynyy AynAynAyn + 3rfZyt AynAyn + 3r2 Ayn + r3 f^t),

/( n4) = f^y AynAynAynAyn + 4r /ynyyt AynAyn Ayn + 6r2 /ynyt tAynAyn+

+ 4 t 3 f^t tA yn + r4 f^t

and Rf is obtained from fifth partial derivatives, to be determined for y = yn + BAyn,t = tn + Ot. The matrix B and the scalar d admit the corresponding intermediate values of derivatives.

Theorem 5. If a solution of (35) used a difference scheme of the fourth order of accuracy (34) and for all values of the variables y = yn, t = tn,n = 0,1, 2,..., N, the values t,5 > 0, f0, the matrix K(y, t) and the remainder in Taylor series expansion T5Rk5 satisfy the inequalities:

||/°|| < e, T5 11Rks || < (1 — S)e,

T2 T3 T4

I + rKn + -(Kn)2 + — (Kn)3 + ~.(Kn)4 2 6 24

< 1,

then || fn\\ for all n = 1,2,...,N.

Proof. Substitute (34) in (35) and rearranging terms gives:

T2 T3

fn+1 = fn + rFn + —(f^n'Fn + Gn) + —((f^)2Fn + f£Gn + 3FnIn + Hn)+ 2 6

4

+ 2^(( fy)3Fn + (4n)2Gn + f^Hn + 4GnIn + 6FnNn + 7f;FnIn + 3 f^ (Fn)2 + Mn)+

+ r5Rk5, (36)

where N = v2fyyy+2vfyyt+fytt and M = v4fyyyy+4v3fyyyt+6v2fyytt+4vfyttt+ftttt.

90 Bulletin of PFUR. Series Mathematics. Information Sciences. Physics. No 3, 2013. Pp. 81-91 Let us consider the relation f = K(y, t)f and respective derivatives

f = K(y, t)f = fyv + ft = F, / = (K + K2)/ = fyF + G,

f = (K + 3KK + K 3)f = f2F + fyG + 3FI + H,

"f = (K + 4KK + 3K2 + 6KK2 + K 4)f =

= f3F + f?G + /y (H + 7FI ) + 3 /yyF+4G/ + 6F^ + M

In relation with (36),

T2 ■ T3

jn+1 = jn + TKnjn + _(Kn + (Kn)2)fn + _(Kn + 3^nKn + (f")3)/"+

2 6

4

+ 2^(Kn + 4KnKn + 3(K)n2 + 6KKn(Kn)2 + (Kn)4 )fn + r5^fc5. (37)

Taking into account the given conditions and the fact that K as constant, we conclude

rn+1 \

<

T2 T3 —4

n , f Tsn\2 i f Tsn\3 I / T;"n\4

I + rKn + -(Kn)2 + -(Kn)3 + ^(Kn) 2 6 24

11/1 +

+ T° \\Rk51 < fe + (1 - S)e = e.

Example 3. Consider the problem, which we have constructed in example 1 and take = 1, = 1, = 2,

_ k cos 01(sin d1 — 2sin d2) 2k cos 6>2(sin d1 — 2sin d2)

v1 = —2cos02 +—--^ i -^ , v2 = — cost)1 — -

(cos di)2 + (2 cos Ô2)2

(cos di)2 + (2 cos é»2)2

We can easily see that the solution using MATLAB is stable for large values of the constant k as shown below in the figure.

Figure 2. Solution using MATLAB

References

1. Arabyan A., Wu F. An Improved Formulation for Constrained Mechanical Systems // Kluwer Academic Pub. — 1998. — Pp. 49-69.

2. Muharlyamov R. Equations of Motion of Mechanical Systems. — PFUR, 2001. — (in russian).

3. Gonze D. Numerical Methods to Solve Ordinary Differential Equations // P. J. Brief Bioinform. — 2009. — Pp. 53-64.

4. Lakoba T. Runge-Kutta Methods. — University of Vermont, 2006. — Pp. 15-20.

УДК 531.3

Решение дифференциальных уравнений движения для механических систем со связями

Р. Г. Мухарлямов*, А. В. Бешау^

* Кафедра теоретической механики Российский университет дружбы народов ул. Миклухо-Маклая д. 6, Москва, 117198, Россия

^ Кафедра математики Бахрдарский университет Эфиопия, Бахрдар

В работе рассматривается задача построения систем дифференциальных уравнений по известным частным интегралам. Приводится метод определения правых частей систем дифференциальных уравнений, основанный на определении общего решения системы линейных алгебраических уравнений с прямоугольной матрицей коэффициентов. Предлагается использовать для численного решения построенной системы дифференциальных уравнений метод Рунге-Кутта. Для рассматриваемой задачи ранее были использованы простейшие разностные схемы первого порядка и метод Рунге Кутта для случая линейных дифференциальных уравнений возмущений связей с постоянными коэффициентами. В статье получены ограничения на коэффициенты уравнений возмущений связей, зависящие от фазовых координат системы, при решении дифференциальных уравнений методом Рунге-Кутта. Подробно рассмотрены случаи разностных уравнений первого порядка, состоящих из нескольких стадий. Получена общая форма условий стабилизации уравнений связей. Метод иллюстрируется на примере решения кинематической задачи кривошипно-шатунного механизма.

Ключевые слова: дифференциальные уравнения, численное интегрирование, кинематические ограничения, стабилизация, ряд Тейлора.