Difficulty when trying to find roots to high accuracy

Question

I'm working with the logistic map $f(x,\lambda)=4\lambda x(1-x)$, and iterations of the logistic map $f^{(2^n)}(x,\lambda)=f^{(2^{n-1})}(f^{(2^{n-1})}(x,\lambda),\lambda)$. There are some special values $\lambda_n$ which have a $2^n$ cycle, have $f^{(2^n)}(1/2,\lambda_n)=1/2$, and which also have $\frac{d}{dx} f^{(2^n)}(x,\lambda_n)|_{x=1/2}=0$ by symmetry. For example, here are three of the functions $f^{(2)}(x,\lambda_1)$, $f^{2^2)}(x,\lambda_2)$, and $f^{(2^3)}(x,\lambda_3)$,

I've found $\lambda_0$ through $\lambda_{10}$, and I've found $\lambda_{11}$ but it's not very accurate. I'd like to push things a bit further and get $n=11,12,13,\ldots$ to a higher degree of accuracy but whenever I try to add accuracy and plug things into FindRoot I get Overflow[] errors! I've also tried using FindMinimum instead of FindRoot to no avail. This is confusing since f is a function from the unit interval to the unit interval, so it's impossible for the iterated $f$ to diverge no matter what value of $\lambda\in[0,1]$ is plugged in. I imagine that the gradients of the functions can get large, but the gradients near $x=1/2$ should be on the order of $\alpha^{11}\approx 24000$ where $\alpha\approx 2.5$ is Feigenbaum's second constant. So it feels like this problem is solvable even with machine precision, and it should be easily solved by an arbitrary precision arithmetic library.

I've also restricted FindRoot and FindMinimum to the domain $[0.89,0.9]$, which is where all the rest of the $\lambda_n$ should be.

Is there any way to fix this code? How is this code producing an overflow even though the function is restricted to the interval [0,1]?

Note that $\lambda_{n+1}-\lambda_n \approx (\lambda_{n}-\lambda_{n-1})/\delta$, where $\delta\approx 4.6692016$ is Feigenbaum's first constant, so the lambda values do get very close to each other very quickly.

I'm working with Mathematica 11.3

(* Define the logistic map and iterated logistic map *)
f[x_?NumericQ,lambda_?NumericQ]:=4 lambda x(1-x);
f[n_,x_?NumericQ,lambda_?NumericQ]:=Nest[f[#,lambda]&,x,n];
(* starting estimates for the roots whose precision I'd like to improve. l[8] through l[10] are accurate to about 16 digits, and l[11] is less accurate but very close. *)
{l[8],l[9],l[10],l[11]}={0.8924846935583266`60,0.8924860486520165`60,0.8924863388716187`60,0.8924864027916384`60};
(* Try to find better approximations to the root using FindRoot *)
Table[
FindRoot[f[2^k,1/2,lambda]-1/2,{lambda,l[k],0.89,0.9},WorkingPrecision->80,PrecisionGoal->60,AccuracyGoal->60],
{k,8,11}]
(* Trying the same with FindMinimum *)
Table[
FindMinimum[(f[2^k,1/2,lambda]-1/2)^2,{lambda,l[k],0.89,0.9},WorkingPrecision->80,PrecisionGoal->60,AccuracyGoal->60],
{k,8,11}]

I should also note that without the precision arguments, things converge just fine (although lambda[11] is still very inaccurate):

{l[8], l[9], l[10], l[11]} = {0.8924846935583266, 0.8924860486520165, 0.8924863388716187, 0.8924864027916384};
Table[FindRoot[f[2^k, 1/2, lambda] - 1/2, {lambda, l[k], 0.89, 0.9}], {k, 8, 11}]

Answer 1

Increase the precision:

Table[
  With[{k = k},
   FindRoot[f[2^k, 1/2, lambda] - 1/2,
    {lambda,
     SetPrecision[l[k], Infinity],  (* changed *)
     89/100, 9/10},                 (* not necessary, just style *)
    WorkingPrecision -> 5000,       (* changed *)
    PrecisionGoal -> 60, AccuracyGoal -> 60]
   ],
  {k, 8, 11}] // SetPrecision[#, 80] &

(*
{{lambda -> 
   0.89248469355832637194836151689074673159028753656083109926743491269675432307574010},
 {lambda -> 
   0.89248604865201623331090833596804333411752419940218509930716741401006794017472219},
 {lambda -> 
   0.89248633887161714522320059557704106220409625694389351599123512202281630080794540},
 {lambda -> 
   0.89248640102776960953352931534540018803062500172009213180415107613512736663560703}}
*)

---

Response to comment: Another approach.

Fixing the precision is faster and prevents the loss of precision that leads to overflow. The algorithm in FindRoot corrects in the next step any loss of precision in a previous step.

Block[{$MinPrecision = 80, $MaxPrecision = 80},
 Table[
  With[{k = k},
   FindRoot[f[2^k, 1/2, lambda] - 1/2,
    {lambda,
     l[k],
     89/100, 9/10},
    WorkingPrecision -> 80, PrecisionGoal -> 60, AccuracyGoal -> 60]
   ],
  {k, 8, 11}]
 ]

The output is the same as above, except for the last digit in the solution for k == 8.

Update/addendum: I should add that finding the root would be futile if f[] is not calculated accurately with 80-digit-precision bignums. The reason I know it worked is by comparing with the 5000-digit calculation. What is happening is that the uncertainty bound is being magnified in Nest[f[#,lambda]&,x,n] and is much bigger than the actual error. The Accuracy loses a little over a half digit at each iteration until the Accuracy becomes negative; at that point it soon grows exponentially, doubling at each step until Overflow[] is reached. [More precisely: The Accuracy of a arbitrary-precision "bignum" is equal to -Log10[uncertaintybound]. When one starts a calculation, the uncertainty bound is prescribed by the precision of the input; see this guide, which points out several related tutorials and commands. This bound is calculated and propagate through a computation. In this case, Log10[uncertaintybound] increases by 0.55 at each iteration of f[] until uncertaintybound > 1 (the Accuracy becomes negative); at that point Log10[uncertaintybound] soon doubles at each step, which means uncertaintybound grows proportional to Exp[2^n] until Overflow[] is reached.]