Three body problem using Runge Kutta on JS

53 Views Asked by At

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

0

There are 0 best solutions below