I am trying to simulate the three body problem, my current code is:
// Mass of planets
const mass = {m_1: 1, m_2: 1, m_3: 1}
function vector_sum(p1, p2){
let vec = p1.slice()
for(let i=0; i< p1.length; i++){
vec[i] += p2[i]
}
return vec
}
function const_multiplication(c, p){
let vec = p.slice()
for(let i=0; i<p.length; i++){
vec[i] *= c
}
return vec
}
function solve(){
// Starting coordinates for planets
var p1_start = {x: -1, y: 0}
var v1_start = {x: 0.306893, y: 0.125507}
var p2_start = {x: 1, y: 0}
var v2_start = {x: 0.306893, y: 0.125507}
var p3_start = {x: 0, y: 0}
var v3_start = {x: -0.613786, y: -0.251014}
let y = [p1_start.x, p1_start.y, v1_start.x, v1_start.y,
p2_start.x, p2_start.y, v2_start.x, v2_start.y, p3_start.x, p3_start.y, v3_start.x, v3_start.y]
y = [-1, 0, 0.306893, 0.125507, 1, 0, 0.306893, 0.125507, 0, 0, -0.613786, -0.251014]
function f(yprime, y, t){
// We copy the values know from state into yprime
yprime[0] = y[2]
yprime[1] = y[3]
yprime[4] = y[6]
yprime[5] = y[7]
yprime[8] = y[10]
yprime[9] = y[11]
yprime[2] = (y[8]-y[0]) / (Math.sqrt((y[8]-y[0])**2 + (y[9]-y[1])**2))**3 * mass.m_3 +
(y[4]-y[0]) / (Math.sqrt((y[4]-y[0])**2 + (y[5]-y[1])**2))**3 * mass.m_2
yprime[3] = (y[9]-y[1]) / (Math.sqrt((y[8]-y[0])**2 + (y[9]-y[1])**2))**3 * mass.m_3 +
(y[5]-y[1]) / (Math.sqrt((y[4]-y[0])**2 + (y[5]-y[1])**2))**3 * mass.m_2
yprime[6] = (y[0]-y[4]) / (Math.sqrt((y[0]-y[4])**2 + (y[1]-y[5])**2))**3 * mass.m_1 +
(y[8]-y[4]) / (Math.sqrt((y[8]-y[4])**2 + (y[9]-y[5])**2))**3 * mass.m_3
yprime[7] = (y[1]-y[5]) / (Math.sqrt((y[0]-y[4])**2 + (y[1]-y[5])**2))**3 * mass.m_1 +
(y[9]-y[5]) / (Math.sqrt((y[8]-y[4])**2 + (y[9]-y[5])**2))**3 * mass.m_3
yprime[10] = (y[0]-y[8]) / (Math.sqrt((y[0]-y[8])**2 + (y[1]-y[9])**2))**3 * mass.m_1 +
(y[4]-y[8]) / (Math.sqrt((y[4]-y[8])**2 + (y[5]-y[9])**2))**3 * mass.m_2
yprime[11] = (y[1]-y[9]) / (Math.sqrt((y[0]-y[8])**2 + (y[1]-y[9])**2))**3 * mass.m_1 +
(y[5]-y[9]) / (Math.sqrt((y[4]-y[8])**2 + (y[5]-y[9])**2))**3 * mass.m_2
}
// Parameters
const delta_t = 0.01
const n = 10000
function runge_Kutta(t, y, f, delta_t){
let k1 = Array(12).fill(0);
let k2 = Array(12).fill(0);
let k3 = Array(12).fill(0);
let k4 = Array(12).fill(0);
// K1 = f(t, p1, p2, p3)
f(k1, y)
// K2 = f(t + delta_t/2, (p1, p2, p3) + delta_t/2*k1)
f(k2, vector_sum(y, const_multiplication(delta_t/2, k1)))
// K3 = f(t + delta_t/2, (p1, p2, p3) + delta_t/2*k2)
f(k3, vector_sum(y, const_multiplication(delta_t/2, k2)))
// K4 = f(t + delta_t, (p1, p2, p3) + delta_t*k3)
f(k4, vector_sum(y, const_multiplication(delta_t, k3)))
// v(t + delta_t) = v(t) + delta_t/6*(k1 + 2*k2 + 2*k3 + k4)
return vector_sum(y, const_multiplication(delta_t/6, vector_sum(k1, vector_sum(const_multiplication(2, k2), vector_sum(const_multiplication(2, k3), k4)))))
}
// Trajectory arrays
p1 = [p1_start]
v1 = [v1_start]
p2 = [p2_start]
v2 = [v2_start]
p3 = [p3_start]
v3 = [v3_start]
// Evolution of the system
for(let i = 0; i < n; i++){
y = runge_Kutta(i, y, f, delta_t)
p1[i+1] = {x: y[0], y: y[1]}
p2[i+1] = {x: y[4], y: y[5]}
p3[i+1] = {x: y[8], y: y[9]}
}
}
The code runs fine but the results don’t make sense. The initial conditions set in the code should yield a periodic orbit representing a butterfly. Instead when I check p1[1000], I get around {x: -8, y: 19} which means is not periodic because according to the solution: https://observablehq.com/@rreusser/periodic-planar-three-body-orbits#computeShape
x and y never leave the range of (-2,2).
I believe the error is withing Runge Kutta but I don’t seem to find it.
Any help would be appreciated
