Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Im Trying to model a cloud of point masses that act according to gravity, what im strugeling with is the exclusion of cases where euclidian distance = 0, If I try to do it with an If statement in my forces function, NDsolve complains about too many dependant variables I think, but I think its rather a fault in my usage/understanding of NDsolve

Any help is appreciated

edit: im aware that the equations arent yet in working order

n = 20;
g = 1;
pos = RandomReal[{-10, 10}, {n, 3}];
masses = RandomReal[{1, 50}, n];
Tmax = 20;

Forces[xp_, yp_, zp_] := 
   EuclideanDistance[{xp, yp, zp}, {x[i][t], y[i][t], z[i][t]}]^2
  , {i, 1, n}]

  Sphere[pos[[i]], ((masses[[i]]*4)/(3*Pi))^(1/3)/5]
  , {i, 1, n}]
 , PlotRange -> 11.5]

sol = NDSolve[
     (x[i])''[t] == Forces[x[i][t], y[i][t], z[i][t]][[1]]/masses[[i]],
     (y[i])''[t] == Forces[x[i][t], y[i][t], z[i][t]][[2]]/masses[[i]],
     (z[i])''[t] == Forces[x[i][t], y[i][t], z[i][t]][[3]]/masses[[i]],
     (x[i])'[0] == 0,
     (y[i])'[0] == 0,
     (z[i])'[0] == 0,
     x[i][0] == pos[[i, 1]],
     y[i][0] == pos[[i, 2]],
     z[i][0] == pos[[i, 3]]
    , {i, 1, n}]
     x[i][t], y[i][t], z[i][t]
    , {i, 1, n}]
  , {t, 0, Tmax}]
share|improve this question
You might be interested in this notebook, which contains code to deal with a fairly large number of objects specified in {mass,pos,vel} form. When two objects collide, however, a singularity arises and you simply cannot expect to integrate past that point. You can use InterpolatingFunctionDomain to deal with it gracefully but that's about it. –  Mark McClure Jun 23 '14 at 0:50

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.