# Plotting All Possible Points Belonging to a Group Orbit

Given that $$X = \{(x,y,z) \in \mathbb{R}^3 |\, x^2 + y^2 + z^2 - 2(xy + xz + yz) = k\}\,,$$ where $k$ is a constant. Also given that a group $G$ is represented by $$\langle g_1,g_2,g_3|\, g_1^2 = g_2^2 = g_3^2 = 1_G\rangle\,.$$ $G$ acts on $X$ such that $$g_1 \cdot (x,y,z) = (2(y+z) - x,y,z)\,,$$ $$g_2 \cdot (x,y,z) = (x,2(x+z)-y,z)$$ and $$g_3 \cdot (x,y,z) = (x,y,2(x+y)-z)\,.$$ So given $(x_0,y_0,z_0) \in X$, how do I go about plotting the orbit of the point $(x_0,y_0,z_0)$, ie $$\{g \cdot (x_0,y_0,z_0) \,|\, g \in G\}\,,$$ on a 3D graph? Would ListPointPlot3D be useful to plot all the points? How about other plotting functions?

• @Yves This problem does not appear to have anything to do with equations of motion: the term "orbit" refers to the group action, not to any actual astronomical body :-). The image of this group representation is infinite (e.g., $g_1g_2$ has infinite order), yet discrete. Thus, plotting all the points in an orbit is rather a challenge! – whuber Oct 29 '12 at 11:19
• @whuber drat! Did not read properly - I am in aerospace, so the word orbit + {x,y,z} kind of triggered celestial reflexes. Either way, it is out of my hands (and over my head) now :-) – Yves Klett Oct 29 '12 at 11:54
• Did you try something on your own already? How would you normally go about deriving the relevant points? – Yves Klett Oct 29 '12 at 11:57

Because the image of the group under this (linear) representation is infinite, we will need to limit the orbits.

### Working in the abstract group

Presuming it may eventually be of interest to depict multiple orbits, let's compute a large number of group elements once and for all. It seems efficient to do this abstractly, in terms of the given presentation, before performing the matrix multiplications, because (a) the representation appears to be faithful (it introduces no new relations among the group elements) and (b) we can eliminate many unnecessary matrix multiplications at the outset. To this end, let's create a new object--word--to represent an abstract element of any group all of whose generators have order 2.

Unprotect[NonCommutativeMultiply]; Clear[NonCommutativeMultiply];
word[g___] ** word[h___ ] := word[g, h];
word[g1___, g_] ** word[g_, h2___] := word[g1, h2];
Protect[NonCommutativeMultiply];


The first line is generic--that's how multiplication works in a free group--and the second line expresses the relations $g_i^2=1$.

NestList will create all words involving up to n products of generators, starting with the identity (as expressed by the empty word). Applying Union at each stage eliminates duplicates:

twoGroup[generators_List, n_Integer] :=
Flatten[NestList[(Flatten[Outer[NonCommutativeMultiply, #, generators]] // Union) &,
{word[]}, n]] // Union


(If you do not flatten the list, it will be partitioned into sublists corresponding to the word length.)

### Computing with a linear representation of the group

The group action can now be computed by converting these abstract words into products of matrices and performing the multiplications once and for all. Because matrix multiplication (Dot) does not know the dimension of the representation for the empty word, we need to make provision for that special case.

rep = {Subscript[g, 1] -> {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}},
Subscript[g, 2] -> {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}},
Subscript[g, 3] -> {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}},
word[] -> IdentityMatrix, word -> Dot};


Applying these replacement rules to the result of twoGroup will produce a list of matrices corresponding to all the abstract group elements it outputs (including the identity). (When the group action is not faithful, there may be duplicates among these matrices.)

### Plotting point orbits in a surface

We're all set. But before showing the orbits, let's plot the surface as a reference.

f[x_, y_, z_] := (x - y - z)^2 - 4 y z;
surface[{x0_, y0_, z0_}, cmax_] :=
ContourPlot3D[f[x, y, z], {x, -cmax, cmax}, {y, -cmax, cmax}, {z, -cmax, cmax},
Contours -> {f[x0, y0, z0]}, ContourStyle -> Opacity[0.5], Mesh -> None];


This has been formulated to draw the surface passing through a specified point $x_0$, limited within a specified cube. Now we can specify any point $x$, the maximum word length $n$, and let it fly (using, say, ListPointPlot3D):

orbit[x_List, g_List] := #.x & /@ g;

Module[{x = {-1/6, -1/3, 1/4}, n = 9, points},
points = orbit[x, twoGroup[word /@  {Subscript[g,1], Subscript[g,2], Subscript[g,3]}, n] /. rep];
Show[surface[x, Max[Abs[points]]],
ListPointPlot3D[points, PlotStyle -> Directive[Black, PointSize[0.01]]]]
] Using these functions, generalizations are now easy--plotting multiple orbits, symbolizing the points by the lengths of the words they correspond to, etc. Just be a little careful: the orbit sizes grow exponentially; there are $3(2^{n})-2$ words of length through $n$.

