CC BY
36
UDC 531.3
On Solving Differential Kinematic Equations for Constrained
Mechanical Systems
A. W. Beshaw
Department of Mathematics Bahr Dar University Ethiopia, Bahr Dar
This paper proposes a method for constructing the kinematic equations of the mechanical system, which imposed geometric constraints. The method is based on the consideration of kinematic constraints as particular integrals of the required system of differential equations. Runge-Kutta method is used for the numerical solution of nonlinear differential equations. The developed methods allow us to estimate the range of variation of the parameters during the numerical solution which determine conditions for stabilization with respect to constraint equations. The numerical results illustrate the dependence on the stabilization of the numerical solution is not only due to the asymptotic stability with respect to the constraint equations, but also through the use of difference schemes of higher order accuracy. To estimate the accuracy of performance of the constraint equations additional parameters are introduced that describe the change in purpose-built perturbation equations. It is shown that unstable solution, with respect to constraint equations, obtained by the Euler method can be stable by using Runge-Kutta method.
Key words and phrases: kinematic constraints, pseudo-inverse, approximate solution, stability, Euler's method, Runge-Kutta methods.
1. Introduction
Let the kinematic state of a mechanical system be determined by the position coordinates qi, q2,..., qn and dependent velocities qi, q2, ■ ■ ■ , qn described by the system of ordinary differential equations [1]:
q = v(q, t), q= q<ERn, (1)
where v is a vector valued function and we shall assume that the function v together with its first partial derivatives are continuous in some domain D c Rn x R.
The functions q(t) satisfy the equations (1) and take prescribed values q° at the initial time t0:
q(to) = q°. (2)
And assume that the system is subjected to m holonomic constraints:
f(q, t) = 0, where, feRm, m (3)
then it is usually assumed that the initial values of the coordinates also exactly satisfy the constraint equations
/( Q0, to) = 0. (4)
As stated in [2] relations (1) should be constructed according to the constraint equations (3), which under relations (4), form a set of particular integrals of the system of differential equations (1).
Moreover, on all solutions q = q(t),t > t0, of system (1) satisfying condition (3), the vector ( , ) of the right-hand side should be chosen so as to satisfy the conditions
fqV + ft =0, (5)
Qfi Qfi
where fq = (fij), ft = (fit), fij = ^, fit = , i = ^2,...,m, j = l,2,...,n
which form a system of m linear equations with n unknowns v1,v2,... ,vn forming the vector v.
2. Solution of systems of differential equations
The equation (5); the system of m ordinary differential equations with n unknowns, can be rewritten as
Mv = b, (6)
where M = fq is an m x n matrix and b = —ft is a column matrix. As we can see from the dimension of matrix M, we cannot apply inverse (the usual method) to solve the system in (6).
The structure of the general solution of equation (6) is described by the following theorem [3]:
Theorem 1. The set of all solutions of the linear systems (6), in which the matrix M has rank m is determined by the relation
v = c[MC ] + M+b, (7)
where c is an arbitrary scalar quantity, [MC] = [M\ ... MmCm+\ ... Cn-\] is vector product, Mi = (Mij), and arbitrary CT = (CTj), t = m + l,... ,n — l, j = l,... ,n, M+ = MT (MMT )-1.
Note that M + = MT(MMT)-1 is pseudo-inverse (Moore-Penrose inverse) of matrix M. Given an m by n matrix M, the n by m matrix M + is a unique matrix that satisfies the following four relations [4]:
1. MM+M = M;
2. M+MM+ = M+;
3. (MM+)T = MM+; and, (8)
4. (M+M)T = M+M.
Since M+M and MM + are both symmetric matrices,
MT (MMT )-1 = MT (M+)T M + = (M+M )T M + = M+MM + = M+.
In our case replace fq instead of M in (7), accordingly fq = M, M + = (fq)+ and
V = c[fq C] + (fq )+(—ft). (9)
In order to compensate the deviation of (5) during numerical integration it is necessary to use a correcting term and restate equation (5) as [3]:
fq v + ft = F, (l0)
where F = F(f, q, t), F(0, q, t) = 0.
Definition 1. The function f (q,t) is a particular integral of (l) if the equality (l0) and F(f,q,t) vanishes when the condition f (q,t) = 0.
It is usually assumed that the equality (4) is satisfied exactly. However, if the initial conditions q(t0) = q° do not satisfy the constraints (3), depending on how you
will be determined by the function F(/; q, t) in equation (10), the solution will either approach the manifold (3), or away from it.
Having a general solution of equation (10), we can obtain the desired set of systems of differential equations that determine the kinematic relations between the coordinates of constrained mechanical systems and their derivatives [3]
Q = c[fqC] + (fq)+(F - ft). (11)
In the expression (11), [ fqC] = [ fiq f2q ... fmqCm+i ...Cn-i] is vector product Q fi
fq = (fij), fij = QT^ (i = 1,... ,m, j = 1,... ,n) is row vector matrix and Cm+i(q), Q Qj
Cm+2(q),... ,Cn-i(q) are arbitrary vectors, (fq)+ = (fq)T(fqfq)T)-i. Equation (3) represents a collection of particular integrals of the system (11).
3. Numerical Solutions
The physical phenomena of the differential equations (11) are highly nonlinear and usually solved by numerical methods (numerical approximations). Let us assume that the initial value q(t0) = q° be given and suppose F can be expressed as a linear product of /; F = Kf, where K is m x m matrix. Then we may use one of different numerical methods to solve (11). Among these methods we use the Euler's method, second order, third order and fourth order Runge-Kutta methods.
3.1. Euler Method
For the system of differential equation in (11) we construct the equation:
q%+1 = q% + t q%,
(12)
where q% = q(tí), ti+i = U + t and (f = c[fqC](r,U) + (fq) + (Kf - ft)(q^ti), i = 0,1,2,....
In relation to (12) we state the following statement [5]:
(i) ||/°|| < e,
(ii) t < ^, and for all q = ql, t = U, i = 0,1,2,... fulfilled the inequality,
(iii) ||/ + tK|| < a < 1,
(13)
(iv) f(2)
<
2(1 - a)e
where f(2) = vT fqTqv+2fqtv+ftt, then the solution of the difference equation (12) will satisfy the condition || f( q1, ij)|| < e, for any i = 1, 2, ■■■, where a,e, n are constants.
i
3.2. Second Order Runge-Kutta Method
In Euler method there was no attempt made to include information from the point ti+\ at which the solution is being obtained. Heun (Second order Runge-Kutta) method uses Euler method but takes the result qi + t<f to calculate the gradient there and then, whereas Euler method just uses qi for the gradient over the interval, Heun method uses the mean of that value and the approximately calculated value at ii+i.
qi+1 =qi + 2>( U) + v(qi+ r, qi+Tv(qi, U))], (14)
where q1 = v(qz,ti). For our system (11):
v(q\ti) = c[fq C]{tM) + (fq )+(Kf - ft)(q,A). (15)
3.3. Third Order Runge-Kutta Method
The solution of (ll) can also be solved using Runge-Kutta method in more efficient order:
ql+1 = ql + r (ciki + C2 k2 + c3k3), (l6)
where the multipliers kj(j = l, 2,3) and Cj(j = l, 2, 3) are obtained as in [6]:
qi+1 = ql + \ (ki +4k2 + k3), (l7)
6
and
k1 = v(qi,ti); k2 = v (qi + 2h,U + 2); = — rh +2rk2,U + r).
Now to solve (ll) in this method, expand each expressions of k1,k2,k3 up to order three:
v
(q\U) = c[fqC](9Vi) + (fq)+(Kf - ft){qi>U)
k2 = v (,f + 2v,ti + 2) ;
fa = v (qi - 2rv(qi,ti) +2tv (qi + 2v(qi,ti),ti + 2) ,ii + r) .
The second and third order Runge-Kutta methods are available in MATLAB as built in functions called ODE solvers ode23. It is also possible to construct user defined functions and solve ODE problems with efficiency as desired. The approximate solution will be close to the actual solution if n is taken sufficiently large enough.
3.4. Fourth Order Runge-Kutta Method
The fourth order Runge-Kutta method is the popular version of the members, for which reference can be made to any book on numerical methods. They require the ability to solve the initial value problem for a system of ordinary differential equations at a number of points intermediate between ti and ti+1 [7]. If the function v and the data q°,t0 are given, then
qi+1 = qi + 2 (k1 +2k2 + 2k3 + k4), (l8)
6
where
k1 = v(qi, U) = c[fqC](q,,u) + (fq)+(Kf — ft)(q,,u);
k2 = v (ti + 2, q* + ^v(qi, t^ ; h = v (ti + 2, qi + 2V + 2, qi + 2v(ti, qi))) ; fa = v (ti + t, qi + tv (ti + 2, qi + (ti + 2, q* + 2v(q\ U))) ) .
Here ql+1 is the Runge-Kutta approximation of q(ti+1), and the next value (ql+1) is determined by the present value (qz) plus the weighted average of corresponding (2,3,4) increments, where each increment is the product of the size of the interval,
t, and an estimated slope specified by the function G on the right-hand side of the differential equation [8]:
- k1 is the increment based on the slope at the beginning of the interval, using q ;
- k2 is the increment based on the slope at the mid-point of the interval, using q+ 2^1;
- fa is the increment based on the slope at the mid-point of the interval, but now using q + 2 fa;
- fa is the increment based on the slope at the end of the interval, using q + rfa. Note that the 4th, 5th order Runge-Kutta methods are combined together and
termed as 'ode45', which is available in MATLAB. It is also possible to build user defined functions.
Example 1. Let us consider systems of ordinary differential equations of two variables:
{X = v1(x, y),
1(( , ^ (l9)
y = V2(x, У),
and satisfy a constraint equation:
f(x, y) = l/2(9x2 + y2 — 9) = 0. (20)
Differentiating (20) with respect to time and assigning the result to a function F we
get
9xX + yy = F, F = F(f ,x, y). (2l)
We assume this function F as a linear combination of the given constraint function f; F = kf.
The solution for this equation (2l) can be found using Theorem l, it gives
9 x
x = °y+ 2 i 2 (kf),
81x2 + M (22)
y = —9 cx + 81x2+7( kf).
Substitute in relation (22) we obtain
x = Cy + 2(8llk2X+ y2)(9x2 + V2 — 9), y = —9 cx +2(81xfc2 + ,2)(9x2 + ^ — 9).
Next we solve (22) numerically and investigate the stability of the solution as k varies. For simplicity we consider k as constant and c = 1.
case I: Let k > 0, say k = 2, and initial values (x0, y°) = (1,0). Therefore, relation (22) becomes
9 x
x = <9x2+y'—9»- ,23)
i = —9x + <9x2 + ^ — 9>.
As we can see from Fig. l, the numerical solution of (23) is unstable, neither it is solved in Euler method nor in Runge-Kutta methods. case II: let k = 0, the system reduces to
fx = y,
• 9 (24)
= —9 x.
a Position versus time k = 2; unstable solution
b Phase portrait graph for k = 2; unstable solution
Figure 1. The numerical solution of (23)
This system has exact solution: u(t) = (c1 cos 3t, c2 sin 3i) and take x0 = 1, y0 = 0 we get u(t) = (cos3t, —3sin3t).
case III: Let k < 0, say k = —2, the system (22) becomes
9x
x = y
81x2 + y2
y = —9x —
81x2 + y2
(9x2 + y2 — 9), (9x2 + y2 — 9).
(25)
The numerical solution of (25) is stable in any of the methods discussed above, see Fig. 2. But the stability of the solution depends on not only the value of the parameter k but also the method we use to solve it. For instance, taking k = —400 the solution is unstable when we use Euler method but it is stable in Runge-Kutta methods as shown in Fig. 3.
a Position-time graph for k = -2; stable solution
b Phase portrait graph for k = -2; stable solution
Figure 2. The numerical solution of (25)
From the above cases the solutions are
- Unstable for all values of k in case 1 (k > 0);
- Stable for case 2 ( k = 0) or it gives more sense if we say F = 0;
- Stable for case 3 ( k < 0) for certain domain of k values.
a Position versus time graph using Euler method for k = — 400; unstable solution
b Position versus time graph using Runge-Kutta methods for k = — 400; stable solution
Figure 3. The numerical solution of (25)
To determine the possible values of k so that the solutions obtained are stable we proceed as follows. To approximate equation (22) numerically we guess the possible values of the parameter k
Zi+1 =Zi + TZi, zi = z(ti), ti+i = ti + T, z=(x, y),
9 kx , 2 2 ^ ky , 2 2
x = °y+ ofai 2 ,—2\(9x +y — 9), y= —9cx + Ofai 2 X—2\(9x +y — 9). 2(81x2 + y2) 2(81x2 + y2)
Referring these with the statements of (22) [3]:
f(2) = Sx2 + g y2 = 9x2 + y2.
dx2 dy2
From the last condition of (13), we have |9x2 + y2\ <
2(1 — a)e
Letting k = c and using (22), condition 3 can be rewritten as
|2
2 + 2) + 90A Ixllyl |9x2 + y2 - 91 + X2(729x2 + y2) [Q^2 + y2 - 9|2 < 2(1 - a)e 9( +V ) + 81x2 + y2 + 4(81x2 + y2)2 " c2t2 .
Now we take |x| < 1.1, |y| < 3.1, 81x2 + y2 > x2 + y2 > 0.81, and substitute these to variables in the above inequality results
81.9 + 37.9 A + 3.4 A2 <
. 2(1 — a)e
22
2
'1
To have possible range of the discriminates of the inequality must be nonnegative:
^ — ^81.9 — ^^U 0, this leads, r^ .
V c2T1 ) ' ' V 48.3
Now choose c = 1, a = 0.99, e = 0.1, r < 0.0045, say r = 0.002. Then from the quadratic inequality we get -359.6 < A < 6.
1
On the other side from condition 2 above; |1 + rfc| < a < 1, we get
<k < —, -995 <k 0.5.
Hence —359.6 €k = X <-0.5.
4. Conclusion
The paper considers a method of constructing kinematic equation of mechanical system that contains holonomic constraints. It is shown that the sign and size of the factor (parameter) in the general equation affects the stability of its solution. Experimental investigations using MATLAB are made to find the range and sign of the factor especially for Euler numerical method.
A remarkable result is obtained that stability depends not only on the size and sign of factor but also on the numerical method. In the demonstrative example, the Euler method results in unstable solution for k = —400 whereas Runge-Kutta method produces stable solution, see Fig. 3. This paper can be used as an initial for further research on the factor (indicated in the paper) in other numerical methods.
References
1. R. G. Muharlyamova, N. V. Abramov, J. K. Kirgiyaev, Equations of Dynamical Systems with Program Constraints, URGU, 2013, in Russian.
2. R. G. Mukharlyamov, On the Construction of Differential Equations of Motion of Constrained Mechanical Systems, Differential Equations, Differential Equations 39 (3) (2003) 343-353.
3. R. G. Mukharlyamov, Equations of Motion of Mechanical Systems, PFUR, 2001, in Russian.
4. F. E. Edwadia, Recent Advances in Multi-body Dynamics and Nonlinear Control, ABCM 28 (3) (2006) 311-315.
5. R. G. Mukharlyamov, On the Equations of Kinematics and Dynamics of Constrained Mechanical Systems, Multibody System Dynamics (6) (2001) 17-28.
6. R. G. Mukharlyamov, A. W. Beshaw, Solving Differential Equations of Motion for Constrained Mechanical Systems, Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics" (3) (2013) 81-91.
7. J. D. Fenton, Numerical Methods, Karlsplatz (3) (2010) 19-22.
8. L. Xie, A Criterion for Hurwitz Polynomials and its Applications, Modern Education and Computer Science (1) (2011) 38-44.
УДК 531.3
Численное решение уравнений кинематики механических
систем
А. В. Бешау
Кафедра математики Бахрдарский университет Эфиопия, Бахрдар
Работа посвящена решению задачи стабилизации связей при численном решении дифференциальных уравнений, описывающих кинематические соотношения в механической системе. В статье предлагается метод построения системы дифференциальных уравнений, соответствующих кинематическим соотношениям в механической системе, на которую наложены геометрические связи. Предлагаемый метод основан на преставлении
кинематических связей в качестве частных интегралов соответствующей системы дифференциальных уравнений. Для определения численного решения нелинейных дифференциальных уравнений используется метод Рунге-Кутта. Разработанный метод позволяет в процессе численного решения дифференциальных уравнений оценить границы изменения параметров управляющих воздействий, которые соответствуют условиям стабилизации решения по отношению к заданным уравнениям связей. Результаты вычислений показывают, что стабилизация численного решения зависит не только от асимптотической устойчивости по отношению к уравнениям связей, но также от точности используемой той или иной разностной схемы. Для оценки точности выполнения уравнений связей вследствие стабилизации связей вводятся дополнительные параметры, изменение которых определяется специально построенными дифференциальными уравнениями возмущений связей. Показано, что численное решение, полученное методом Эйлера, которое оказывается неустойчивым, может оказаться устойчивым при использовании метода Рунге-Кутта.
Ключевые слова: кинематические ограничения, приближенное решение, псевдообратный, устойчивость, метод Эйлера, методы Рунге—Кутты.
Литература
1. Мухарлямов Р. Г., Абрамов Н. В., Киргираев Ж. К. Уравнение днамики системы с программными связями. — Юргу, 2013.
2. Mukharlyamov R. G. On the Construction of Differential Equations of Motion of Constrained Mechanical Systems, Differential Equations // Differential Equations. — 2003. — Vol. 39, No 3. — Pp. 343-353.
3. Мухарлямов Р. Г. Уравнение движения механических систем. — РУДН, 2001.
4. Edwadia F. E. Recent Advances in Multi-body Dynamics and Nonlinear Control // ABCM. — 2006. — Vol. 28, No 3. — Pp. 311-315.
5. Mukharlyamov R. G. On the Equations of Kinematics and Dynamics of Constrained Mechanical Systems // Multibody System Dynamics. — 2001. — No 6. — Pp. 1728.
6. Мухарлямов Р. Г., Бешау А. В. Решение дифференциальных уравнений движения для механических систем со связями // Вестник рУДН. Серия «Математика. Информатика. Физика». — 2013. — № 3. — С. 81-91.
7. Fenton J. D. Numerical Methods // Karlsplatz. — 2010. — No 3. — Pp. 19-22.
8. Xie L. A Criterion for Hurwitz Polynomials and its Applications // Modern Education and Computer Science. — 2011. — No 1. — Pp. 38-44.
