Asymptotic behavior of iterative sequences

Although I don't have a "final" answer, this suggestion may help.

For polynomials or for analytic functions (having a power series representation with nonzero radius of convergence) I'd employ the concept of Carleman-matrices.
Assume a vectorfunction $V(x) = [1,x,x^2,x^3,...]$ as rowvector and $F$ as carleman-matrix (transposed) for your function $f(x)$ and I for the identity-matrix then we could in principle write

$\qquad V(a_1) \cdot I = V(a_1)$
$\qquad V(a_1) \cdot F = V(f(a_1)) $
but in
$\qquad V(a_1) \cdot (I+F) = V(a)+ V(f(a_1)) \ne V(a_2)$
the sum of two V()-vectors is not a V()-vector.

Instead we define first the Carleman-matrix $G$ for the function $g(x)=x+f(x)$

Then we can iterate:

$\qquad V(a_0) = V(a_0) \cdot I\\ \qquad V(a_1) = V(a_0) \cdot G \\ \qquad V(a_2) = V(a_1) \cdot G = V(a_0) \cdot G^2 \\ \qquad \cdots \\ \qquad V(a_k) = V(a_0) \cdot G^k \\ $
as long as taking the k'th power $G^k $ makes sense (requires only convergent (or as generalization for certain divergent cases for instance Euler-summable) series).

If $G$ is triangular, the formal power series for your iterated expression $a_k$ can exactly be given to any power (even for fractional powers!) and with your initial value $a_o$ it might be convergent up to some highest power of G . For instance, the transposed Carleman-matrix for the function $ f(x) = \ln(1+x)$ is triangular, and also that for $g(x) = x + \ln(1+x) $, however the range of convergence of the formal power series decreases with the iterations...
If $G$ is triangular and $g(x)' = 1$ we can express the power series for the k'th iteration with coefficients, which are polynomials in k and are thus especially easy to compute.
If $G$ is triangular and $g(x)' \notin (0,1) $ we can apply matrix-diagonalization, which gives again exact coefficents for the power series, but are likely more complicated.
(See two examples below)

For non-triangular Carleman-matrices (where $g(0) \ne 0$) this is even more delicate and exceeds the focus of this answer...

---

Two examples. For $f(x) = \ln(1+x)$ the power-series for the *h*-fold iterate ($f°^h(x)$) begins with

$ \small \begin{eqnarray} f°^h(x) =& x \cdot &(1) \\ &+x^2 \cdot &-( 1/2 h) \\ &+x^3 \cdot & (1/12 h&+1/4 h^2 ) \\ &+x^4 \cdot & -( 1/48 h &+ 5/48 h^2& +1/8 h^3) \\ &+x^5 \cdot & ( 1/180 h& + 1/24 h^2& + 13/144 h^3&+1/16 h^4 ) \\ &+ O(x^6) \end{eqnarray} $

and for $g(x)= x + \ln(1+x) $ the powerseries for the h-fold iterate begins with

$ \small \begin{eqnarray} g°^h(x) =& x \cdot &(1 u) \\ &+x^2 \cdot &(1/4 u-1/4 u^2) \\ & +x^3 \cdot &(1/36 u-1/8 u^2+7/72 u^3) \\ & +x^4 \cdot &(1/672 u-17/576 u^2+7/96 u^3-181/4032 u^4) \\ & +x^5 \cdot &(11/75600 u-17/4032 u^2+91/3456 u^3-181/4032 u^4+13687/604800 u^5)\\ &+O(x^6) \\ \end{eqnarray} $
where we have to replace $u$ by $2^h$ and h is the iteration-height-parameter.