# Help with findroot optimization

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$

• Npart=100, \lambda=5, \beta=0, \epsilon=1, \gamma final = -1. The program starts finding the solution with \gamma=0 and then decrease the value of \gamma in steps of -0.01. From \gamma = 0 to \gamma = -0.08 the program ran well. – Dieff Nov 25 '15 at 16:03

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

ListPlot[%[], DataRange -> {0, .1}, AxesLabel -> {γ, ν1}] where [%[] 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.

• Thanks for your answer. Actually it was expected that the values of v1 go to zero as γ become more negative. But even very, very small, the value of v1 must be positive. For example, when I run the routine in the interval -0.09 < γ < - 0.089 the program works well until γ = -0.0899, with v1 = 0.00004453031... I can't understand why the program fail in the next step. What did you mean by "ceases to exist"? – Dieff Nov 25 '15 at 22:03
• @Dieff My answer shows that ν1 > 0 is not true for Npart = 3. How large must Npart, be before it is true? My plot suggests, larger than 100. – bbgodfrey Nov 25 '15 at 22:49
• In principle that should be true for any value of N. Note that the solutions of the equation (1) satisfy the condition v(N-j+1) = 1/vj when \beta=0. If N is an odd number, one root must be 1 because v(N/2+1) = 1/v(N/2+1). So it's kinda weird that mathematica found a negative root in your example. – Dieff Nov 25 '15 at 23:19
• @Dieff Note that ν1 approximately equal to 1/νN does not guarantee that ν1 is a positive real, because νN could be complex. This is exactly what happens for Npart = 3 at γ >= 3/2. – bbgodfrey Nov 26 '15 at 1:07
• Interesting... I'll consider complex roots and see what happens.. thanks! – Dieff Dec 1 '15 at 6:51