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.