# Errors using FindRoot on slow numerical function

I would like to find the solution of an equation containing a slow function

FindRoot[func[x, onsideS, hoppingS, 0.01] - (1/E) == 0, {x, 1,10}]


for different values of onsideS and hoppingS. func is given by

func[tau_?NumericQ, onsideS_?NumericQ, hoppingS_?NumericQ, eta_?NumericQ] :=
func[tau, onsideS, hoppingS, eta] =
Min[Exp[eta*tau]*Abs[NIntegrate[Exp[-I*x*tau]*(1/(2 Pi))*((1/(x - onsideS + I*eta))
+ ((1/(x - onsideS + I*eta))^2)*(hoppingS^2)*localGF[x,eta]*
(1/(1 - (hoppingS^2) * (1/(x - onsideS + I*eta))*localGF[x, eta]))), {x, -Infinity, Infinity},
Method -> {"LevinRule", "Kernel" -> Exp[-I*x*tau]}], 1]


with

localGF[energy_?NumericQ, eta_?NumericQ] := localGF[energy, eta] =
NIntegrate[ldos[x]/(energy - x + I*eta), {x, -hbw, hbw},
Method -> {"GaussKronrodRule", "Points" -> 25,"SymbolicProcessing" -> 0}]
ldos[energy_?NumericQ] := (1/(Pi))*Sqrt[1 - (energy/2)^2]
hbw = 2;


However Mathematica always returns a bunch of error messages containing of

NIntegrate: Integrand is not a Levin function

FindRoot: Failed to converge to the requested accuracy or precision within 100 iterations

and FindRoot never gives a valid result, even though calling func manually returns something useful. The same happens when trying to Plot this function. I tried using NSolve aswell

NSolve[func[x, onsideS, hoppingS, 0.01] - 1/E == 0, x, Method -> {Automatic, "SymbolicProcessing" -> 0}]


and it works atleast for certain values (i.e. onsideS = 0 and hoppingS = 1) quite fast, but evaluation time increases very quickly when going to other values (i.e. onsideS = 0 and hoppingS = 0.05).

Why does this happen and how could I fix this?

• what value have hbw? – Mariusz Iwaniuk Jul 7 '18 at 17:52
• @MariuszIwaniuk Sorry, I forgot adding it into the question. hbw = 2; – StirriX Jul 7 '18 at 18:17
• The error message you get suggests that localGF has not been passed into the FindRoot properly, I can't see why immediately though. – KraZug Jul 8 '18 at 6:03

This is not a complete answer, but some hints.

First, do as much integrations, you can, analytically.

hbw = 2;

ldos[energy_] = (1/(Pi))*Sqrt[1 - (energy/2)^2] // Simplify;

localGF[energy_, eta_] =
Integrate[ldos[x]/(energy - x + I*eta), {x, -hbw, hbw},
Assumptions -> eta > 0 && energy \[Element] Reals]

(*   -(1/2) (-1 + Sqrt[1 - 4/(energy + I eta)^2]) (energy + I eta)   *)


Then regard the integrand for NIntegrate at different parameters.

integrand1[tau_, onsideS_, hoppingS_, eta_, x_] =
Exp[-I*x*tau]*(1/(2 Pi))*((1/(x - onsideS +
I*eta)) + ((1/(x - onsideS + I*eta))^2)*(hoppingS^2)*
localGF[x,
eta]*(1/(1 - (hoppingS^2)*(1/(x - onsideS + I*eta))*
localGF[x, eta]))) // Simplify;

Manipulate[
int = integrand1[tau, onsideS, hoppingS, 1/100, x]; {Plot[
Re@int, {x, -10, 10}], Plot[Im@int, {x, -10, 10}]}, {tau, 1,
10}, {onsideS, 0, 4}, {hoppingS, 0, 2}]


For some parameter combinations, NIntegrate works quite good

NIntegrate[integrand1[1, 0, 1, 1/100, x], {x, -Infinity, Infinity}]

(*   2.77556*10^-17 - 0.57145 I   *)

NIntegrate[integrand1[2, 0, 1, 1/100, x], {x, -Infinity, Infinity}]

(*   0. + 0.0321347 I   *)


For others, like shown in picture, there is a non integrable singularity at zero or other x-values, you have to use Method -> "PrincipalValue" an explicitly give the x-value of it in integration range.

NIntegrate[integrand1[1, 0, 0, 1/100, x], {x, -Infinity, 0, Infinity},
Method -> "PrincipalValue"]

(*   0. - 0.99005 I   *)


For other parameter values, integral failes to converge.

NIntegrate[
integrand1[1, 0, 2/10, 1/100, x], {x, -Infinity, 0, Infinity},
Method -> "PrincipalValue"]


It's a tricky job to further examine integrand behavior and positions of singularities.