• Thanks very much for your solution. Looks like it's not trivial, especially to those who are not too familiar with the Mathematica software. However, the surface can't be plotted. I'd like to list down the problem, but it's going to be very long, so I'll reply with a new answer to my own question. – Markeur Oct 30 '12 at 13:10
• You forgot to define orbit and twoGroup. – rcollyer Oct 30 '12 at 13:33
• @rc Sorry; I left out a line when pasting: orbit is now defined in the text. My apologies for being so long to respond--there's no power at home or work and I don't expect any for a week or two. – whuber Oct 31 '12 at 16:29
• No apologies necessary. By some miracle, we are outside of the areas in my county where there were outages, and having lived in a more hurricane prone zone, I empathize with anyone who was hit hard. Best of luck. – rcollyer Oct 31 '12 at 16:51
• @whuber: Thanks very much for the function. It works fine. Much appreciated. Anyway, I'm sorry to hear about the recent hurricane as well as your encounter. I hope the hurricane subsides as soon as possible so that everything can return back to normal. – Markeur Oct 31 '12 at 18:01

This answer is actually to reply to whuber, who offered the solution. I'd like to reply in the comment section, however, as this reply is going to be very long and it might take more than 5 comments, I think it'd be better to reply here instead. I'd like to apologise for violating the rules.

Anyway, I've tried to plot the points as well as the surface on the software, but to no avail. This is the error when n = 1 (same goes for any other value of n):

ContourPlot3D::plln: Limiting value -Abs[orbit[{-(1/6),-(1/3),1/4},{{{1,0,0},{0,1,0},{0,0,1}},{{-1,2,2},{0,1,0},{0,0,1}},{{1,0,0},{2,-1,2},{0,0,1}},{{1,0,0},{0,1,0},{2,2,-1}}}]] in {x,-Abs[orbit[{-(1/6),-(1/3),1/4},{{{1,0,0},{0,1,0},{0,0,1}},{{-1,2,2},{0,1,0},{0,0,1}},{{1,0,0},{2,-1,2},{0,0,1}},{{1,0,0},{0,1,0},{2,2,-1}}}]],Abs[orbit[{-(1/6),-(1/3),1/4},{{{1,0,0},{0,1,0},{0,0,1}},{{-1,2,2},{0,1,0},{0,0,1}},{{1,0,0},{2,-1,2},{0,0,1}},{{1,0,0},{0,1,0},{2,2,-1}}}]]} is not a machine-sized real number. >>

ListPointPlot3D::arrayerr: orbit[{-(1/6),-(1/3),1/4},{{{1,0,0},{0,1,0},{0,0,1}},{{-1,2,2},{0,1,0},{0,0,1}},{{1,0,0},{2,-1,2},{0,0,1}},{{1,0,0},{0,1,0},{2,2,-1}}}] must be a valid array or a list of valid arrays. >>

Show::gcomb: Could not combine the graphics objects in Show[ContourPlot3D[f[x,y,z],{x,-Abs[orbit[{-(1/6),-(1/3),1/4},{{<<3>>},{<<3>>},{<<3>>},{<<3>>}}]],Abs[orbit[{-(1/6),-(1/3),1/4},{{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}}}]]},{y,-Abs[orbit[{-(1/6),-(1/3),1/4},{{<<3>>},{<<3>>},{<<3>>},{<<3>>}}]],Abs[orbit[{-(1/6),-(1/3),1/4},{{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}}}]]},{z,-Abs[orbit[{-(1/6),-(1/3),1/4},{{<<3>>},{<<3>>},{<<3>>},{<<3>>}}]],Abs[orbit[{-(1/6),-(1/3),1/4},{{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}},{{<<3>>},{<<3>>},{<<3>>}}}]]},Contours->{f[-(1/6),-(1/3),1/4]},ContourStyle->Opacity[0.5],Mesh->None],ListPointPlot3D[orbit[<<1>>],<<1>>]]. >>

And this is the output:

Show[ContourPlot3D[ f[x, y, z], {x, -Abs[ orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}]], Abs[orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}]]}, {y, -Abs[ orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}]], Abs[orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}]]}, {z, -Abs[ orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}]], Abs[orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}]]}, Contours -> {f[-(1/6), -(1/3), 1/4]}, ContourStyle -> Opacity[0.5], Mesh -> None], ListPointPlot3D[ orbit[{-(1/6), -(1/3), 1/ 4}, {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{-1, 2, 2}, {0, 1, 0}, {0, 0, 1}}, {{1, 0, 0}, {2, -1, 2}, {0, 0, 1}}, {{1, 0, 0}, {0, 1, 0}, {2, 2, -1}}}], PlotStyle -> Directive[GrayLevel, PointSize[0.01]]]]

I surmise there should be some problem with the function "orbit[x, twoGroup[word /@ {Subscript[g,1], Subscript[g,2], Subscript[g,3]}, n] /. rep]" in the "Module" function. I've tried to search for "orbit" function in the internet, but can't get anything. Is it because of the compatibility issue? I'm using Mathematica v8.0.4. And was wondering whether the "orbit" function is predefined, or how do we go about defining such a function?

• There are two functions that are missing: orbit and twoGroup. – rcollyer Oct 30 '12 at 13:33
• You are right, I had overlooked twoGroup. Incidentally, you should look at the editing help. – rcollyer Oct 30 '12 at 15:22