Help with findroot optimization

Question

I'm trying to solve the set of coupled equations

$$\frac{-N -2( \lambda + N(\frac{\beta}{\epsilon}-\lambda))\upsilon_l + N \upsilon_l^2-2(N-1)\gamma\upsilon_l^3}{\gamma-2\lambda\upsilon_l^2+\gamma\upsilon_l^4}=\sum_{j\neq l}^N \frac{2}{\upsilon_j-\upsilon_l}\,\,\,\,\,(1)$$

With the solutions $\upsilon_l$ I can compute the quantity

$$E=\frac{N \epsilon \lambda}{2}+\frac{\beta N}{2}-\frac{\epsilon}{2}\sum_{j=1}^N\upsilon_j-\frac{\epsilon \gamma}{N}\sum_{j=1}^{N-1}\sum_{l=j+1}^N\upsilon_j\upsilon_l\,\,\,\,\,(2)$$

Since the solutions of equation (1) are unique, each set $\{\upsilon_j\}$ gives me just one value to $E$ in (2), and I can check this value by another methods. I want to set $N=100,\,\lambda=5,\,\beta=0,\,\epsilon=1$ and vary $\gamma$ in the interval $[0,1]$. and find the solutions using FindRoot. The solutions are very sensitive to the initial conditions, so I wrote a simple routine that uses the previous solutions as starting points to find the new ones. The program is the following:

Equations

ν[Npart_] := Table[Symbol["ν" <> ToString[i]], {i, 1, Npart}];
u[Npart_, j_] := ν[Npart][[j]];
BAE[Npart_, λ_, β_, ϵ_, γ_] := 
  SetPrecision[
   ParallelTable[(-Npart - 
      2 ( λ + Npart (β/ϵ - λ)) u[Npart,
         l] + Npart  u[Npart, l]^2 - 
      2 (-1 + Npart) γ  u[Npart, l]^3)/(γ - 
      2 λ u[Npart, l]^2 + γ u[Npart, l]^4) == 
     Sum[2/(u[Npart, j] - u[Npart, l]), {j, 
       Complement[Range[Npart], {l}]}], {l, 1, Npart}], 200];

Inicial points

BAEroots[Npart_, λ_, β_, ϵ_, γ_] := 
 FindRoot[BAE[Npart, λ, β, ϵ, γ], 
  Table[{u[Npart, j], 0.01 + (j - 1)/Npart}, {j, 1, Npart}], 
  WorkingPrecision -> 200, AccuracyGoal -> 100, PrecisionGoal -> 100, 
  MaxIterations -> 1000]

Routine

NewRoots[Npart_, λ_, β_, ϵ_, γfinal_] :=
  Block[{SP = BAEroots[Npart, λ, β, ϵ, 0]},
  Do[NSP = 
    FindRoot[BAE[Npart, λ, β, ϵ, i], 
     Table[{u[Npart, j], u[Npart, j] /. SP}, {j, 1, Npart}], 
     WorkingPrecision -> 200, AccuracyGoal -> 100, 
     PrecisionGoal -> 100, MaxIterations -> 1000];
   Print["*****************************************"];
   Print["γ = ", N[i]];
   Print["roots : ", SetPrecision[NSP, 20]];
   Print["Energy : ", 
    SetPrecision[((Npart ϵ λ)/2 + (β Npart)/
        2 + -ϵ/
         2 Sum[u[Npart, j], {j, 1, Npart}] - (ϵ i)/
         Npart Sum[
          u[Npart, j] u[Npart, l], {j, 1, Npart - 1}, {l, j + 1, 
           Npart}]) /. NSP, 20]];
   SP = NSP, {i, 0, γfinal, γfinal/100}]]

The routine run well, but at some point an error occour

I know from theory that all roots must be real and positive. Have you any ideia about what's wrong? Any help will be much appreciable.

Regards, Dieff

P.S.: In the program, $Npart=N=100$

Answer 1

For the code in the question, ν1 is not always positive, as can be seen by plotting it as a function of γ.

ListPlot[%[[2]], DataRange -> {0, .1}, AxesLabel -> {γ, ν1}]

where [%[[2]] is obtained from Sow[ν1 /. NSP] within NewRoots. Evidently, the computation breaks down when ν1 decreases to zero. I have run similar computations for smaller Npart, with similar results, although at more negative values of γ. In the extreme case of Npart = 3, values of ν can be obtained with Solve in a reasonable amount of time.

Solve[BAE[3, 5, 0, 1, -149/100] // Rationalize, {ν1, ν2, ν3}, Reals] // N;
Union[Sort[#] & /@ ({ν1, ν2, ν3} /. %), SameTest -> (Norm[#1 - #2] < 10^-5 &)]
(* {{-1., 0.249264, 4.01181}, {0.000769822, 1., 1299.}} *)

Solve[BAE[3, 5, 0, 1, -150/100] // Rationalize, {ν1, ν2, ν3}, Reals] // N;
Union[Sort[#] & /@ ({ν1, ν2, ν3} /. %), SameTest -> (Norm[#1 - #2] < 10^-5 &)]
(* {{-1., 0.249668, 4.00533}} *)

Evidently, the set {ν1, ν2, ν3} with all elements positive ceases to exist between γ = -1.49 and γ = -1.50. Probably the same happens for Npart = 100, and the negative ν1 values appearing in the plot above are spurious, associated with the error reported in the question.