Animating Chaotic Systems Using Matplotlib - Alpha Centauri 3-Body Problem

Chaos theory resurfaces as a popular modeling tool for complex systems from time to time. Physics makes good use of it; financial market modeling tries to adapt it and finds the constant feedback and higher-order interactions that may give rise to the apparent randomness of the market. In this post, we will visualize this chaotic behavior using matplotlib. The three-body problem, where Newton´s law cannot predict the motion of three interacting point masses, will serve as an initial problem. This problem can be easily extended to any other system type by altering the descriptive parameters.

For this demonstration, we need NumPy, scipy, and matplotlib 3D plotting tools:

import numpy as np
import scipy as sci
# Matplotlib plotting tools:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation

We define the physical guiding constants, in this case, the universal gravitational constant:

# Define universal gravitation constant
G = 6.67408e-11 #N-m2/kg2

We also need reference constants to model the motion of a three-star system; we use the mass of our sun as the reference value for star mass and the parameters for a known three-star system as Alpha-Centauri. We take the chance to celebrate the awesome Sid Meier´s Alpha-Centauri game, which, in part, inspires everything we do.

# Mass of the sun.
m_nd = 1.989e30

# Define masses:
# Rigil Kentaurus:
m1 = 1.1
# Toliman:
m2 = 0.907
# Mass of Proxima Centauri:
m3 = 0.256

# Distance between stars in the Alpha-Centauri System:
r_nd = 5.326e12 
# Orbital Period of Alpha Centauri:
t_nd = 79*365*24*3600
# Orbital velocity:
v_nd = 23000

# Net constants
K1 = G*t_nd*m1*m_nd/((r_nd**2)*v_nd)
K2 = v_nd*t_nd/r_nd

The condition for the starting positions for the stars in our 3-body start system will be randomly set. We will call the random position generation function for each of the three stars. The velocities are set to fixed values as these have to be sufficiently small not to cause the stars to escape each other´s gravitational pull. With these data points, we can also find the center of mass of the system and its velocity:

# Define initial position vectors
def roll():
  r = 2 * np.random.random_sample(3) - 1
  return r

r1 = roll()
r2 = roll()
r3 = roll()

# Define initial velocities
v1 = np.array([0.01, 0.01, 0])
v2 = np.array([-0.05, 0, -0.01])
v3 = np.array([0, -0.01, 0])

# Find Centre of Mass
r_com = (m1*r1+m2*r2+m3*r3)/(m1+m2+m3)

# Find velocity of COM
v_com = (m1*v1+m2*v2+m3*v3)/(m1+m2+m3)

The equations for the motion of the stars are defined in the function we will solve. The equations represent the models shown here.

We are trying to minimize the repetition of expressions by defining nested functions. This is not a good style for production code; in this case, it facilitates understanding the individual steps of the calculation.

def ThreeBodyEquations(w, t, G, m1, m2, m3):
    DIMS = 3

    # We can insert a 0 value at "V0" for readability.
    # This is done to access v values with a
    # 1-based index only and is not a good practice.

    v = [0]

    for i in range(DIMS*2):
      v.append(w[i*DIMS:(i+1)*DIMS])

    norm = np.linalg.norm
    r12 = norm(v[2] - v[1])
    r13 = norm(v[3] - v[1])
    r23 = norm(v[3] - v[2])

    def dvdt(K1, m1, m2, v1, v2, v3, r1, r2):
      A = m1*(v2-v1)/r1**3 
      B = m2*(v3-v1)/r2**3
      C = K1
      return C*(A+B)

    dv1dt = dvdt(K1, m2, m3, v[1], v[2], v[3], r12, r13)
    dv2dt = dvdt(K1, m1, m3, v[2], v[1], v[3], r12, r23)
    dv3dt = dvdt(K1, m1, m2, v[3], v[1], v[2], r13, r23)

    dr1dt=K2*v[4]
    dr2dt=K2*v[5]
    dr3dt=K2*v[6]

    concat = np.concatenate

    r12_d = concat((dr1dt, dr2dt))
    r_d = concat((r12_d, dr3dt))
    v12_d = concat((dv1dt, dv2dt))
    v_d = concat((v12_d, dv3dt))
    d = concat((r_d, v_d))

    return d

With the initial parameters set, we solve the position and velocity vector equations using scipy ordinary differential equations solver.

# Package initial parameters
init_params = np.array([r1, r2, r3, v1, v2, v3])

#250 orbital periods and 20000 points
time_span = np.linspace(0, 250, 20000) 

# Run the ODE solver
import scipy.integrate
three_body_sol = sci.integrate.odeint(ThreeBodyEquations,
                                    init_params.flatten(),
                                    time_span,
                                    args=(G, m1, m2, m3))

The coordinates in (x,y,z) format can be extracted now for plotting purposes:

def extract_coords(r):
  l = []
  for i in range(3):
    l.append(r[:, i])
  return np.array(l)

coords_1 = extract_coords(three_body_sol[:, :3])
coords_2 = extract_coords(three_body_sol[:, 3:6])
coords_3 = extract_coords(three_body_sol[:, 6:9])

To generate an animated plot, we have to define a canvas plot, shown as an example, and an update function. Matplotlib will generate updated plots for each frame we want to show in the animation. We are not doing anything more complex than plotting an ever-increasing set of data, in this case, showing a time series:

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation

def update_lines(num, dataLines, lines):
    for line, data in zip(lines, dataLines):
        line.set_data(data[0:2, :num])
        line.set_3d_properties(data[2, :num])
    return lines

# Attaching 3D axis to the figure
fig = plt.figure()
ax = p3.Axes3D(fig)

# Data with position coordinates:
data = [coords_1, coords_2, coords_3]

# Creating line objects:
lines = [ax.plot(dat[0, 0:1], dat[1, 0:1], dat[2, 0:1])[0] for dat in data]

# Setting the axes properties
ax_limits = [-1, 1]
ax.set_xlim3d(ax_limits)
ax.set_xlabel('X')

ax.set_ylim3d(ax_limits)
ax.set_ylabel('Y')

ax.set_zlim3d(ax_limits)
ax.set_zlabel('Z')

ax.set_title('3D 3-Body Problem')

# Creating the Animation object
line_ani = animation.FuncAnimation(fig,
                                   update_lines,
                                   fargs=(data, lines),
                                   frames = 500,
                                   interval = 10,
                                   blit=False,
)

plt.show()

Finally, the animation can be displayed as an HTML animation; initial parameters are random; when running the notebook, the trajectories will change. Here we have some examples:

from IPython.display import HTML
HTML(line_ani.to_html5_video())

Do not hesitate to contact us if you require quantitative model development, deployment, verification, or validation. We will also be glad to help you with your machine learning or artificial intelligence challenges when applied to asset management, automation, or intelligence gathering from satellite, drone, or fixed-point imagery.

The notebook with this model is here.