Edit Interestingly, I realized solving my problem with $2n$ equations appeared to be more efficient than using mathematical results on the solutions to reduce the problem to $n$ equations (the $n$ equations given here). I have asked a different question related to the set of $2n$ equations here. Depending on the outcome, I will delete or answer this question.
There are many questions about solving systems of nonlinear equations. I think this one has its own specificity.
The mathematical problem is to find two vectors $x$ and $t$ in $\mathbb{R}^n$ such that:
$$A(t)x=0 \quad\text{and}\quad B(t)x=j$$
where $A$ is a given skew-symmetric $n\times n$ matrix, $B$ a given symmetric $n\times n$ matrix and $j^\top=[1,\dots,1]$. $A,B$ are given in the code at the end of the post (aMat is $A$ and bMat is $B$).
Odd case (OK) When $n=2p+1$, $x$ is in the kernel of $A(t)$ and can be explicitly expressed in terms of the elements of $A(t)$ (see this mathoverflow question). Let's write this as $x=\phi(A(t))$. So $x$ can be eliminated and the equations in $t$ can be written as $B(t)\phi(A(t))=j$ which gives $n-1$ equations ($n$ equations, but one is redundant). This works very smoothly and I managed do what I wanted in this case.
Even case (not OK) When $n=2p$, $A(t)$ is in general of full rank. A first approach is to consider that, in general, $B(t)$ is invertible, so $x$ can be eliminated by solving for $t$ the equations: $$A(t)B(t)^{-1}j=0$$
Necessarily, for this system to have non trivial solutions, $\det(A(t))=0$. Actually, $A$ is skew-symmetric hence it is less costly to compute ${\rm Pf}(A(t))=0$ (see Pfaffian). Then, it follows from its definition that $A(t)$ is of rank $n-2$. So in total I have $k-1$ equations, again.
This is what I implemented as
eqs[t_] := ({Pf[aMat[t]]}~Join~
((aMat[t]).Inverse[bMat[t]].ConstantArray[1, Length[t]])[[1 ;;
Length[t] - 2]])
It "works", but it is a bit slow:
FindRoot[eqs[{t1, t2, t3, 6}],
Transpose[{{t1, t2, t3}, RandomReal[{.5, 5}, 3]}]]
takes about 8 seconds on my computer, and corresponds to $n=4$ (I'd like to increase $n$ to 8 or more, if possible).
My question is about solving the system when $n$ is even : is there a better way of solving it?
Additionally, I have a result on the kernel of $A(t)$ which I am not using in the previous equation. Since $A(t)$ is skew-symmetric, when ${\rm Pf}A(t)=0$, the columns of its copfaffian are of rank 2. The copfaffian can be computed in Mathematica by:
copf[t_]:=Block[{com=Table[If[i>j,0,(-1)^(i+j)*Pf@Transpose@Delete[
Transpose@Delete[aMat[t],{{i},{j}}],{{i},{j}}]],{i,1,4},{j,1,4}]},
com-Transpose[com]]
I tried using it to speed up the resolution but I was not successful, see function eqs2 below.
Full MWE:
(********* GENERAL DEFINITIONS *********)
phi[p_, T_, t_] := -1/Sqrt[m[[p]]]*Sum[Cos[omegas[[j]]*t/2]/(omegas[[j]]*Sin[omegas[[j]]*T/2]) q[[p, j]]*q[[n,j]], {j, 1, n}];
psi[p_, sigma_, t_] := 1/Sqrt[m[[p]]]* Sum[Sin[omegas[[j]]*t/2]/(Sin[omegas[[j]]*sigma/2]) q[[p, j]]* q[[n,j]], {j, 1, n}];
n=2;
m={.5,.5};
q={{-0.8506508083520401`,-0.5257311121191336`},{0.5257311121191336`,-0.8506508083520401`}};
omegas={3.23606797749979`,1.2360679774997898`};
aMat[t_]:=Block[{tmp,T=Total[t],k=Length[t]},
tmp=Table[0,{i,1,k},{j,1,k}];
Table[tmp[[i,j]]=Sqrt[m[[n]]]*psi[n,T,2(Total[t[[1;;j]]]-Total[t[[1;;i]] ])-T],{i, 1, k}, {j, i+1,k}];
Table[tmp[[j,i]]=-tmp[[i,j]],{i,1,k},{j,i,k}];tmp
];
bMat[t_]:=Block[{tmp,T=Total[t],k=Length[t]},
tmp=Table[0,{i,1,k},{j,1,k}];
Table[tmp[[i,j]]=phi[n,T,2(Total[t[[1;;j]]]-Total[t[[1;;i]] ])-T],{i, 1, k}, {j, i,k}];
Table[tmp[[j,i]]=tmp[[i,j]],{i,1,k},{j,i,k}];tmp
];
Pf[A_]:=If[Length[A]==0,1,Module[{L,A1,MatrixDelete},MatrixDelete[M_,i_]:=Delete[#,i]&/@Delete[M,i];
L=Length[A];A1=MatrixDelete[A,1];
Sum[(-1)^i (A[[1]][[i]] Pf[MatrixDelete[A1,i-1]]),{i,2,L}]]]
eqs[t_]:=({Pf[aMat[t]]}~Join~((aMat[t]).Inverse[bMat[t]].ConstantArray[1,Length[t]])[[1;;Length[t]-2]])
(************* END OF GENERAL DEFINITIONS *********)
(* works, but slow *)
FindRoot[eqs[{t1,t2,t3,6}],Transpose[{{t1,t2,t3},RandomReal[{.5,5},3]}]]//AbsoluteTiming
(* output: {7.558096`,{t1->3.8374275530641726`,t2->3.9630506992218164`,t3->3.8374275530641713`}} *)
(* uses info of kernel but does not converge *)
copf[t_]:=Block[{com=Table[If[i>j,0,(-1)^(i+j)*Pf@Transpose@Delete[Transpose@Delete[aMat[t],{{i},{j}}],{{i},{j}}]],{i,1,4},{j,1,4}]},com-Transpose[com]]
eqs2[t_]:=Block[{aa=aMat[t],CP=copf[t],alpha,beta},{alpha,beta}={alpha,beta}/.Flatten@Solve[(bMat[t].(alpha*CP[[1]]+beta*CP[[2]]))[[{1,2}]]=={1,1}];
{Pf[aa]}~Join~((bMat[t].(alpha*CP[[1]]+beta*CP[[2]]))[[3;;-1]]-ConstantArray[1,Length[t]-2] )]
eqs2[{t1,t2,t3,6}/.{t1->3.8374275530641726`,t2->3.9630506992218164`,t3->3.8374275530641713`}]//Chop;
(* output: -> {0,0,0}: this means a solution of `eqs` is apparently also a solution of `eqs2` *)
(* NEVER CONVERGES! *)
FindRoot[eqs2[{t1,t2,t3,6}],Transpose[{{t1,t2,t3},RandomReal[{.5,5},3]}]]//AbsoluteTiming
