> restart:
with(ListTools): with(LinearAlgebra): with(Groebner):
>
Q:=Matrix([[x[1,1],x[1,2],x[1,3],x[1,4]],[x[1,2],x[2,2],x[2,3],x[2,4]],[x[1,3],x[2,3],x[3,3],x[3,4]],[x[1,4],x[2,4],x[3,4],x[4,4]]]):
>
U:=Transpose(Matrix([[1,-1/2,1,1],[1,5/2,2,1/2],[0,1,-1/4,-1/4],[1,4/5,-2/5,-4/5],[1,1/4,5/4,3/4],[1,-4/3,-4/3,-1]])):
>
> for t from 1 to 6 do
>
NewU[1..4,t]:=U[1..4,t]/VectorNorm(U[1..4,t],2):
> end do:
> L:=[seq(x[1,1]+x[2,2]+x[3,3]+x[4,4]-simplify(Transpose(NewU[1..4,k]).Q.NewU[1..4,k]),k=1..6)]:
> S:=combinat[choose](4,3):
>
>
mjk:=[seq(seq(Determinant(SubMatrix(Q,r,c)),r=S),c=S)]:
> Eq:=[op(L),op(mjk)]:
>
Gb:=Basis(Eq,plex(x[1,1],x[1,2],x[1,3],x[1,4],x[2,2],x[2,3],x[2,4],x[3,3],x[3,4],x[4,4])):
> f0:=subs({x[3,4]=1,x[4,4]=y},Gb[1]):
> s := sturmseq(f0,y):
> NumberZero:=sturm(s, y, -infinity,
+infinity):
> printf(" The number of nonzero real
roots of f0 is %d ", NumberZero):
>
The number of nonzero real roots
of f0 is 0
> for j from 1 to 4 do
> for k from j to 4 do
>
NewEq:=[op(Eq),x[j,k]-1,x[3,4],x[4,4]]:
> NewGb:=Basis(NewEq,plex(x[1,1],x[1,2],x[1,3],x[1,4],x[2,2],x[2,3],x[2,4],x[3,3],x[3,4],x[4,4])):
> if
NewGb[1]=1 then
>
printf(" If x[%d,%d]=1 then the polynomial system has no roots
\n", j,k):
> end
if
> end do:
> end do:
>
If x[1,1]=1 then the polynomial system
has no roots
If x[1,2]=1 then the polynomial system
has no roots
If x[1,3]=1 then the polynomial system
has no roots
If x[1,4]=1 then the polynomial system
has no roots
If x[2,2]=1 then the polynomial system
has no roots
If x[2,3]=1 then the polynomial system
has no roots
If x[2,4]=1 then the polynomial system
has no roots
If x[3,3]=1 then the polynomial system
has no roots
If x[3,4]=1 then the polynomial system
has no roots
If x[4,4]=1 then the polynomial system
has no roots
>