Why is FindRoot throwing a “singular Jacobian” error when function has small values?

Running the code:

FindRoot[PolyLog[2, -E^-(10 + Sqrt[1 + x^2])] == 10^-30, {x, 10}]


Gives me the error:

FindRoot::jsing: Encountered a singular Jacobian at the point {x} = {27.0291}.
Try perturbing the initial point(s).


As far as I can tell, the Jacobian (derivative?) is not singular; perhaps something is going wrong because some very small values are involved?

Mathematica also returns the solution x → 27.0291, which is just wrong. What's the issue here and how can I avoid it?

Update

The above equation is unsolvable as is, because it should have a negative sign in front. With the negative sign, the exact same issue occurs:

FindRoot[-PolyLog[2, -E^-(10 + Sqrt[1 + x^2])] == 10^-30, {x, 10}]

>> FindRoot::jsing: Encountered a singular Jacobian at the point {x} = {27.0291}.
Try perturbing the initial point(s).

{x -> 27.0291}

• I agree it is confusing. help says This message is generated when the derivatives of the first argument in FindRoot are singular at the indicated point. At first I thought the indicated point means the starting root finding point, which is 10 in your case. But it is be the ending point (i.e. root found) (as message says). Try this root = x /. FindRoot[PolyLog[2, -E^-(10 + Sqrt[1 + x^2])] == 10^-30, {x, 10}]; der = D[-E^-(10 + Sqrt[1 + x^2]), x]; (der /. x -> root) // Chop which gives zero. So derivative of the function is zero at the root found (not at the starting point !) – Nasser Apr 17 '20 at 22:19
• Mathematica also returns the solution x → 27.0291, which is just wrong can you please explain why is this wrong? A quick plot indicates it seems to be correct. – Nasser Apr 17 '20 at 22:22
• When I plot your PolyLog function, it is negative. It can't be equal to 10^-30. You can increase the WorkingPrecision, but you still need a solvable equation. – Bill Watts Apr 17 '20 at 22:41
• 1. Expressing the Fermi-Dirac integral in terms of the polylogarithm is usually not numerically advisable, since evaluations can get unstable. 2. Are you sure your polylogarithm can become positive? Plotting PolyLog[2, -E^-(10 + Sqrt[1 + x^2])] indicates that this function is always negative for all real x. – J. M.'s ennui Apr 18 '20 at 4:19
• @Nasser -PolyLog[2, -E^-(10 + Sqrt[1 + 27.0291^2])] returns 8.13645*10^-17 – WillG Apr 18 '20 at 15:21

If we examine the Jacobian, we can see numerical problem:

D[-PolyLog[2, -E^-(10 + Sqrt[1 + x^2])], x]
(*  -(x Log[1 + E^(-10 - Sqrt[1 + x^2])])/Sqrt[1 + x^2]  *)


The argument 1 + E^(-10 - Sqrt[1 + x^2]) of Log will be equal to 1 and the Jacobian equal to zero, when 0 < E^(-10 - Sqrt[1 + x^2]) < $MachineEpsilon. This happens for all x > 26.0244: FindRoot[E^(-10 - Sqrt[1 + x^2]) ==$MachineEpsilon, {x, 1}]
(*  {x -> 26.0244}  *)


Therefore the root, which is about 59, cannot be found.

A workaround is to use InternalLog1p[x], which is equal to Log[1 + x] but avoids computing 1+x:

Block[{Log = InternalLog1p[-1 + #] &},
FindRoot[-PolyLog[2, -E^-(10 + Sqrt[1 + x^2])] == 10^-30, {x, 10}]
]
(*  {x -> 59.0691}  *)


Another, simpler way is to use arbitrary precision:

FindRoot[-PolyLog[2, -E^-(10 + Sqrt[1 + x^2])] == 10^-30, {x, 10},
WorkingPrecision -> 16]
(*  {x -> 59.06908873204437}  *)


Note we don't have to use a high precision because the precision of 1 + 1.16*^-30 is 46 by the rules of arbitrary precision. Therefore the Jacobian will be calculated accurately.

• Thanks. Both answers were good, but I like that this one pointed out the simple fix of increasing WorkingPrecision. – WillG Apr 21 '20 at 16:54

As I noted in a comment, this is one of those situations where a closed form is not helpful. In this case, forcing the Fermi-Dirac integral to be expressed as the polylogarithm of an exponential results in numerical grief.

Thus, by instead applying the formula in item 4 of this Wiki article (which can be derived using the Abel-Plana formula), we have the following computation:

prec = 50;
tmp[x_?NumberQ] := NIntegrate[Csch[π u] Sin[2 ArcTan[u] + u (10 + Sqrt[1 + x^2])]/(1 + u^2),
{u, 0, ∞}, Method -> "DoubleExponential",
WorkingPrecision -> (prec - 5)]

With[{ε = 1*^-30},
FindRoot[tmp[x] == ε Exp[10 + Sqrt[1 + x^2]] - 1/2,
{x, 50, 60}, WorkingPrecision -> prec]]
{x -> 59.069088732044366955570311882437432292701716044411}

-PolyLog[2, -Exp[-(10 + Sqrt[1 + x^2])]] /. First[%]
1.000000000000000000000000000000000000003352846600*^-30

% - 1*^-30
3.352846600*^-69


Note the use of NumericQ[] in the definition for tmp[], as well as the use of a higher WorkingPrecision setting in FindRoot[] than in NIntegrate[]`.