Coupled Equations of Motion for Two-Body Gravitational Interaction

Welcome back. Last lesson, you scaled a fixed-centre gravitational model so that its equations could be interpreted without committing to kilometres, seconds, or a particular central mass. That model treated the central body as stationary. Now we remove that approximation: both bodies move, and each feels the gravity of the other.

This is the starting point for a genuine two-body orbital simulation. We will formulate the coupled equations in a common inertial frame, check their directions and units, and express them as the kind of initial-value problem you can eventually integrate in Python. The next lesson will show how a change of coordinates simplifies these coupled equations; for now, the aim is to understand precisely what is being simplified.

---

One separation, two forces

Consider two isolated Newtonian point masses $m_1>0$ and $m_2>0$. Their positions, measured from the same origin in a non-rotating inertial frame, are ${\mathbf{r}}_1(t)$ and ${\mathbf{r}}_2(t)$. Their masses are constant, and mutual gravity is the only force in the model.

Define the vector pointing from body 1 to body 2 as

$$\mathbf{d}={\mathbf{r}}_2-{\mathbf{r}}_1,d=\Vert \mathbf{d}\Vert .$$

The scalar $d$ is their separation. We require $d>0$: point-mass gravity is undefined when both positions coincide.

Newton’s gravitational law gives the magnitude of the force on either body:

$$F=\frac{G{m}_1{m}_2}{d^2}.$$

A magnitude is not yet an equation of motion. Gravity is attractive, so the force on body 1 points along $\mathbf{d}$, towards body 2. The unit vector in that direction is $\mathbf{d}/d$. Thus

$${\mathbf{F}}_1=\frac{G{m}_1{m}_2}{d^2}\frac{\mathbf{d}}{d}=G{m}_1{m}_2\frac{{\mathbf{r}}_2-{\mathbf{r}}_1}{\Vert {\mathbf{r}}_2-{\mathbf{r}}_1{\Vert }^3}.$$

The third power in the denominator does not make gravity an inverse-cube law. The vector in the numerator contributes one power of distance, leaving the force magnitude proportional to $1/d^2$.

Newton’s third law gives the force on body 2:

$${\mathbf{F}}_2=-{\mathbf{F}}_1.$$

It has the same magnitude and points towards body 1. These are forces on different bodies, not two forces that cancel on either individual body.

A reliable way to check a gravitational force formula is to temporarily place the bodies on a line. Suppose body 1 is to the left of body 2. Then ${\mathbf{r}}_2-{\mathbf{r}}_1$ points right, and the formula for ${\mathbf{F}}_1$ points right too. Its negative, ${\mathbf{F}}_2$, points left. Both bodies are pulled towards one another.

---

Deriving the coupled equations

Apply Newton’s second law to each body separately:

$$m_1{\ddot{\mathbf{r}}}_1={\mathbf{F}}_1,m_2{\ddot{\mathbf{r}}}_2={\mathbf{F}}_2.$$

Substituting the gravitational forces and dividing by the respective masses yields the governing equations:

$$\boxed{{\ddot{\mathbf{r}}}_1=G{m}_2\frac{{\mathbf{r}}_2-{\mathbf{r}}_1}{\Vert {\mathbf{r}}_2-{\mathbf{r}}_1{\Vert }^3},}$$ $$\boxed{{\ddot{\mathbf{r}}}_2=G{m}_1\frac{{\mathbf{r}}_1-{\mathbf{r}}_2}{\Vert {\mathbf{r}}_2-{\mathbf{r}}_1{\Vert }^3}.}$$

Equivalently, using $\mathbf{d}={\mathbf{r}}_2-{\mathbf{r}}_1$,

$${\ddot{\mathbf{r}}}_1=G{m}_2\frac{\mathbf{d}}{d^3},{\ddot{\mathbf{r}}}_2=-G{m}_1\frac{\mathbf{d}}{d^3}.$$

