I have programmed a 3 body problem, that can in theory work for more than 3 bodies. I apologize in advance for the length of the code that I posted, I don't really know how much one would need to help.

My problem in a nutshell is that when I switched the program from 2D to 3D by adding a component to my initial vectors and making my ListPlot 3D, it didn't work any more and I couldn't see what was even happening. (I just want a way for me to see the path each of my particles takes in 3D)

On a side note... how do I create my Orbit1, Orbit2 and Orbit3 in a loop? so that I don't have to just keep adding terms when I add a particle.

Off[General::spell, General::spell1, General::infy];
dt = 0.01;
bodies = 3;
rnum := RandomReal[{-1.00, 1.00}];
rnum2 := RandomReal[{-0.07, 0.07}];

PosVel = Table[0, {0}, {j, 2}];
   Pos = {rnum, rnum, rnum},
   Vel = {rnum2, rnum2, rnum2},
   AppendTo[PosVel, {Pos, Vel}]
   }, {j, 1, bodies}];

r = Table[0, {i, bodies}, {j, bodies}];
rmag = Table[0, {i, bodies}, {j, bodies}];
rsq = Table[0, {i, bodies}, {j, bodies}];

   r[[i, j]] = (PosVel[[i, 1]] - PosVel[[j, 1]]),
   rmag[[i, j]] = Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])],
   rsq[[i, j]] = (Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])])^2
   }, {j, 1, bodies}, {i, 1, bodies}]

f = -1/rmag^3;
  f[[i, i]] = 0
  }, {i, 1, bodies}]

Force = Table[0, {i, bodies}];
  Force[[i]] = Sum[f[[i, j]]*r[[i, j]], {j, 1, bodies}]
  }, {i, 1, bodies}]

n = 3000;
Orbit = Table[0, {i, n}, {j, bodies}];

  PosVel[[i, 2]] = PosVel[[i, 2]] + Force[[i]]*dt,
  PosVel[[i, 1]] = PosVel[[i, 1]] + PosVel[[i, 2]]*dt,

    r[[i, j]] = (PosVel[[i, 1]] - PosVel[[j, 1]]),
    rmag[[i, j]] = Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])],
    rsq[[i, j]] = (Norm[(PosVel[[i, 1]] - PosVel[[j, 1]])])^2
    }, {j, 1, bodies}, {i, 1, bodies}],
  f = -1/rmag^3,
    f[[i, i]] = 0
    }, {i, 1, bodies}],
    Force[[i]] = Sum[f[[i, j]]*r[[i, j]], {j, 1, bodies}]
    }, {i, 1, bodies}],
  Orbit[[j, i]] = PosVel[[i, 1]]
  }, {i, 1, bodies}, {j, 1, n}]

