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.

Let's say I have five four dimensional vectors $p_i^\mu=(p_i^1,p_i^2,p_i^3,p_i^4)$ with $i=1,2,3,4,5$. Now, I want to fill the entries of these vectors ($4\cdot5=20$ in total) with some numbers, but such that the resulting vectors satisfy the following $9$ constraints:

$$(p_i^1)^2-(p_i^2)^2-(p_i^3)^2-(p_i^4)^2=0~~~\text{for each }i\in\{1,2,3,4,5\}$$ $$p_1^\mu+p_2^\mu+p_3^\mu+p_4^\mu+p_5^\mu=0~~~\text{for each }\mu\in\{1,2,3,4\}$$

What is a convenient way to create such constrained vectors in Mathematica? Maybe there is a built in function that does just that?

share|improve this question
1  
You can't mean to include $\mu=0$. –  bill s Feb 20 at 0:56
    
If you do not like that the labels of the four dimensions for each vector start with $\mu=0$ and go to $\mu=3$ you are welcome to relabel them to $\mu\in\{1,2,3,4\}$. Actually, let me do that above, to avoid confusion. –  Kagaratsch Feb 20 at 1:05
    
I suppose you've at least a range for the coordinates values ... –  belisarius Feb 20 at 1:10
    
Sure! Let's take each coordinate out of a box with length $-L<p<L$ and $L=$ your preferred non-zero real number. –  Kagaratsch Feb 20 at 1:15
1  
@Kagaratsch Maybe you should remove random from the title in that case. Generating truly uniformly distributed random numbers under general constraint like these seems to be far from a trivial problem, and the current title implies that this is what you want. –  Szabolcs Feb 20 at 3:16

2 Answers 2

up vote 4 down vote accepted

Update - general function, description at the end

ClearAll[LightlikeVectorsThatSumUpToZero];
LightlikeVectorsThatSumUpToZero[n_Integer: 5] := Module[
 {m, p, eq, set, sol, res, s = 4},
   While[
    m = Array[p, {n, s}];
    m[[2 ;;, 2 ;;]] = RandomReal[{-1, 1}, {n - 1, s - 1}];
    m[[2, 2]] = p[2, 2];

    (# = -Total[{##2}]) & @@@ Transpose[m][[3 ;;]];
    (# = Sqrt[{##2}.{##2}]) & @@@ m[[3 ;;]];

    eq = Join[(# == -Total[{##2}]) & @@@ Transpose[m][[ ;; 2]], 
              (#^2 == {##2}.{##2}) & @@@ m[[ ;; 2]]];
    set = Flatten[Array[p, {2, 2}], 1];
    sol = Quiet@NSolve[eq, set, Reals];
    Length[sol] == 0
    ,
    ClearAll[p, m];
    ];
   (res = m /. sol[[1]])
   ] /; n > 3

LightlikeVectorsThatSumUpToZero[7]
{{-2.71877, 2.05619, 0.176773, 1.7699}, {-2.17312, -2.00133, -0.762688, -0.368013}, 
 {1.18436, -0.346408, 0.618285, -0.948914}, {1.16035, -0.625387, 0.82161, 0.529396}, 
 {0.99203, 0.339333, 0.329197, -0.872127}, {0.393379, -0.112851, -0.262639, -0.270245},
 {1.16177, 0.69045, -0.920539, 0.160005}}

The following method is very strightforward. We are generating entries till the point there are as many unknowns as we can create quations. The we are using NSolve to get last.

With RandomReal it is highly unlikely we can get dependent vectors.

ClearAll[p, m]

(*1*) m = Array[p, {5, 4}];
      m[[ 2 ;;, 2 ;;]] = RandomReal[{-1, 1}, {4, 3}];
(*2*) m[[ 2, 2]] = p[2, 2];

      (# = -Total[{##2}]) & @@@ Transpose[m][[ 3 ;;]];

      (# = Sqrt[{##2}.{##2}]) & @@@ m[[ 3 ;;]];
(*3*)

       eq = Join[
                 (# == -Total[{##2}]) & @@@ Transpose[m][[ ;; 2]],
                 (#^2 == {##2}.{##2}) & @@@ m[[ ;; 2]]
                ];
       set = Flatten[Array[p, {2, 2}], 1]; (*{p[1, 1], p[1, 2], p[2, 1], p[2, 2]}*)

       sol = NSolve[ eq, set, Reals];

       (res = m /. sol[[ 1]]); 
(*4*)

       (#^2 - {##2}.{##2}) & @@@ res // Chop
       Total /@ Transpose[res] // Chop
{0, 0, 0, 0, 0}
{0, 0, 0, 0}

enter image description here

share|improve this answer
    
+1, very, very nice. –  rasher Feb 20 at 9:55
    
@rasher Thanks ;) It seems that sometimes system can't be solved so one may want to use Check or something and redo this. –  Kuba Feb 20 at 10:01
    
Can it do more or less than 5 vectors also? My only problem with it is that I do not understand all the symbolic shortcut gibberish... –  Kagaratsch Feb 21 at 3:06
1  
@Kagaratsch done. –  Kuba Feb 21 at 9:03

OK, so here is a really scetchy solution (up to numerical accuracy). It can compute the requested vectors not only in the case of 5 vectors, but for arbitrary number (set variable point to the desired value).

Remove["Global`*"]
point = 7;
SeedRandom;
rnd[x_] := RandomReal[{-10, 10}];
mat = {};
Do[mat = Append[
 mat, {Subscript[p0, i], sg Subscript[p1, i], sg Subscript[p2, i],
   sg Subscript[p3, i]}];, {i, 1, point}];
vec = {};
Do[vec = Append[vec, 1];, {i, 1, point}];
vec1 = Diagonal[(mat /. sg -> 1).(Transpose[mat] /. sg -> -1)];
vec2 = Transpose[mat /. sg -> 1].vec;
vecf = vec1;
Do[vecf = Append[vecf, vec2[[i]]];, {i, 1, Length[vec2]}];
res = vecf - vecf;
now = 4 point - 4 - point;
now1 = now - point + 1;
If[now1 > -1, cnt3 = point - 1;, If[now > -1, cnt3 = now, cnt3 = 0]];
now = now1;
now1 = now - point + 1;
If[now1 > -1, cnt2 = point - 1;, If[now > -1, cnt2 = now, cnt2 = 0]];
now = now1;
now1 = now - point + 1;
If[now1 > -1, cnt1 = point - 1;, If[now > -1, cnt1 = now, cnt1 = 0]];
sol = mat[[All, 1]];
Do[Subscript[p3, i] = rnd[1];, {i, 1, cnt3}];
Do[sol = Append[sol, Subscript[p3, i]], {i, cnt3 + 1, point}];
Do[Subscript[p2, i] = rnd[1];, {i, 1, cnt2}]; 
Do[sol = Append[sol, Subscript[p2, i]], {i, cnt2 + 1, point}];
Do[Subscript[p1, i] = rnd[1];, {i, 1, cnt1}];
Do[sol = Append[sol, Subscript[p1, i]], {i, cnt1 + 1, point}];
sub = FindInstance[vecf == res, sol, Reals][[1]];
Do[Subscript[p0, i] = Subscript[p0, i] /. sub; 
   Subscript[p1, i] = Subscript[p1, i] /. sub; 
   Subscript[p2, i] = Subscript[p2, i] /. sub; 
   Subscript[p3, i] = Subscript[p3, i] /. sub;, {i, 1, point}];

The variable vecf containes all the constraints. After the solution has been found all the entries of vecf are properly numerically zero.

share|improve this answer

Your Answer

 
discard

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

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