Notice the apparent exchange of masses: body 1’s acceleration depends on $m_2$, while body 2’s acceleration depends on $m_1$. Each body’s own mass cancels between its gravitational force and $m\mathbf{a}$. That cancellation does not make the two accelerations equal. The forces have equal magnitude, but the less massive body accelerates more.

The equations are coupled because neither acceleration can be evaluated from that body’s position alone. To compute ${\ddot{\mathbf{r}}}_1$, you need both ${\mathbf{r}}_1$ and ${\mathbf{r}}_2$; the same is true of ${\ddot{\mathbf{r}}}_2$. As the bodies move, their shared separation changes, changing both accelerations.

This also clarifies what “isolated” means here. It does not mean that either body experiences zero force. It means no external forces act on the two-body system. Mutual gravity remains in the model. Atmospheric drag, other gravitating bodies, and departures from the point-mass approximation do not.

---

Checks that catch modelling mistakes

Before trusting an ODE implementation, test the equations against simple physical requirements. These checks are particularly useful because a sign error can still produce a smooth-looking trajectory.

First, the force magnitude must have the correct units:

$$[G]=\mathrm{lengt}{\mathrm{h}}^3\mathrm{mas}{\mathrm{s}}^{-1}\mathrm{tim}{\mathrm{e}}^{-2}.$$

Consequently, $G{m}_2\mathbf{d}/d^3$ has units of length divided by time squared, as an acceleration should. Using consistent units matters: positions in kilometres require a value of $G$ expressed using kilometres, not its numerical value in SI units.

Second, the forces must be equal and opposite. Multiplying the two acceleration equations by their respective masses gives

$$m_1{\ddot{\mathbf{r}}}_1+m_2{\ddot{\mathbf{r}}}_2=0.$$

This identity follows directly from the paired force formulas. It is a powerful implementation check: if your computed gravitational forces do not sum to zero, you have broken the action–reaction symmetry of this isolated model. In the next lesson, this identity will help us understand the motion of the system as a whole.

Third, compare the equations with the fixed-centre approximation used earlier. If $m_1$ is much larger than $m_2$, then body 1’s acceleration, proportional to $m_2$, is very small, while body 2’s acceleration, proportional to $m_1$, is not. Treating body 1 as fixed can then be an excellent approximation over a suitable interval. Small is not the same as zero, however: the coupled equations allow body 1 to respond.

At the opposite extreme, imagine equal masses initially separated by a distance $d$, with no other forces. Each experiences a force of magnitude $G{m}_1{m}_2/d^2$, and their acceleration magnitudes are equal. Declaring either mass permanently stationary would impose an artificial constraint absent from the two-body model.

These cases test different aspects of the formulation. The large-mass-ratio case explains why the earlier model often works; the equal-mass case exposes why it is not the general model.

---

From two acceleration equations to one initial-value problem

You have already converted a second-order gravitational equation into a first-order state equation. The same construction now includes both bodies. Define their velocities by

$${\mathbf{v}}_1={\dot{\mathbf{r}}}_1,{\mathbf{v}}_2={\dot{\mathbf{r}}}_2,$$

and stack the variables in a specified order:

$$\mathbf{y}=\left[\begin{array}{c}{\mathbf{r}}_1\\ {\mathbf{r}}_2\\ {\mathbf{v}}_1\\ {\mathbf{v}}_2\end{array}\right].$$

In three spatial dimensions, $\mathbf{y}$ has twelve scalar components. Its first-order equation is

$$\boxed{\dot{\mathbf{y}}=\left[\begin{array}{c}{\mathbf{v}}_1\\ {\mathbf{v}}_2\\ G{m}_2({\mathbf{r}}_2-{\mathbf{r}}_1)/d^3\\ -G{m}_1({\mathbf{r}}_2-{\mathbf{r}}_1)/d^3\end{array}\right],d=\Vert {\mathbf{r}}_2-{\mathbf{r}}_1\Vert .}$$

The initial-value problem is completed by specifying, at some initial time $t_0$,

