Encountered a singular Jacobian while solving nonlinear equations

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?

• Try to use the affine covariant Newton method Commented Dec 10, 2022 at 21:22

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.