Celestial OpenSeesing

Original Post - 14 Sep 2025 - Michael H. Scott

The three-body problem is an application of Newton’s laws of motion where each of three bodies moves through space according to gravitational forces exerted by the other two bodies. Despite its simplicity, there is no general closed-form solution for the trajectories of the three bodies.

However, you can approximate the trajectories with numerical methods. And since it’s so good at numerical solutions of Newton’s laws, why not use OpenSees to simulate the three-body problem?

The force body j applies to body i is based on Newton’s law of universal gravitation

where m_i and r_i are the mass and position (vector) of body i and G is the gravitational constant, whose true value is G=6.674e-11 m³/(kg sec²). For demonstration of the three-body problem, G can take on arbitrary values, typically G=1.

We can put the gravitational force calculations for three bodies into a Python function.

import numpy as np

def three_body_forces(m1,m2,m3,r1,r2,r3):
    '''
    m1,m2,m3 - scalar masses
    r1,r2,r3 - position vectors [x,y,z]
    '''
    d12 = np.linalg.norm(r1-r2)
    d13 = np.linalg.norm(r1-r3)
    f1 = -G*m1 * (m2*(r1-r2)/d12**3 + m3*(r1-r3)/d13**3)

    d23 = np.linalg.norm(r2-r3)
    f2 = -G*m2 * (m1*(r2-r1)/d12**3 + m3*(r2-r3)/d23**3)

    f3 = -G*m3 * (m1*(r3-r1)/d13**3 + m2*(r3-r2)/d23**3)

    return f1,f2,f3

In some implementations of the three-body problem, it makes more sense for this function to return the accelerations rather than the forces. With OpenSees, I found it easier to work directly with the forces rather than the accelerations.

At each simulation time step, we can apply the gravitational forces in a load pattern, perform an analysis, then remove the load pattern and repeat this process at the next time step. This OpenSees simulation is one of the rare cases where the Diagonal solver is OK to use because the equations of motion are uncoupled.

import openseespy.opensees as ops

G = 1.0 # gravitational constant

# Masses
m1 =  98.0
m2 = 196.0
m3 = 294.0

# Initial conditions
r1o = np.array([-10.0,10.0,-11.0])
r2o = np.array([  0.0, 0.0,  0.0])
r3o = np.array([ 10.0,10.0, 12.0])

v1o = np.array([-3.0,0.0,0.0])
v2o = np.array([ 0.0,0.0,0.0])
v3o = np.array([ 3.0,0.0,0.0])

# Time stepping
dt = 0.001
Tfinal = 200
Nsteps = int(Tfinal/dt)

ops.wipe()
ops.model('basic','-ndm',3,'-ndf',3)

# Nodes
ops.node(1,*r1o); ops.mass(1,m1,m1,m1)
ops.node(2,*r2o); ops.mass(2,m2,m2,m2)
ops.node(3,*r3o); ops.mass(3,m3,m3,m3)

# Set initial velocity
for dof in range(3):
    ops.setNodeVel(1,dof+1,v1o[dof],'-commit')
    ops.setNodeVel(2,dof+1,v2o[dof],'-commit')
    ops.setNodeVel(3,dof+1,v3o[dof],'-commit')

ops.system('Diagonal')
ops.integrator('Newmark',0.5,1.0/6)
ops.algorithm('Linear')
ops.analysis('Transient','-noWarnings')

ops.timeSeries('Constant',1)

# Body locations
r1 = r1o
r2 = r2o
r3 = r3o

for i in range(Nsteps):

    f1,f2,f3 = three_body_forces(m1,m2,m3,r1,r2,r3)

    ops.pattern('Plain',1,1)
    ops.load(1,*f1)
    ops.load(2,*f2)
    ops.load(3,*f3)

    ops.analyze(1,dt)

    r1 = r1o + ops.nodeDisp(1)
    r2 = r2o + ops.nodeDisp(2)
    r3 = r3o + ops.nodeDisp(3)

    ops.remove('loadPattern',1)

The masses and initial conditions for this demonstration are from this page, which was very helpful to my understanding of the three-body problem and its implementation.

The trajectories of the three bodies computed with OpenSees are shown below.

Let’s be honest–this simulation is only using OpenSees for time integration. No actual changes to the OpenSees source code were required. This approach is in the same superficial-to-OpenSees ballpark as doing state determination and assembly on CPUs, moving A and b to a GPU to solve for x, then shipping x back to the CPUs to update the model for the next time step.

Anyway, numerical solution of the three-body problem is very sensitive to initial conditions and the chosen time integration method–there’s no stiffness or damping to rein things in. The sensitivity to numerical approximations makes it difficult to simulate motions over long periods of time.

For example, below, the motion on the left uses Newmark time integration with the linear acceleration approximation (\(\beta\)=1/6, same method used for the animation above) while the motion on the right uses the constant acceleration approximation (\(\beta\)=1/4). The computed trajectories are similar until about half way through the simulation. Then, masses 1 and 3 spiral off in one direction while mass 2 escapes in the opposite direction.

Numerical solutions to the three-body problem are based on anything from forward Euler to fourth order Runge-Kutta and beyond. There’s like 50 other time integration methods, both explicit and implicit, implemented in OpenSees, so you can play around with those.

And if you want to model open seas of stars, you can extend the script to solve the N-body problem. But beware, there are about 200 billion stars in the Milky Way galaxy alone.

This catnip of a post topic was suggested by a reader of the blog. If you have any classical non-structural mechanics problems you think OpenSees can solve, let me know.

OpenSees Cloud

OpenSees AMI

Celestial OpenSeesing