2
$\begingroup$

Has someone an idea how to solve the 1. question for ANY number of mass points and a clue to the 2. question?

center of mass system, momentum

$\endgroup$
  • 1
    $\begingroup$ Didn't you tell me to do it like this in main question?? $\endgroup$ – Tom Jan 14 at 9:24
  • $\begingroup$ Yes, but I didn't mean something like the above which just references the previous question. I meant a fully written out question with the full problem statement so the new question can stand alone. $\endgroup$ – m_goldberg Jan 14 at 9:34
  • $\begingroup$ @Tom Indeed, a bit more self-containedness of this question would be desirable. When I read it the first time I thought "What the ... do they mean with '1. question'?" $\endgroup$ – Henrik Schumacher Jan 14 at 9:38
  • $\begingroup$ ok, in future I will try to do it more correctly! thank you for patience ! $\endgroup$ – Tom Jan 14 at 9:42
  • $\begingroup$ Near duplicate: mathematica.stackexchange.com/questions/189362/… $\endgroup$ – David G. Stork Jan 14 at 18:43
2
$\begingroup$

The data for many particles are best stored in arrays.

n = 10000;
masses = RandomReal[{0.1, 1.}, {n}];
positions = RandomReal[{-1, 1}, {n, 3}];
velocities = RandomReal[{-1, 1}, {n, 3}];

Now, one obtains

centerofmass = masses.positions;
momenta = masses velocities;
totalmomentum = masses.velocities;
angularmomenta = MapThread[Cross, {positions, momenta}];
totalangularmomentum = Total[angularmomenta];

Performance Tuning

There is an issue with Cross: For some reason, it is way slower than it should be. Here is a compiled, listable, and parallelized version of it for 3-vectors:

cross = With[{
   code = Cross[
     Table[Compile`GetElement[X, i], {i, 1, 3}],
     Table[Compile`GetElement[Y, i], {i, 1, 3}]
     ]
   },
  Compile[{{X, _Real, 1}, {Y, _Real, 1}},
   code,
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ]
  ]

It is about 500 times faster than Cross on my machine and it is already Listable so that we don't need MapThread:

angularmomenta2 = cross[positions, momenta];

Up to machine precision, the result is the same:

Max[Abs[angularmomenta - angularmomenta2]]

2.22045*10^-16

$\endgroup$
  • $\begingroup$ you´re genious Henrik! thank you! $\endgroup$ – Tom Jan 14 at 9:43
  • $\begingroup$ Huh, thanks for the flowers. You're welcome. $\endgroup$ – Henrik Schumacher Jan 14 at 9:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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