1
$\begingroup$

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

enter image description here

enter image description here

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$

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Dieff
    Nov 25, 2015 at 16:03

1 Answer 1

1
$\begingroup$

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}]

enter image description here

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.

$\endgroup$
5
  • $\begingroup$ 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"? $\endgroup$
    – Dieff
    Nov 25, 2015 at 22:03
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Nov 25, 2015 at 22:49
  • $\begingroup$ 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. $\endgroup$
    – Dieff
    Nov 25, 2015 at 23:19
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Nov 26, 2015 at 1:07
  • $\begingroup$ Interesting... I'll consider complex roots and see what happens.. thanks! $\endgroup$
    – Dieff
    Dec 1, 2015 at 6:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.