# Inconsistent results solving system of equations with FindRoot

I have been trying to solve the following problem (Is basically a physics problem: an expansion of equations of state up to 1st order of the order parameter in order to solve for the critical temperature, or, to be more precise, the inverse temperature $$\beta$$ for our case, represented in the code by the variable beta).

So basically the system of equations that I want to solve (I put directly the Mathematica code) is

Delta = 3/10
phi = 1 + Delta^2/12
FindRoot[{Etn[Delta, 1, beta, l1, l2, phi, 0]^2 == beta^2(Etn[Delta, 1, beta, l1, l2, phi, 2]^2 - Etn[Delta, 1, beta, l1, l2, phi, 4]*Etn[Delta, 1, beta, l1, l2, phi, 0]),
phi*Etn[Delta, 1, beta, l1, l2, phi, 0] == Etn[Delta, 1, beta, l1, l2, phi, 2],
Etn[Delta, 1, beta, l1, l2, phi, 0] == Etn[Delta, 1, beta, l1, l2, phi, 1]},
{{l1,1}, {l2,1}, {beta,2.4}}, PrecisionGoal->5, WorkingPrecision->6]


Where I have defined the functions

f[t_, phi_, l1_, l2_, beta_, J_, Delta_]:= Exp[1/2 beta^2 J^2 (phi + 2 l2/(beta J^2))t^2 + beta l1 t - beta (l1 + l2 + l2 Delta^2/12)]

Etn[Delta_, J_, beta_, l1_, l2_, phi_, n_] := (1/Delta*Integrate[t^n f[t, phi, l1, l2, beta, J, Delta], {t, 1-Delta/2, 1+Delta/2}])


Here I inserted as a seed value for the beta 2.4, this yields in the following result

{l1 -> -5467.98, l2 -> 2735.66, beta -> 3.02056*10^-6}

With the following output message

FindRoot: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than 6. digits of working precision to meet these tolerances

The problem is, when I change the seed value of beta from 2.4 to any other number, the result changes, sometimes yielding negative numbers, other numbers with an order of magnitude away from other values arrived at with other seed values, sometimes no result at all because some singularity in the Jacobian, etc...

I want to know how to obtain consistent results as beta should be a fixed and positive quantity. Moreover I want to find this beta for several values of Delta (that in the code I'm fixing to 3/10)

I'll be extremely grateful if someone can help me with this!

• [Phi] is not defined. Commented Oct 20, 2022 at 18:51
• The value of $\phi$ (see Daniel Huber's comment) is still not defined! Commented Oct 21, 2022 at 6:37
• Sorry it was a mistake transcribing the system of equations, [Phi] was simply what I have defined as phi and I set equal to 1 + Delta^2/12 Commented Oct 21, 2022 at 14:41
• I do not see a Mathematica question here yet. You may have strong reasons to believe that there is a solution somewhere near the initial values that you provide, but you do not give enough details for people here to come to this conclusion themselves. Are you expecting people to just take your word? Maybe it is not a Mathematica problem, and your system simply does not have a solution... Commented Oct 21, 2022 at 15:23
• It is easy to make examples without solutions that yield other errors, e.g. FindRoot[{1+Sqrt[x]+x==0},{x,1}]. You say it yields different "solutions": Do you mean solutions without error message, or some numbers with an error message (I would not call that a solution)? Commented Oct 21, 2022 at 16:31

This is just an extended comment that suggests what you might want to look at to determine the cause of the instability of the result.

Consider evaluating all of the pieces separately as there is an exact result for the integration:

eq1Left = Etn[Delta, 1, beta, l1, l2, phi, 0]^2;
eq1Right = beta^2 (Etn[Delta, 1, beta, l1, l2, phi, 2]^2 -
Etn[Delta, 1, beta, l1, l2, phi, 4]*Etn[Delta, 1, beta, l1, l2, phi, 0]);
eq2Left = phi*Etn[Delta, 1, beta, l1, l2, phi, 0];
eq2Right = Etn[Delta, 1, beta, l1, l2, phi, 2];
eq3Left = Etn[Delta, 1, beta, l1, l2, phi, 0];
eq3Right = Etn[Delta, 1, beta, l1, l2, phi, 1];


Now evaluate each piece with respect to the answer you got (after rationalizing the results):

sol = Rationalize[{l1 -> -5467.98, l2 -> 2735.66, beta -> 3.02056*10^-6}, 0];
N[{eq1Left, eq1Right, eq2Left, eq2Right, eq3Left, eq3Right} /. sol, 50]

(* {1.0000000030827008893968408115631043435979982475800,
-2.7413765436049050266371405996912199223608668941690*10^-13,
1.0075000015529105718368688595084319080537124433646,
1.0075005247423274415176177009375394837002435886736,
1.0000000015413504435105398109264832834280024251757,
1.0000000772101989945419900197176287358518885269469} *)


It looks like eq1Right doesn't fit. Either something is wrong with that equation or the solution is not adequate.

If eq1Right` is correct, then setting values of one of the parameters and then solving two of the equations for the other two parameters as a search procedure might be considered.