# Problem with FindRoot function

I have set up a code to solve simultaneous integral equations. These result essentially from taking integrals of the Fermi-Dirac distribution. Although there are no obvious problems with the code it does not give the correct answers for the roots. I have checked that in the limiting case when I could evaluate the integrals analytically and find the roots exactly. Here is the code :

dv1 = 1
dv2 = 1
u1 = 5
u2 = 8
η1 = 1
η2 = 0.5
N1 = 2 * η1 * dv1
N2 = 2 * η2 * dv2
K = η1 * (dv1 ^ 3 / 12 + u1 ^ 2 * dv1) + η2 * (dv2 ^ 3 / 12 + u2 ^ 2 * dv2)
f[b_, x_, y_] := 1/Sqrt[2*η1]*NIntegrate[b^(-1/2)*z^(-1/2)*η1*Exp[-z + b*x*η1]/(1 + Exp[-z + b*x*η1] + Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}]
g[b_, x_, y_] := 1/Sqrt[2*η2]*NIntegrate[b^(-1/2)*z^(-1/2)*η2*Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] + Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]
h[b_, x_, y_] :=
2^(-1/2)*(η1)^(-3/2)*NIntegrate[b^(-3/2)*z^(1/2)*η1*Exp[-z + b*x*\[Eta]1]/(1 + Exp[-z + b*x*η1] + Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}] +
2^(-1/2)*(η2)^(-3/2)*NIntegrate[b^(-3/2)*z^(1/2)*η2*Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] + Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]
sol = FindRoot[{N1 - f[b, x, y] == 0, N2 - g[b, x, y] == 0, K - h[b, x, y] == 0}, {{b, 0.0349599}, {x, -54.4214}, {y, -128.67}}, WorkingPrecision -> 40] /. Rule -> List
b = Re[sol[[1]][[2]]]
x = Re[sol[[2]][[2]]]
y = Re[sol[[3]][[2]]]


I am trying to solve the equations for the variables b,x,y. The problem seems to be with the FindRoot part of the code. Even if I specify as the initial values for these variables the values that are the true solutions I do not get them in return. As I am only a beginner in Mathematica I don't have an idea how to proceed now. There might be more than one root to the system of equations but I don't know for sure.

Need to change N1 to n1 as N is a built-in function name. Also need to change function arguments to x_?NumericQ as I did below. Then change 0.5 to 1/2 to set up arbitrary-precision calculations. Then this seems to work:

dv1 = 1
dv2 = 1
u1 = 5
u2 = 8
η1 = 1
η2 = 1/2
n1 = 2*η1*dv1
n2 = 2*η2*dv2
K = η1*(dv1^3/12 + u1^2*dv1) + η2*(dv2^3/12 + u2^2*dv2)
f[b_?NumericQ, x_?NumericQ, y_?NumericQ] :=
1/Sqrt[2*η1]*
NIntegrate[
b^(-1/2)*z^(-1/2)*η1*
Exp[-z + b*x*η1]/(1 + Exp[-z + b*x*η1] +
Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}]

g[b_?NumericQ, x_?NumericQ, y_?NumericQ] :=
1/Sqrt[2*η2]*
NIntegrate[
b^(-1/2)*z^(-1/2)*η2*
Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] +
Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]

h[b_?NumericQ, x_?NumericQ, y_?NumericQ] :=
2^(-1/2)*(η1)^(-3/2)*
NIntegrate[
b^(-3/2)*z^(1/2)*η1*
Exp[-z + b*x*η1]/(1 + Exp[-z + b*x*η1] +
Exp[-(η2/η1)*z + b*y*η2]), {z, 0, Infinity}] +
2^(-1/2)*(η2)^(-3/2)*
NIntegrate[
b^(-3/2)*z^(1/2)*η2*
Exp[-z + b*y*η2]/(1 + Exp[-z + b*y*η2] +
Exp[-(η1/η2)*z + b*x*η1]), {z, 0, Infinity}]

sol = FindRoot[{n1 - f[b, x, y] == 0, n2 - g[b, x, y] == 0,
K - h[b, x, y] == 0}, {{b,
0.0349599}, {x, -54.4214}, {y, -128.67}},
WorkingPrecision -> 40] /. Rule -> List
b = Re[sol[[1]][[2]]]
x = Re[sol[[2]][[2]]]
y = Re[sol[[3]][[2]]]
theSol = #[[2]] & /@ sol
Apply[{(n1 - f[#1, #2, #3]), n2 - g[#1, #2, #3],
K - h[#1, #2, #3]} &, theSol]*emphasized text*


{0.,3.33067*10^-16,1.42109*10^-14}

• This gives me the same (presumably wrong) values as the OP for b, x, and y... Sep 4, 2019 at 21:06
• To clarify, the correct values for the variables are : {b,x,y}={0.0349599,-54.4214,-128.67}. Those were computed analytically and are the correct values for the specified values of the input parameters (dv1, dv2, u1, u2,). I have just checked the above code but it also gives those incorrect values. Nevermind what values you start with as the input values for the {b,x,y} (even if you use the correct values) you still always get the same incorrect values {0.0433179826,-14.02830453,-47.91130641}. Sep 5, 2019 at 12:10

We have the three simultaneous equations:

\begin{align} n1-f(b,x,y)&=0\\ n2-g(b,x,y)&=0\\ K-h(b,x,y)&=0 \end{align}

Now when I run the code above, I get:

{0.04331798263222832287097452988786750625601,
-14.02830453956710576253245993437615631186,
-47.91130641078542149075012315062585520063}


and back-substituting those values into the three equations returns:

{0., 3.33067*10^-16, 1.42109*10^-14}


However, when I back-substitute your solutions:

{0.0349599,-54.4214,-128.67}


I do not get answers close to zero:

{0.841642, 0.432857, 26.3829}


so must conclude the problem is not with FindRoot but rather before it, that is, the integrations perhaps but do not receive any type of error message with the integrations.