Simpler way to compute gravity?

Asked Oct 13 '16 at 16:23

Active Oct 13 '16 at 16:23

Viewed 117 times

I have code for a simple n-body problem:

GravityForce[pos1_, pos2_] := (pos1 - pos2)/
   EuclideanDistance[pos1, pos2]^2;
GravityForce[p_, p_] := {0, 0, 0};
Forces[pos_List] := With[{n = Length[pos]},
   Mean /@ Table[GravityForce[pos[[i]], pos[[j]]], {i, n}, {j, n}]
   ];
n = 10;
{vel, pos} = RandomReal[{-20, 20}, {2, n, 3}];
vel += -7 Forces[pos] - 10 (GravityForce[#, sun] & /@ pos);
pos += vel

Is there some way to make this simpler?

asked Oct 13 '16 at 16:23

user5601

3,573
2
24
56

2

pos += vel??! What units are you working in here? – march Oct 13 '16 at 16:35
Why not use NDSolve in first place to solve the differential equations? – Julien Kluge Oct 13 '16 at 16:45
3

I believe your Table can be done as Outer[GravityForce, pos, pos]. I think you need to divide by r in GravityForce if you intend to direct your force along a unit vector by the way (although studying your physics problem is outside the scope of this forum..) – george2079 Oct 13 '16 at 16:56
I believe it should be EuclideanDistance[pos1, pos2]^3 not ^2... – BlacKow Oct 13 '16 at 18:05
@george2079 f[p] === Outer[g, p, p] – user5601 Oct 13 '16 at 18:35
@BlacKow you're right. – user5601 Oct 13 '16 at 18:35
@march unitless – user5601 Oct 13 '16 at 18:36
sorry @george2079 I meant to say f[p] != Outer[g, p, p] – user5601 Oct 13 '16 at 19:07
I meant just the Table , so Mean/@Outer[] – george2079 Oct 13 '16 at 19:22
@george2079 Ah I see, also Outer[g,p,p,1] works – user5601 Oct 13 '16 at 20:05
2

yes, the 1 arg is needed, sorry I forgot that. As a physics matter shouldn't it be Total not Mean ? – george2079 Oct 13 '16 at 20:11
1

NDSolve would be a better option here. I've done stuff like this (https://github.com/barrycarter/bcapps/blob/master/ASTRO/bc-trajectory.m but it's hideously ugly), and the coefficient 3/2 seems to come up naturally. I think your equations are oversimplified, and, as others note, pos2-pos1 is fine for direction, but you need to divide it by its own norm to get the unit vector. – Oct 14 '16 at 13:50

Simpler way to compute gravity?

0 Answers0