## 2 Models and Algorithms
### Model 1: Aperiodic Motion on Conic Section  
Adding to both sides of Equation (1) the same term \(\vec{F}\), we rewrite it in the form  
\[\frac{d^{2}\vec{r}}{dt^{2}}+Gm(t)\frac{\vec{r}}{r^{3}}+\frac{\gamma(t)}{2\gamma (t)}\frac{d\vec{r}}{dt}=\vec{F}. \tag{6}\]  
where  
\[\vec{F}=\frac{\gamma(t)}{2\gamma(t)}\frac{d\vec{r}}{dt}. \tag{7}\]  
Note that the perturbing force (7) coincides with (2) if \(m(t)=m_{0}(t)+m_{1}(t)\).  
Obviously, Equation (6) coincides with Equation (1), but it becomes integrable in the case of \(\vec{F}=0\). Indeed, one can readily check that at \(\vec{F}=0\), it has an exact solution that describes aperiodic motion of the body \(P_{1}\) on a conic section; the corresponding Cartesian coordinates are given by  
\[x =a\Big{(}(\cos E-e)(\cos\omega\cos\Omega-\sin\omega\sin\Omega\cos i )-\] \[\qquad\sqrt{1-e^{2}}\sin E(\sin\omega\cos\Omega+\cos\omega\sin \Omega\cos i)\Big{)},\] \[y =a\Big{(}(\cos E-e)(\cos\omega\sin\Omega+\sin\omega\cos\Omega\cos i )-\] \[\qquad\sqrt{1-e^{2}}\sin E(\sin\omega\sin\Omega-\cos\omega\cos \Omega\cos i)\Big{)},\] \[z =a\Big{(}(\cos E-e)\sin\omega+\sqrt{1-e^{2}}\sin E\cos\omega\Big{)} \sin i. \tag{8}\]  
The constants \(a,e,i,\Omega\), and \(\omega\) in (8) are the well-known Keplerian orbital elements, which are uniquely determined from the initial conditions of motion (see [3]). The eccentric anomaly \(E\) is determined by the well-known Kepler equation  
\[E-e\sin E=M, \tag{9}\]  
where the mean anomaly \(M\) is given by  
\[M=\frac{\sqrt{\kappa}}{a^{3/2}}(\Psi(t)-\Theta),\ \ \Psi(t)=\int_{0}^{t}\frac{dt}{ \sqrt{\gamma(t)}}\,\ \ \Theta=\Psi(\tau), \tag{10}\]  
where \(\kappa=G(m_{00}+m_{10})\). By \(\tau\) in (10), we denote the time when the body \(P_{1}\) passes through the pericenter.  
One can readily check that in the case of \(\gamma(t)=1\), solution (8) reduces to the corresponding solution to the two-body problem with constant masses, and the mean anomaly \(M\) becomes a linear function of time as \(\Psi(t)=t\) (see (10)). If the laws of mass variation \(m_{0}(t),m_{1}(t)\) are given, the function \(\gamma(t)\) defines the mean anomaly \(M(t)\) and Equation (9) enables one to find the eccentric anomaly \(E(t)\) as functions of time. Therefore, solution (8) defines the unperturbed motion of the body \(P_{1}\) on a conic section in terms of the time and 6 constants of integration \(a,e,i,\Omega,\omega\), and \(\Theta\), which are the Keplerian orbital elements (see [3]). However, the unperturbed motion of the body \(P_{1}\) becomes aperiodic because of nonlinear dependence of the mean anomaly \(M(t)\) on time.  
Using Equations (9) and (10), one can write the total time derivatives of the coordinates (8) in the form  
\[\dot{x} =-a\dot{E}\Big{(}\sin E(\cos\omega\cos\Omega-\sin\omega\sin\Omega \cos i)+\] \[\sqrt{1-e^{2}}\cos E(\sin\omega\cos\Omega+\cos\omega\sin\Omega\cos i )\Big{)},\] \[\dot{y} =-a\dot{E}\Big{(}\sin E(\cos\omega\sin\Omega+\sin\omega\cos\Omega \cos i)+\] \[\sqrt{1-e^{2}}\cos E(\sin\omega\sin\Omega-\cos\omega\cos\Omega \cos i)\Big{)},\] \[\dot{z} =a\dot{E}\Big{(}-\sin E\sin\omega+\sqrt{1-e^{2}}\cos E\cos\omega \Big{)}\sin i. \tag{11}\]  
where  
\[\dot{E}=\frac{\sqrt{\kappa}}{a^{3/2}(1-e\cos E)\sqrt{\gamma(t)}}. \tag{12}\]