The Two-body problem

Various implementations of the Two-body problems are described

Available with Zymplectic v.0.1.0

The Two-body problem is about predicting the motion of two interacting objects. Here we will consider two heavenly bodies referred to as objects with masses m₀ and m₁ gravitationally bound by the potential V. The energy of the system is denoted by the constant Hamiltonian H equal to the sum of the kinetic energy T and potential energy V: $$H(\mathbf{q},\mathbf{p}) = T(\mathbf{p}) + V(\mathbf{q})$$ $$T(\mathbf{p}) = \frac{\mathbf{p}_{object 0}^2}{2m_0} + \frac{\mathbf{p}_{object 1}^2}{2m_1}$$ $$V(\mathbf{q}) = -G\frac{m_0m_1}{|\mathbf{r}|}$$ where q and p are vectors of canonical coordinates for position and momenta respectively, G is the gravitational constant and r = q₂ - q₁. Note that this Hamiltonian system is separable because T is only a function of p while V is only a function of q. Four implementations of the two-dimensionel two-body problem are described in the following section.

Two-body problem in Cartesian coordinates (system 1)

First obtain: $$\frac{\partial V}{\partial \mathbf{q}} = G\frac{m_0m_1\mathbf{r}}{|\mathbf{r}|^3}$$ $$\frac{\partial T}{\partial \mathbf{p}} = \frac{\mathbf{p}_{object0}}{m_0} + \frac{\mathbf{p}_{object1}}{m_1}$$ Given a set of appropriate constants, we consider the practical implementation of this system in the code below. The kinetic energy of the system can be split into two parts. The contribution of the first object 0 is introduced by the connector "Z.T(T,dT,0,0)" and the second object 1 by "Z.T(T,dT,1,1)", where T is the energy function and dT is the gradient function. The calculation of the potential energy requires object 0 and object 1 and therefore requires the connector Z.V(V,dV,0,1); The function T (ignore the argument) is a function returning the kinetic energy of either object whereas dT calculates the gradient with components d1 and d2. Note that obj1 refers to the first connected object in Z.T(T,dT,obj1,obj2). The V and dV functions utilize the coordinates of both objects with q and d referring to obj1 and Q and D referring to obj2.

//Copy-paste into Zymplectic or load as file to run the simulation
const double G = 1.0; //gravitational constant
const double m[] = {1.0,5.0}; //mass of object 0 and 1

#include "math.h"
double T(Arg) {
	return 0.5*(p1*p1 + p2*p2)/m[obj1];
}
void dT(Arg) {
	d1 = p1/m[obj1];
	d2 = p2/m[obj1];
}
double V(Arg) {
	double dx = q1 - Q1;
	double dy = q2 - Q2;
	return -G*m[obj1]*m[obj2]/sqrt(dx*dx + dy*dy);
}
void dV(Arg) {
	double dx = q1 - Q1;
	double dy = q2 - Q2;
	double d = sqrt(dx*dx + dy*dy);
	d = G*m[obj1]*m[obj2]/(d*d*d);
	d1 = dx*d; d2 = dy*d; //object 0
	D1 =-dx*d; D2 =-dy*d; //object 1
}

double _q_init[] = { //initial conditions
	5.0, 0.0,	//object 0
	0.0, 0.0};	//object 1
double _p_init[] = {
	0.0,-1.0,
	0.0, 1.0}; 

void main() {
	Z.V(V,dV,0,1);
	Z.T(T,dT,0,0);
	Z.T(T,dT,1,1);
}

Central-force problem in Cartesian coordinates (system 2)

It's possible to reduce the two-body problem to a single object problem by simulating the relative coordinate being the vector between the two objects: $$\mathbf{q} = \mathbf{q}_{object0} - \mathbf{q}_{object1}$$ where object 1 is now stationary. It can be shown that an identical simulation of system 1 is achieved by simulating object 0 with reduced mass mu or μ orbiting a stationary object body with a mass M. The stationary object is chosen to be located at origo and is intentionally excluded from the simulation.

//Copy-paste into Zymplectic or load as file to run the simulation
const double G = 1.0; //gravitational constant
const double m[] = {1.0,5.0}; //mass of object 0 and 1
const double mu = m[0]*m[1]/(m[0] + m[1]); //reduced mass
const double M = m[0] + m[1];
#include "math.h"
double T(Arg) {
	return 0.5*(p1*p1 + p2*p2)/mu;
}
void dT(Arg) {
	d1 = p1/mu; 
	d2 = p2/mu;
}
double V(Arg) {
	return -G*mu*M/sqrt(q1*q1 + q2*q2); 
}
void dV(Arg) {
	double d = sqrt(q1*q1 + q2*q2);
	d = G*mu*M/(d*d*d);
	d1 = q1*d;
	d2 = q2*d;
}

double _q_init[] = {5.0, 0.0}; //initial conditions
double _p_init[] = {0.0,-1.0}; 

void main() {
	Z.V(V,dV,0);
	Z.T(T,dT,0);
}

Central-force problem in polar coordinates (system 3)

