# Generated Numbers under a constraint?

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?

-
You can't mean to include $\mu=0$. –  bill s Feb 20 '14 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 '14 at 1:05
I suppose you've at least a range for the coordinates values ... –  belisarius Feb 20 '14 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 '14 at 1:15
@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 '14 at 3:16

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}


-
+1, very, very nice. –  ciao Feb 20 '14 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 '14 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 '14 at 3:06
@Kagaratsch done. –  Kuba Feb 21 '14 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.

-