3
$\begingroup$

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
$\endgroup$
  • $\begingroup$ It seems to me your question is fundamentally mathematical more than it is Mathematica-specific, so you may have better luck asking on the maths StackExchange. $\endgroup$ – Rahul Nov 26 '15 at 4:41
  • $\begingroup$ How do you have a square matrix A(t) that is both of full rank and has a determinant of zero? There is something fundamentally messed up here at least for the description of the n-even case. $\endgroup$ – Daniel Lichtblau Nov 28 '15 at 17:32
  • $\begingroup$ When $n$ is even, $A(t)$ is generally of full rank. But the solutions of my equations correspond to $A(t)$ which is not of full rank, because the equations include ${\rm Pf}(A)=0$. Is it more clear? $\endgroup$ – anderstood Nov 28 '15 at 17:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.