3
$\begingroup$

I am solving the nonlinear questions in equation 4.1a-4.1c given in the paper. (https://arxiv.org/pdf/cond-mat/0008249.pdf) but no success.

In[26]:= d = 1;
t = 10^-4 d;
ϵf0 = - 0.102 d;
p = 0.875;
q0 = 0.5;
v = 0.5 d;
ds[ϵ_] := 3/(4 Sqrt[2] d) Sqrt[(ϵ + d)/d];

(*we calculate renormalized position of the f-level, chemical \
potential and slave boson amplitude by solving three nonlinear \
equations*)

ζ[ϵ_, μ_] := ϵ - μ;  \
(*ϵ=k^2/(2 m)-d*)
rk[ϵ_, μ_, ϵf_, a_] := 
 Sqrt[(ϵf - ζ[ϵ, μ])^2 + (2 v a)^2];
enkp[ϵ_, μ_, ϵf_, a_] := 
  0.5 (ϵf + ζ[ϵ, μ] + 
     rk[ϵ, μ, ϵf, a]);
enkm[ϵ_, μ_, ϵf_, a_] := 
  0.5 (ϵf + ζ[ϵ, μ] - 
     rk[ϵ, μ, ϵf, a]);
fwp[ϵ_, μ_, ϵf_, a_] := 
  0.5 (1 - Tanh[enkp[ϵ, μ, ϵf, a]/(2 t)]);
fwm[ϵ_, μ_, ϵf_, a_] := 
  0.5 (1 - Tanh[enkm[ϵ, μ, ϵf, a]/(2 t)]);

In[38]:= diff[ϵ_, μ_, ϵf_, a_] := 
  fwm[ϵ, μ, ϵf, a] - 
   fwp[ϵ, μ, ϵf, a];
sum[ϵ_, μ_, ϵf_, a_] := 
  fwm[ϵ, μ, ϵf, a] + 
   fwp[ϵ, μ, ϵf, a];
eq1int[ϵ_, μ_, ϵf_, a_] := 
  diff[ϵ, μ, ϵf, a]/
  rk[ϵ, μ, ϵf, a];
eq2int[ϵ_, μ_, ϵf_, 
   a_] := (ζ[ϵ, μ] - ϵf) \
eq1int[ϵ, μ, ϵf, a];

eq3int[ϵ_, μ_, ϵf_, a_] := 
  sum[ϵ, μ, ϵf, a];

fint1[ϵ_, μ_, ϵf_, a_] := 
  ds[ϵ] eq1int[ϵ, μ, ϵf, a];
fint2[ϵ_, μ_, ϵf_, a_] := 
  ds[ϵ] eq2int[ϵ, μ, ϵf, a];
fint3[ϵ_, μ_, ϵf_, a_] := 
  ds[ϵ] eq3int[ϵ, μ, ϵf, a];

In[46]:= eqn1[μ_?NumericQ, ϵf_?NumericQ, 
  a_?NumericQ] := ϵf - ϵf0 - 
  v^2 NIntegrate[
    fint1[ϵ, μ, ϵf, a], {ϵ, -d, d}, 
    MaxRecursion -> 20, 
    Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]
eqn2[μ_?NumericQ, ϵf_?NumericQ, a_?NumericQ] := 
 2 (q0 - a^2 ) - p - 
  NIntegrate[
   fint2[ϵ, μ, ϵf, a], {ϵ, -d, d}, 
   MaxRecursion -> 20, 
   Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]
eqn3[μ_?NumericQ, ϵf_?NumericQ, a_?NumericQ] := 
 p - NIntegrate[
   fint3[ϵ, μ, ϵf, a], {ϵ, -d, d}, 
   MaxRecursion -> 20, 
   Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]

In[49]:= mysol = 
 FindRoot[{eqn1[μ, ϵf, a], eqn2[μ, ϵf, a], 
   eqn3[μ, ϵf, a]}, {{μ, -0.42628}, {ϵf, 
    0.035940}, {a, -0.37342}}]


During evaluation of In[49]:= FindRoot::jsing: Encountered a singular Jacobian at the point {μ,ϵf,a} = {0.845501,0.196729,-0.730513}. Try perturbing the initial point(s).

Out[49]= {μ -> 0.845501, ϵf -> 0.196729, a -> -0.730513}

In[50]:= {eqn1[μ, ϵf, a], eqn2[μ, ϵf, a], 
  eqn3[μ, ϵf, a]} /. mysol

Out[50]= {0.0631157, -0.293102, -0.125}

Also, I have looked similar problem but couldn't help me. FindRoot::jsing: Encountered a singular Jacobian at the point when solving NONLINEAR EQUATIONS

Could someone help me out?

$\endgroup$
1
  • 1
    $\begingroup$ Try to use the affine covariant Newton method $\endgroup$
    – user21
    Dec 10, 2022 at 21:22

1 Answer 1

3
$\begingroup$

Several settings do find a solution:

sol[opts___] := FindRoot[{eqn1[μ,ϵf,a],eqn2[μ,ϵf,a],eqn3[μ,ϵf,a]},
                   {{μ,-0.42628},{ϵf,0.035940},{a,-0.37342}},opts];

sol[DampingFactor -> 0.3]
(* {μ->0.446947,ϵf->0.190651,a->-0.540233} *) 

sol[Method -> {"Newton", "StepControl" -> "TrustRegion"}]
(* {μ->0.446947,ϵf->0.190651,a->-0.540233} *)

sol[Method -> "Secant"]
(* {μ->0.446947,ϵf->0.190651,a->-0.540233} and warnings *)

This point does satisfy the equations.

For a good discussion of options, such as damping, see this Mathematica SE thread.

$\endgroup$

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.