Orbit1 = Table[Orbit[[i, 1]], {i, 1, n}];
Orbit2 = Table[Orbit[[i, 2]], {i, 1, n}];
Orbit3 = Table[Orbit[[i, 3]], {i, 1, n}];

 ListPlot3D[{Orbit1[[1 ;; j]], Orbit2[[1 ;; j]], Orbit3[[1 ;; j]]},
  Mesh -> None,
  Epilog -> {PointSize -> 0.05, Blue, 
    Point[{Orbit1[[j]], Orbit2[[j]], Orbit3[[j]]}]}], {j, 1, n, 1}]
  • 1
    $\begingroup$ I think the plot command you are looking for is ListPointPlot3D. $\endgroup$ Oct 19, 2014 at 22:28
  • 1
    $\begingroup$ As for the second part of your question. You could define the creation of one of those orbits like: orbit[x_] := Table[Orbit[[i, x]], {i, 1, n}]; And then do it for every body using map: Map[Orbit[#] &, Range[bodies]] I think in general you could speed and clean up your code by using constructs like Map instead of Do loops. I'd write this as an actual answer if I were more confident. But I happen to have recently written some NBody code of my own, if you would like to see it to get some ideas or inspiration feel free to send me a message. $\endgroup$ Oct 19, 2014 at 22:52
  • $\begingroup$ What you are saying kind of makes sense, but I am a little unfamiliar with the # & construct. I know what it does, but I don't quite know how to use it properly. Do you think you could expand a little on the second question? Also the ListPointPlot3D works, thank you. $\endgroup$
    – Karl
    Oct 20, 2014 at 13:15
  • 1
    $\begingroup$ & is shorthand notation for Function[], so I basically said Function[Orbit[#]]. Within a pure function # represents a slot for input, if you have multiple inputs you can have #1, #2 etc. So Orbit[#]& is a function with 1 input slot. Using Map, you can apply the function to a set of data. Although I later realised you don't even need to use the # and &. Map[Orbit, Range[bodies]] also generates output {Orbit[1],Orbit[2],Orbit[3]..Orbit[bodies]} $\endgroup$ Oct 20, 2014 at 14:06
  • $\begingroup$ I used your idea of the function, (which works great, so thank you) but implementing map didn't actually need to happen. I am curious though about how to use Map instead of Do? Do you think you could give an example with regards to my code? Lastly, how do I make animate and ListPlot work in this new case? Should I just use Plot since I now have a function? Or do I need to use one inside the other? $\endgroup$
    – Karl
    Oct 24, 2014 at 15:09

1 Answer 1


in response to your comment, I wrote a small example using map:

First create some data:

bodies = 4;
createBody[size_] := Module[
  {rp := RandomReal[{-1.00 size, 1.00 size}],
   rv := RandomReal[{-0.07 size, 0.07 size}]},
  {{rp, rp, rp}, {rv, rv, rv}, 1./bodies}]

data = Map[createBody, Range[bodies]]


{{{0.353151, -0.362113, -0.177183}, {-0.033645, -0.0480111, -0.0448151}, 0.25},
 {{0.377195, 0.53147, 1.16556}, {0.00865239, -0.0489291, -0.0692621}, 0.25}, 
 {{0.123244, 2.74798, 1.83201}, {0.142131, -0.0563642, -0.185028}, 0.25}, 
 {{-1.05087, 3.08883, -2.36998}, {-0.0721776, 0.187565, -0.0197556}, 0.25}}

Now we tell mathematica where it can find positions, velocities etc. I'm not sure how much this slows the program down, but it makes the code easier to read.

positions := Table[data[[n, 1]], {n, 1, Length[data]}]
velocities := Table[data[[n, 2]], {n, 1, Length[data]}]
masses := Table[data[[n, 3]], {n, 1, Length[data]}]

Now we define the force on one body:

force[body_] := G *
  Sum[If[body != others, 
     Norm[positions[[others]] - 
       positions[[body]]]^3*(positions[[others]] - positions[[body]]),
     0], {others, Length[data]}]

And map this computation over all the bodies:

Map[force, Range[bodies]]


{{-0.00403823 G, 0.0791318 G, 0.0830302 G}, {-0.0101503 G, -0.00272245 G, -0.0759115 G}, 
 {0.00265276 G, -0.058108 G, -0.0355709 G}, {0.0115358 G, -0.0183013 G, 0.0284522 G}}

I'll be gone for a few hours now, but I'll gladly expand this some more later on if you like.

  • $\begingroup$ So for the Createbody function, you essentially establish the function and then the Map function is the one that actually puts the values in for it, correct? Secondly, what is masses set to do? Would it matter if they were all just 1 for theoretical purposes? $\endgroup$
    – Karl
    Oct 28, 2014 at 2:39
  • $\begingroup$ basically yes. masses just defines where mathematica can find the mass of each body, so that I can easily and clearly refer to it in my code. You can make all the masses 1 and just leave it out if you want to, I prefer to make the total mass of the system 1, so each body has mass 1/bodies. $\endgroup$ Oct 28, 2014 at 9:24
  • $\begingroup$ I have mapped out all my orbits as points, thanks to you. My last question is about how to get them plotted. Here's what I have: Animate[ ListPointPlot3D[Orbits, Epilog -> {PointSize -> 0.05, Blue, Point[{Orbits[[i, j]], {i, 1, Particles}}]}], {j, 1, n, 1}] Orbits is currently an array of the list of positions of my 3 particles. So 3 columns, 1000 positions. $\endgroup$
    – Karl
    Nov 3, 2014 at 17:00
  • $\begingroup$ You could either use some form of Graphics3D[Map[Point,yourdata]] or just use listpointplot3d straight up I think. $\endgroup$ Nov 5, 2014 at 22:23

Not the answer you're looking for? Browse other questions tagged or ask your own question.