Mass in gravitational Field

One of the simplest problem in classical physics is that of the motion of a small mass influenced by the gravitational field generated by another heavy mass. The latter approximately remains stationary under the assumption that it is very large. The exact solution of this one-body system for negative total energy predicts that the small mass will follow an elliptical orbit with ne of its foci located at the heavy mass location. The orbit remains in a plane whose direction is specified by the initial conditions.

This seemingly simple problem can become a challenge when solved numerically. This stems from the very singular nature of the gravitational potential which diverges as $1/R$ as the distance between the masses, $R$ , becomes small. When the small mass approaches close from the heavy mass, the gravitational field can cause extremely large curvature in its orbit. Said in another way, the orbits possessing high ellipticity will likely be hard to compute numerically, even when using a method as accurate as the Runge-Kutta 4th order scheme.

This system is conservative, namely the total energy $E = K + V$ is a constant of the motion. A warning of the poor performance of the ODE numerical integration will be the lack of constancy of what should be constants of the motion, like the total energy $E$ .

This has led to the development of new methods for ODE numerical integration which are specially derived in order to maintain the constants of the motion really constant. We will see some of these methods in the following sections.

First, let's illustrate the problem. The gravitational model above will be used to do so. Let the mass $m$ be in the gravitational potential of a large mass $M$

$$V(r) = \frac{G M m }{r}$$

The product $G M$ is a constant that sets the strength of the gravitational field. $r$ is the distance between the masses, $r=\sqrt{x^2+y^2}$ .

The mass $m$ follows trajectories in a 4 dimensional phase space (not 6), $x$ , $y$ , $v_x$ and $v_y$ , since the motion is in a plane (exact). Newton equations lead to

$$\frac{d x}{d t} v_x$$

$$\frac{d y}{d t} v_y$$

$$m \frac{d v_x}{d t} -\frac{ G M m }{x^2+y^2} \frac{x}{\sqrt{x^2+y^2}}$$

$$m \frac{d v_y}{d t} -\frac{ G M m }{x^2+y^2} \frac{y}{\sqrt{x^2+y^2}}$$

The program orbit.c implements the solution of the this problem using the Runge-Kutta 4th order ODE integration scheme. The mass $M$ is assumed at some arbitrary position $[XC,YC]=[1.0,1.5]$ ( $[0,0]$ should be used to simplify coding) as an example of how to code a real two body system. The mass $m$ initial position is taken at $x0 = XC+d$ and $y0 = YC$ where $d$ is arbitrary and its velocity as $v_x=v$ and $v_y = 0$ where $v$ is arbitrary.

• Exercise Run the code over some range of values for $v$ .

Larger values of $v$ produce small excentricity orbits for which RK4 is very accurate. Check the orbit in $x$ and $y$ for closure and the energy as a time series for constancy.

Smaller values of $v$ produce higher excentricity orbits for which RK4 ceases to be accurate. Check the orbit in $x$ and $y$ for closure and the energy as a time series for constancy.

Of course, smaller $dt$ would improve accuracy. But the point is that larger excentricity orbits would eventually lead to inaccuracies .

These concerns will be addressed in the next sections.