0
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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}

$\endgroup$
  • $\begingroup$ This gives me the same (presumably wrong) values as the OP for b, x, and y... $\endgroup$ – MelaGo Sep 4 at 21:06
  • $\begingroup$ 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}. $\endgroup$ – user1 Sep 5 at 12:10
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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