Question: i can't find our where my error is...

Hi,

Sorry to ask such a stupid question but I can't find out where my error is. Probably it's so huge it blinds me!

The double loop and the matrix product F^+ . F should give the same result, no? (it seems that F^+ . F has its rows reordered ?)


 

restart:

N   := 3:
P   := 2:
niv := [seq(Z[i], i=1..N)];
f   := Matrix(N^P, P, (i,j) -> `if`(j=P, niv[(i mod 3)+1], niv[iquo(i-1,3)+1]));

niv := [Z[1], Z[2], Z[3]]

 

f := Matrix(9, 2, {(1, 1) = Z[1], (1, 2) = Z[2], (2, 1) = Z[1], (2, 2) = Z[3], (3, 1) = Z[1], (3, 2) = Z[1], (4, 1) = Z[2], (4, 2) = Z[2], (5, 1) = Z[2], (5, 2) = Z[3], (6, 1) = Z[2], (6, 2) = Z[1], (7, 1) = Z[3], (7, 2) = Z[2], (8, 1) = Z[3], (8, 2) = Z[3], (9, 1) = Z[3], (9, 2) = Z[1]})

(1)

ds := subs(niv =~ [$0..N-1], f);

ds := Matrix(9, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 0, (2, 2) = 2, (3, 1) = 0, (3, 2) = 0, (4, 1) = 1, (4, 2) = 1, (5, 1) = 1, (5, 2) = 2, (6, 1) = 1, (6, 2) = 0, (7, 1) = 2, (7, 2) = 1, (8, 1) = 2, (8, 2) = 2, (9, 1) = 2, (9, 2) = 0})

(2)

vs := [ seq(V__||i, i=1..P)]:
es := unapply( sort( [ seq( mul(vs ^~ [entries(ds[i,..], nolist)]), i=1..N^P) ] ), vs);
 

proc (V__1, V__2) options operator, arrow; [1, V__1, V__2, V__1^2, V__2^2, V__1*V__2, V__1*V__2^2, V__1^2*V__2, V__1^2*V__2^2] end proc

(3)

ff := convert([ seq(es(entries(ffd[i,..], nolist)), i=1..N^P) ], Matrix);


UnityRoots := [solve(z^3=1, z)]:
F := simplify(subs(niv =~ UnityRoots, ff)) /~ sqrt(N^P):

ff := Matrix(9, 9, {(1, 1) = 1, (1, 2) = Z[1], (1, 3) = Z[2], (1, 4) = Z[1]^2, (1, 5) = Z[2]^2, (1, 6) = Z[1]*Z[2], (1, 7) = Z[1]*Z[2]^2, (1, 8) = Z[1]^2*Z[2], (1, 9) = Z[1]^2*Z[2]^2, (2, 1) = 1, (2, 2) = Z[1], (2, 3) = Z[3], (2, 4) = Z[1]^2, (2, 5) = Z[3]^2, (2, 6) = Z[1]*Z[3], (2, 7) = Z[1]*Z[3]^2, (2, 8) = Z[1]^2*Z[3], (2, 9) = Z[1]^2*Z[3]^2, (3, 1) = 1, (3, 2) = Z[1], (3, 3) = Z[1], (3, 4) = Z[1]^2, (3, 5) = Z[1]^2, (3, 6) = Z[1]^2, (3, 7) = Z[1]^3, (3, 8) = Z[1]^3, (3, 9) = Z[1]^4, (4, 1) = 1, (4, 2) = Z[2], (4, 3) = Z[2], (4, 4) = Z[2]^2, (4, 5) = Z[2]^2, (4, 6) = Z[2]^2, (4, 7) = Z[2]^3, (4, 8) = Z[2]^3, (4, 9) = Z[2]^4, (5, 1) = 1, (5, 2) = Z[2], (5, 3) = Z[3], (5, 4) = Z[2]^2, (5, 5) = Z[3]^2, (5, 6) = Z[2]*Z[3], (5, 7) = Z[2]*Z[3]^2, (5, 8) = Z[2]^2*Z[3], (5, 9) = Z[2]^2*Z[3]^2, (6, 1) = 1, (6, 2) = Z[2], (6, 3) = Z[1], (6, 4) = Z[2]^2, (6, 5) = Z[1]^2, (6, 6) = Z[1]*Z[2], (6, 7) = Z[1]^2*Z[2], (6, 8) = Z[1]*Z[2]^2, (6, 9) = Z[1]^2*Z[2]^2, (7, 1) = 1, (7, 2) = Z[3], (7, 3) = Z[2], (7, 4) = Z[3]^2, (7, 5) = Z[2]^2, (7, 6) = Z[2]*Z[3], (7, 7) = Z[2]^2*Z[3], (7, 8) = Z[2]*Z[3]^2, (7, 9) = Z[2]^2*Z[3]^2, (8, 1) = 1, (8, 2) = Z[3], (8, 3) = Z[3], (8, 4) = Z[3]^2, (8, 5) = Z[3]^2, (8, 6) = Z[3]^2, (8, 7) = Z[3]^3, (8, 8) = Z[3]^3, (8, 9) = Z[3]^4, (9, 1) = 1, (9, 2) = Z[3], (9, 3) = Z[1], (9, 4) = Z[3]^2, (9, 5) = Z[1]^2, (9, 6) = Z[1]*Z[3], (9, 7) = Z[1]^2*Z[3], (9, 8) = Z[1]*Z[3]^2, (9, 9) = Z[1]^2*Z[3]^2})

(4)

Scalar products of pairs of comumn vectors

F must be an orthogonal array

for i1 from 1 to N^P do
  for i2 from 1 to N^P do
    printf("%a ", simplify(add(F[..,i1] . F[.., i2])))
  end do:
  printf("\n"):
end do:
printf("\n");

1 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0
0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 1
 

 

or more simply:

simplify(F^+ . F)

Matrix([[1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0]])

(5)

 


 

Download Too_Blind_To_Find_My_Mistake.mw

Please Wait...