Group Abstract Group Abstract

Message Boards Message Boards

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

GROUPS:
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[m1]]) 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.


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 
POSTED BY: Mark Henwood
Answer
1 year ago
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.
POSTED BY: Daniel Lichtblau
Answer
1 year ago
Doesn't the Outer[Plus,-m1,m1,1] guarntee that the diagonal will be 0?

Thanks,
Mark  
POSTED BY: Mark Henwood
Answer
1 year ago
One would think so, yes. So I'm as confused as anyone as to why that multiplication is needed at all.
POSTED BY: Daniel Lichtblau
Answer
1 year ago
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
POSTED BY: Mark Henwood
Answer
1 year ago
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]]]]

POSTED BY: Sam Carrettie
Answer
1 year ago