It can be shown that the central-force problem (system 2) can be expressed in polar coordinates by the Hamiltonian: $$H = \frac{p_1^2 + \frac{p_2^2}{q_1^2}}{2\mu} - G\frac{\mu M}{q_1}$$ where q₀ is now the radius i.e. the length of vector from object 1 to object 0 and q₁ is the angle of this vector. It can be seen that the Hamiltonian is non-separable, implying that conventional explicit symplectic integrators are not applicable. A different algorithm is used which requires the gradient of H rather than the partial derivatives of T and V. $$\nabla H = \left(\begin{matrix}\frac{\partial H}{\partial q_1}\\\frac{\partial H}{\partial q_2}\\\frac{\partial H}{\partial p_1}\\\frac{\partial H}{\partial p_2}\end{matrix}\right) = \left(\begin{matrix}\frac{G q_1 \mu^2 (m_0 + m_1) - p_2^2}{q_1^3 \mu}\\ 0 \\ \frac{p_1}{\mu} \\ \frac{p_2}{q_1^2 \mu} \end{matrix}\right) $$ Since the gradient of the Hamiltonian H is utilized instead of T and V, the required connector is Z.H(H,dH,0) where H is the Hamiltonian function and dH defines the gradient with components named dq1, dq2, dp1 and dp2.

//Copy-paste into Zymplectic or load as file to run the simulation
const double G = 1.0; //gravitational constant
const double m0 = 1.0; //mass of object 1
const double m1 = 5.0; //mass of object 2
const double mu = m0*m1/(m0 + m1); //reduced mass

double H(Arg) {
	return 0.5*(p1*p1 + (p2*p2/(q1*q1)))/mu - G*mu*(m0 + m1)/q1; //for positive q1
}
void dH(Arg) {
	dq1 = (G*q1*mu*mu*(m0 + m1) - p2*p2)/(q1*q1*q1*mu);
	dq2 = 0;
	dp1 = p1/mu;
	dp2 = p2/(q1*q1*mu);
}

double _q_init[] = {2.0,0.0}; //radius,angle
double _p_init[] = {0.0,2.0}; 

double GLX[2],GLY[2];

#include "math.h"
void DrawFunc(){
	GLX[0] = 0;
	GLY[0] = 0;
	
	GLX[1] = obj_q[0][0]*cos(obj_q[0][1]);
	GLY[1] = obj_q[0][0]*sin(obj_q[0][1]);
	
	draw.points(GLX,GLY,2);
}

void main() {
	Z.OnDraw = DrawFunc; //replaces coordinate plots of polar coordinates with customized Cartesian coordinates
	Z.xmin = -2; //set visuel range
	Z.xmax = 4;
	Z.ymin = -3;
	Z.ymax = 3;
	
	Z.H(H,dH,0);
}

Adaptive central-force problem in Cartesian coordinates (system 4)

If the orbit is highly eccentric, one will notice significant Hamiltonian errors which may break the symplectic properties or contribute to a rotation of the Laplace-Runge-Lenz vector. To address this issue, one may pursue a method to perform more gradient evaluations near the periapsis. Unfortunately, explicit symplectic integrators are not compatible with adaptive time steps unlike other geometric integrators. Altering the time step size during integration results in significant and non-recoverable Hamiltonian errors. However, clever manipulation of the gradient can be utlized to implement the effect of adaptive time steps without changing the integration step size. Consider the following modification to the central-force problem in Cartesian coordinates (system 2): $$\nabla H = \left(\begin{matrix}\frac{\partial H}{\partial q_1}\\\frac{\partial H}{\partial q_2}\\\frac{\partial H}{\partial p_1}\\\frac{\partial H}{\partial p_2}\end{matrix}\right) = \left(\begin{matrix} q_1 \frac{G \mu M}{|\mathbf{r}|^3}\\ q_2 \frac{G \mu M}{|\mathbf{r}|^3} \\ \frac{p_1}{\mu} \\ \frac{p_2}{\mu} \end{matrix}\right) \cdot f(\mathbf{q})$$ where the scalar function f is the only distinction between this system and system 2. Note that f depends explicitly on the distance between the two objects. This implies that the gradient is no longer separable as T is a function of both q and p. An ideal choice of f would be a function proportional with the square of r in response to the inverse square law.

//Copy-paste into Zymplectic or load as file to run the simulation
#include "math.h"

const double G = 1.0;
const double m0 = 1.0;
const double m1 = 1.0;
const double mu = m0*m1/(m0 + m1); //reduced mass

double energy_H(Arg) {
	return -G*mu*(m0 + m1)/sqrt(q1*q1 + q2*q2) + 0.5*(p1*p1 + p2*p2)/mu;
}
void dH(Arg) {
	double d = sqrt(q1*q1 + q2*q2);	
	double f = q1*q1 + q2*q2; //f is the time scale factor
	d = G*mu*(m0 + m1)/(d*d*d);
	
	dq1 = q1*d*f; //time factor is 1 at apoapsis (based on _q_init)
	dq2 = q2*d*f; //time factor becomes small at close approaches
	dp1 = p1/mu*f;
	dp2 = p2/mu*f;
}

double _q_init[] = {1.0,0.0};
double _p_init[] = {0.0,0.2};

void main() { //this function is launched when executing Zymplectic.exe
	Z.H(energy_H,dH,0,0);
	Z.delay = 2.0;
	Z.track = 1;
}

Personal remarks

The two-body problem may not be a particularly interesting system to examine numerically as the results of this superintegrable problem are best presented in analytical contexts. However, The numerical treatment of the two-body problem as a non-separable Hamiltonian has some novelty as it seems unapproached by existing literature. In fact, the application of non-separable symplectic algorithm as mean of artificially controlling the step size for the two-body problem as well as other systems is as far as I know a completely unexplored aspect of explicit symplectic integrators.