Matrix operations in "Inertia Versus Gravity in 3d" Demonstration

Questions
, ,
1

I have been deconstructing Michael Schreiber’s Inertia Versus Gravity In 3D author.nb from the Demo project. The main work is done in the inertiaversusgravity3D function which I have “exploded” below.

inertiaversusgravity3D[m1_, m2_] := Module[
   {tt}, {-m2 + 2 m1 + (
      Plus @@ # & /@ Map[
        With[
          {tt = (#[[1]]^2 + #[[2]]^2 + #[[3]]^2)},
          If[tt > 0., #/tt, {0., 0., 0.}]
          ]
         &,
        (1 - IdentityMatrix[Length[m1]]) Outer[Plus, -m1, m1,1], {2}]) , m1}
   ];


{m = N[pts RandomReal[{-1, 1}, {pts, 3}]]}
(* The first call would be *)
{m1,m2}=inertiaversusgravity3D[m, m]

My question is in the line:
(1 - IdentityMatrix[Length]) Outer[Plus, -m1, m1,1], {2}])

This looks like it is going to multiply two matricies but there is not a default multiplication operation either the dot or cross product must be specified. When I , define m1 manually and look at the results of this line I get Outer[Plus, -m1, m1,1], {2}] unchanged.

InertiaVsGravityMatrix.png
InertiaVsGravityMatrix.png567×482 11.4 KB

My question is should the matricies be dotted together and is the operation Matrix1 Matrix2 (Matrix1 * Matrix2) defined? If so what is the definition.

Any insight would be appreciated.

Mark

2

If m1 is a vector of length n, the Outer makes it an nxn matrix. The multiplication is coordinate-wise. The matrix on the left is in effect a mask that makes diagonal elements of the one on the right zero and leaves others alone.

3

Doesn’t the Outer[Plus,-m1,m1,1] guarntee that the diagonal will be 0?

Thanks,
Mark

4

One would think so, yes. So I’m as confused as anyone as to why that multiplication is needed at all.

5

Again, Thanks for your time and for the validation of my understanding of the code involved. Mathematica is extremely powerful but understanding its code is difficult for me. Kind of like perl, the sometimes turse syntax can be almost imposible to understand on the 1st, 2nd, 3rd … reading.

I am left with the question of how this models inertia or gravity. The whole function seems to be based on the sums of “distances” of all other points on a target point. I am assuming that the square root of the sum of the squares was not taken as a performance boost. But there does not seem to be any modeling of the effects being lessened as the distance increases nor of a points resistance to change (inertia.) This could be typical in this kind of a model but it seems to defeat the point of modeling.

It is neat to watch though.

Thanks
Mark

6

Just FYI - some related information. In the documentation there is an interesting demo on the N-body simulation is a classic Newtonian problem. I think it was built on the following tutorial which may give you some additional information.

Needs["OpenCLLink`"]

srcf = FileNameJoin[{$OpenCLLinkPath, "SupportFiles", "NBody.cl"}];

NBody = OpenCLFunctionLoad[{srcf}, 
  "nbody_sim", {{"Float[4]", _, "Input"}, {"Float[4]", _, "Input"}, _Integer, 
   "Float", "Float", {"Local", "Float"}, {"Float[4]", _, 
    "Output"}, {"Float[4]", _, "Output"}}, 256];

numParticles = 1024;
deltaT = 0.05;
epsSqrt = 50.0;


pos = OpenCLMemoryLoad[RandomReal[512, {numParticles, 4}], "Float[4]"];
vel = OpenCLMemoryLoad[RandomReal[1, {numParticles, 4}], "Float[4]"];
newPos = OpenCLMemoryAllocate["Float[4]", {numParticles}];
newVel = OpenCLMemoryAllocate["Float[4]", {numParticles}];

NBody[pos, vel, numParticles, deltaT, epsSqrt, 256*4, newPos, newVel, 1024];
NBody[newPos, newVel, numParticles, deltaT, epsSqrt, 256*4, pos, vel, 1024];

Graphics3D[Point[
  Dynamic[Refresh[
    NBody[pos, vel, numParticles, deltaT, epsSqrt, 256*4, newPos, newVel, 
     1024];
    NBody[newPos, newVel, numParticles, deltaT, epsSqrt, 256*4, pos, vel, 
     1024];
    Take[#, 3] & /@ OpenCLMemoryGet[pos], UpdateInterval -> 0]]]]

232l3kf0dgfhioweLKJ987fs9d.gif