This paper addresses the description of equations of motion in Jacobi coordinates using osculating elements as an effective method for solving many-body problems. The equations of motion are formulated via osculating elements using Jacobi coordinates for the TRAPPIST-1 seven-planet exoplanetary system, characterized by variable masses and non-isotropic reactive forces. In many-body problems, the use of conventional heliocentric coordinates often complicates calculations. Jacobi coordinates consider the motion of each body relative to the center of mass of the preceding bodies, thereby allowing the problem to be decomposed into independent two-body problems. The bodies interact according to Newton's law of universal gravitation. Due to the variable mass, reactive forces are present, which further complicate the problem. The primary objective of this research is to establish a system of differential equations describing the evolution of osculating elements, treated as time-varying variables. While Jacobi coordinates simplify the mathematical framework, osculating variables provide a geometric representation of the motion (changes in orbital shape and orientation). Based on the equations of motion derived in Jacobi coordinates, a canonical perturbation theory has been developed, and equations for perturbed motion in analogues of Jacobi and Delaunay elements have been obtained. This approach enables high-precision modeling of the dynamics of multi-planetary systems. Based on the derived equations, the problem can be analyzed within a barycentric coordinate system using perturbation theory.