# Improving strategy to solve this constrained root-finding problem

I am looking for roots of the function eqs defined below. It has $2k$ unknowns and returns a vector of $2k-1$ components.

The sole purpose of the following block is to define eqs, no need to try and understand it:

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[sList_]:=Block[{tmp,T=Total[sList],k=Length[sList]},
tmp=Table[0,{i,1,k},{j,1,k}];
Table[tmp[[i,j]]=Sqrt[m[[n]]]*psi[n,T,2(Total[sList[[1;;j]]]-Total[sList[[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[sList_]:=Block[{tmp,T=Total[sList],k=Length[sList]},
tmp=Table[0,{i,1,k},{j,1,k}];
Table[tmp[[i,j]]=phi[n,T,2(Total[sList[[1;;j]]]-Total[sList[[1;;i]] ])-T],{i, 1, k}, {j, i,k}];
Table[tmp[[j,i]]=tmp[[i,j]],{i,1,k},{j,i,k}];tmp
];

eqs[vars_]:=Block[{k=Length[vars]/2},(aMat[vars[[1;;k]]].vars[[k+1;;2k]])[[1;;-2]]~Join~(bMat[vars[[1;;k]]].vars[[k+1;;2k]]-ConstantArray[1,k])]


I am looking for solutions with all positive values and FindRoot cannot deal with constraints. Of course one choice could be to use FindRoot until the $2k-1$ components of solution are positive, but I think this is not an optimal strategy.

As suggested in this SE question, I wrote the problem as a constrained optimisation problem:

subdomain = Fold[And, Table[0 < vars[[i]], {i, 2, 2 k}]];
res = FindMinimum[{Total[eqs[vars0]^2], subdomain},
Transpose[{vars[[2 ;; -1]], RandomReal[{1, 4}, 2 k - 1]}],
Method -> "InteriorPoint"]


I managed to get positive solutions for $k=6$ but it is not very robust and I would like to get solutions for larger $k$ (such as $15$, e.g.).

So my question is if someone see a way of improving my code or strategy in terms of robustness and efficiency.

A second optional question is to force the first $k$ values of the solution to be all different.

Note I have already asked a different question on the same set of equations here. Answers to the present question will probably make the previous one obsolete, so I'll delete it (or answer it if it is relevant).

• So much code often dissuades potential readers. Can you provide a reduced code that still captures the essence of your question? – bbgodfrey Nov 28 '15 at 17:42
• @bbgodfrey I cannot reduce it more. The first block only defines the equations (eqs) and should not be understood in details. I'll edit my question to make it clear, and remove the part on FindRoot to lighten it. – anderstood Nov 28 '15 at 17:49

One possibility is to add an equation that the sum of variables is some constant, say 1. You can actually enforce nonnegativity by augmenting the starting values, per refguide page for FindRoot. I won't repeat the lengthy setup code but just show the steps from there.

k = 9;
len = 2*k;
vars = Array[s, len];
eqns = Join[eqs[vars], {Total[vars] - 1}];
start = Thread[{vars, RandomReal[2, len]/len, 0, 1}];
FindRoot[eqns == 0, start]

During evaluation of In:= FindRoot::reged: The point {0.03140964551874968,0.03874740609445434,0.0743937162451151,0.06308530530957773,0.1144968686002389,0.1358768793893159,0.00637481057165916,<<5>>,0.03178764934909824,0.1038665079060305,0.03872218273437743,0.01836388193562395,0.01507199254295023,0.} is at the edge of the search region {0.,1.} in coordinate 18 and the computed search direction points outside the region. >>

(* Out= {s -> 0.03140964551874968, s -> 0.03874740609445434,
s -> 0.0743937162451151, s -> 0.06308530530957773,
s -> 0.1144968686002389, s -> 0.1358768793893159,
s -> 0.00637481057165916, s -> 0.01906482002753885,
s -> 0.03327004567338578, s -> 0.01749359129816518,
s -> 0.0957320547892758, s -> 0.05456770646648509,
s -> 0.03178764934909824, s -> 0.1038665079060305,
s -> 0.03872218273437743, s -> 0.01836388193562395,
s -> 0.01507199254295023, s -> 0.} *)


The message is not encouraging. My guess is it will be quite difficult to meet the constraints.