$${\mathbf{r}}_1(t_0),{\mathbf{r}}_2(t_0),{\mathbf{v}}_1(t_0),{\mathbf{v}}_2(t_0),$$

along with $m_1$, $m_2$, and a consistent value of $G$. The initial positions must be distinct. There is no requirement that either initial position or velocity be zero; the inertial origin need not be attached to a body.

Here is a direct translation of that right-hand side into Python. It evaluates derivatives at a supplied state; it does not yet advance the state through time.

import numpy as np

def two_body_rhs(t, state, m1, m2, G):
    """Derivative of [r1, r2, v1, v2] in a common inertial frame."""
    r1 = state[0:3]
    r2 = state[3:6]
    v1 = state[6:9]
    v2 = state[9:12]

    d = r2 - r1
    d_squared = np.dot(d, d)
    if d_squared == 0.0:
        raise ValueError("Point-mass gravity is undefined at collision.")

    inverse_d_cubed = 1.0 / (d_squared * np.sqrt(d_squared))
    a1 = G * m2 * d * inverse_d_cubed
    a2 = -G * m1 * d * inverse_d_cubed

    return np.concatenate((v1, v2, a1, a2))

The argument t is present because ODE solvers commonly expect it. This particular right-hand side has no explicit time dependence: time affects acceleration through the changing positions in state.

A practical discipline is to keep the state ordering written beside the code. If positions and velocities are packed in one order but unpacked in another, the resulting numerical trajectory may be wrong without producing a Python error. The paired formulas also make two inexpensive checks available before integration: a1 must point towards r2, and m1 * a1 + m2 * a2 should be approximately zero up to floating-point rounding.

---

Connecting to last lesson’s scaling

The dimensional two-body equations can use the characteristic scales introduced last lesson. Choose a length $L_c$, a mass $M_c$, and

$$t_c=\sqrt{\frac{L_c^3}{G{M}_c}}.$$

Write ${\mathbf{r}}_i=L_c{\mathbf{\rho }}_i$ and $t=t_c\tau$. Substitution into the two acceleration equations gives

$$\frac{d^2{\mathbf{\rho }}_1}{d{\tau }^2}=\frac{m_2}{M_c}\frac{{\mathbf{\rho }}_2-{\mathbf{\rho }}_1}{\Vert {\mathbf{\rho }}_2-{\mathbf{\rho }}_1{\Vert }^3},$$ $$\frac{d^2{\mathbf{\rho }}_2}{d{\tau }^2}=-\frac{m_1}{M_c}\frac{{\mathbf{\rho }}_2-{\mathbf{\rho }}_1}{\Vert {\mathbf{\rho }}_2-{\mathbf{\rho }}_1{\Vert }^3}.$$

For two freely moving bodies, $M_c=m_1+m_2$ is a useful choice. The dimensionless equations then retain two mass fractions whose sum is one. Unlike the fixed-centre model, the individual bodies’ equations need not each have a gravitational coefficient of one: their different accelerations are part of the physics.

This is a useful distinction to preserve in the simulator. Scaling changes the numbers used to describe a system, not which bodies are allowed to move.

---

Key takeaways

A Newtonian two-body model uses two positions in one inertial frame and one shared separation vector. Gravity produces equal-and-opposite forces, but generally unequal accelerations:

$${\ddot{\mathbf{r}}}_1=G{m}_2\frac{{\mathbf{r}}_2-{\mathbf{r}}_1}{d^3},{\ddot{\mathbf{r}}}_2=-G{m}_1\frac{{\mathbf{r}}_2-{\mathbf{r}}_1}{d^3},d=\Vert {\mathbf{r}}_2-{\mathbf{r}}_1\Vert .$$

These are coupled equations: both positions are needed to calculate either acceleration. Four initial vectors—two positions and two velocities—turn them into an initial-value problem ready for numerical integration.

Next, we will use the system’s collective motion and the separation between the bodies to simplify these equations, without reverting to the assumption that one body is fixed.