6
$\begingroup$

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?

$\endgroup$
9
  • 1
    $\begingroup$ You can't mean to include $\mu=0$. $\endgroup$
    – bill s
    Commented Feb 20, 2014 at 0:56
  • $\begingroup$ 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. $\endgroup$
    – Kagaratsch
    Commented Feb 20, 2014 at 1:05
  • $\begingroup$ I suppose you've at least a range for the coordinates values ... $\endgroup$ Commented Feb 20, 2014 at 1:10
  • $\begingroup$ Sure! Let's take each coordinate out of a box with length $-L<p<L$ and $L=$ your preferred non-zero real number. $\endgroup$
    – Kagaratsch
    Commented Feb 20, 2014 at 1:15
  • 1
    $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Commented Feb 20, 2014 at 3:16

2 Answers 2

4
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ @rasher Thanks ;) It seems that sometimes system can't be solved so one may want to use Check or something and redo this. $\endgroup$
    – Kuba
    Commented Feb 20, 2014 at 10:01
  • $\begingroup$ 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... $\endgroup$
    – Kagaratsch
    Commented Feb 21, 2014 at 3:06
  • 1
    $\begingroup$ @Kagaratsch done. $\endgroup$
    – Kuba
    Commented Feb 21, 2014 at 9:03
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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