integrate 3 body system in python
Question: This is more mathemathical question.

We have 3 body system, with known initial parametrs like positions, masses,
velocityes. This make system like  
[![enter image description
here](http://i.stack.imgur.com/VyZjB.gif)](http://i.stack.imgur.com/VyZjB.gif)

i and j = 1,2,3  
So question is how to deal with this, to feed this system to **scipy odeint**
?

**UPDATE**  
I wrote code

    
    
    from scipy import linspace
    from scipy.integrate import odeint
    
    def rf(i, j, r):
        x1, y1, z1 = r[0]
        x2, y2, z2 = r[1]
        x3, y3, z3 = r[2]
        r12 = ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**2
        r13 = ((x1-x3)**2+(y1-y3)**2+(z1-z3)**2)**2
        r23 = ((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)**2
        if i == 1:
            if j == 2:
                r = r12
            elif j == 3:
                r = r13
        if i == 2:
            if j == 1:
                r = r12
            elif j == 3:
                r = r23
        if i == 3:
            if j == 1:
                r = r13
            elif j == 2:
                r = r23
        return r
    
    
    def dotf(r, t0):
        x1, y1, z1 = r[0]
        x2, y2, z2 = r[1]
        x3, y3, z3 = r[2]
        x = [x1, x2, x3]
        y = [y1, y2, y3]
        z = [z1, z2, z3]
    
        k2 = 6.67e-11
        m = [2e6, 2.2e7, 0.1e3]
    
        for i in range(1, 3, 1):
            xs, ys, zs = 0, 0, 0
            for j in range(1, 3, 1):
                if i != j:
                    r = rf(i, j, r)
                    xs = xs + m[j]*(x[i]-x[j])/r
                    ys = ys + m[j]*(y[i]-y[j])/r
                    zs = zs + m[j]*(z[i]-z[j])/r
            x[i] = -1 * k2 * xs
            y[i] = -1 * k2 * ys
            z[i] = -1 * k2 * zs
        return [[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]]
    
    
    t = linspace(0, 50, 1e4)
    r1 = [1, 2, 1]
    r2 = [0.5, 0.1, 3]
    r3 = [0.6, 1, 1.5]
    r = [r1, r2, r3] 
    
    u = odeint(dotf, r, t)
    

  
and get output error:  

    
    
    Traceback (most recent call last):
      File "3b.py", line 73, in <module>
        u = odeint(dotf, r, t)
      File "/usr/local/lib/python2.7/dist-packages/scipy/integrate/odepack.py", line 148, in odeint
        ixpr, mxstep, mxhnil, mxordn, mxords)
    ValueError: object too deep for desired array
    


Answer: I've corrected two obvious bugs in your code; it runs, but I'm not sure it's
correct. Code at the end.

The first bug is that `odeint` wants to deal in vectors for state and state
derivative. By vector I mean rank-1 arrays; you were submitting a matrix
(rank-2 array) as your initial condition. I've changed the `r` assignment,
possibly confusingly, as `r = r1 + r2 + r3` \- the `+` operator on lists is
concatenation, though on numpy arrays it is element-wise addition.

This change meant I had to change the assignment to `x1`, etc., in `dotf` and
`rf`. I also change the return from `dotf` to be a vector.

The second bug is that you used `r` to mean two different things: first, the
system state vector, and second the return from `rf`. I changed the second one
to be `rr`.

The solution appears to be unstable; I don't know if this is plausible. I
imagine there are other bugs, but you at least have something that runs now.

My suggestion is starting with something simpler that you know the solution
to, e.g., stable first-order linear system, or a stable underdamped second-
order system, or Lorenz if you want something chaotic.

See <http://docs.scipy.org/doc/> for Scipy documentation, e.g., I used
<http://scipy-0.13.0/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint>
for `odeint` for Scipy 0.13.0 (which is what comes with my Ubuntu 14.04
system).

    
    
    from scipy import linspace
    from scipy.integrate import odeint
    
    def rf(i, j, r):
        x1, y1, z1 = r[:3]
        x2, y2, z2 = r[3:6]
        x3, y3, z3 = r[6:]
        r12 = ((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**2
        r13 = ((x1-x3)**2+(y1-y3)**2+(z1-z3)**2)**2
        r23 = ((x2-x3)**2+(y2-y3)**2+(z2-z3)**2)**2
        if i == 1:
            if j == 2:
                r = r12
            elif j == 3:
                r = r13
        if i == 2:
            if j == 1:
                r = r12
            elif j == 3:
                r = r23
        if i == 3:
            if j == 1:
                r = r13
            elif j == 2:
                r = r23
        return r
    
    
    def dotf(r, t0):
        x1, y1, z1 = r[:3]
        x2, y2, z2 = r[3:6]
        x3, y3, z3 = r[6:]
        x = [x1, x2, x3]
        y = [y1, y2, y3]
        z = [z1, z2, z3]
    
        k2 = 6.67e-11
        m = [2e6, 2.2e7, 0.1e3]
    
        for i in range(1, 3, 1):
            xs, ys, zs = 0, 0, 0
            for j in range(1, 3, 1):
                if i != j:
                    rr = rf(i, j, r)
                    xs = xs + m[j]*(x[i]-x[j])/rr
                    ys = ys + m[j]*(y[i]-y[j])/rr
                    zs = zs + m[j]*(z[i]-z[j])/rr
            x[i] = -1 * k2 * xs
            y[i] = -1 * k2 * ys
            z[i] = -1 * k2 * zs
        return [x1,y1,z1, x2,y2,z2, x3,y3,z3]
    
    
    t = linspace(0, 50, 1e4)
    r1 = [1, 2, 1]
    r2 = [0.5, 0.1, 3]
    r3 = [0.6, 1, 1.5]
    r = r1 + r2 + r3
    
    u = odeint(dotf, r, t)
    
    
    # uncomment to plot
    # from matplotlib import pyplot
    #
    # pyplot.plot(t,u)
    # pyplot.show()
    

