I need to know how to solve two-body problem by solving a system of first order equation derived from the equation below.

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

How do I go about this, and how would I then use this to output a trajectory in MATLAB?

  • 1
    $\begingroup$ +1 We don't have a Matlab tag but we do have a Python tag and they are similar looking enough languages (superficially) that hopefully the tag will attract programmers of both. I've added MathJax formatting standard for most Stack Exchange sites. $\endgroup$
    – uhoh
    Commented Nov 13, 2020 at 13:08
  • 2
    $\begingroup$ I've posted an answer using free, open source, friendly, popular, fun and well-documented Python. Perhaps it might help sway you into trying to move towards it and away from expensive and black-box MATLAB. $\endgroup$
    – uhoh
    Commented Nov 13, 2020 at 13:43

1 Answer 1


There are several ways to do this. The easiest and most straightforward is to break it into two sets by including velocity as a variable, and solve together.

Instead of a single second order differential equation

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

We can solve the following pair of first order differential equations in parallel

$$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$$ $$\dot{\mathbf{r}} = \mathbf{v}$$

using various simple methods including standard libraries or homegrown implementations of Runge-Kutta including my favorite simple RK4/5 with variable step size.

It's wonderful and really educational to code that once for yourself and appreciate the task first before using standard libraries from then on.

For more on errors using differential equation solvers and for a way to test them yourself, see my question in Math SE: need to understand better ODE solution accuracy vs numerical precision (which I should have asked in Computational Science SE)

For programming purposes you can write $\mathbf{r} / r^3 $ as r * (r**2).sum()**-1.5

From this answer you can see a 2D implementation using not only the monopole $1/r^2$ gravity term but the additional quadrupole $J_2$ term for Earth's oblate shape and field. For more on that see this answer to Trouble deriving rectangular components of acceleration of satellite in orbit around Earth with J2 consideration.

See also a similar solution in this answer to How to determine an orbit of a satellite for a collision detection?.

Also, you can consider going unit-less by using $\mu = 1$ and period $T = 2 \pi$ and semi-major axis $a=1$. Some integrators will handle that a little better.

Blue line is correct, red line has J2 ten times its real value to show exaggerated apsidal precession for fun.

apsidal precession with normal and 10x stronger J2

def deriv(X, t):

    x, v = X.reshape(2, -1)

    acc0  = -GMe * x * ((x**2).sum())**-1.5
    acc2  = -1.5 * GMe * J2 * Re**2 * x * ((x**2).sum())**-2.5

    return np.hstack([v, acc0 + acc2])

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# David Hammen's nice table https://physics.stackexchange.com/a/141981/83380
# See http://www.iag-aig.org/attach/e354a3264d1e420ea0a9920fe762f2a0/51-groten.pdf
# https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere

GMe = 3.98600418E+14  # m^3 s^-2
J2e = 1.08262545E-03  # unitless
Re  = 6378136.3 # meters

X0 = np.hstack([6778000.0, 0.0, 0.0, 10000.])  # x, y, vx, vy

time = np.arange(0, 300001, 100)

J2 = J2e  # correct J2
answerJ2, info = ODEint(deriv, X0, time, full_output=True)

J2 = 10*J2e # 10x larger J2
answer10xJ2, info = ODEint(deriv, X0, time, full_output=True)

if True:
    x, y = answerJ2.T[:2]
    plt.plot(x, y, '-b')
    x, y = answer10xJ2.T[:2]
    plt.plot(x, y, '-r')
    plt.plot([0], [0], 'or')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.