How to get $f(x)$, if we know that $f(f(x))=x^2+x$?
Is there an elementary function $f(x)$ that satisfies the equation?
2026-03-28 20:08:14.1774728494
On
On
On
On
To have a function with another fractional iteration height I use like in the logic of fractional powers of scalars:
We can get this even symbolically, where $h$ is left indeterminate, all in he same small framework:
It might be interesting, that we can partially reproduce Will Jagy's computation with the resources which we have computed so far.
The $\lambda()$ - function can be taken from the coefficients of the second column of the $\log$-matrix $L$, and then it needs two more steps to arrive at the Abel-function $\alpha()$:
On
On
How to obtain $f(x)$, if it is known that $f(f(x))=x^2+x$?
5.2k Views Asked by Bumbble Comm https://math.techqa.club/user/bumbble-comm/detail At
7
There are 7 best solutions below
8
On
There is no such $f:\mathbb C\to\mathbb C$. See this paper:
When is $f(f(z)) = az^2 + bz + c$ ?
by R. E. Rice, B. Schweizer and A. Sklar
The American Mathematical Monthly, vol. 87, no. 4 (Apr., 1980), pp. 252–263
More generally, they prove that a quadratic polynomial has no iterative roots of any order.
0
On
Firstly, let $g(x)$ be equal to $x^2+x$. Now we can say that $g(x)=g(y)$ is equivalent to $x=y$ or $x+y=-1$. So $g(g(x))=g(g(y))$ means that $g(x)=g(y)$ or $g(x)+g(y)=-1$. But $g(x)=\left(x+\frac{1}{2}\right)^2-\frac{1}{4}\ge-\frac{1}{4}$, so the second case can't take place. Thus, $g^n(x)=g^n(y)$ iff $g(x)=g(y)$ for every positive integer $n$.
Also $f$ returns each number not less than $-\frac{1}{4}$.
Well, I don't have a full solution and I print it on my little mobile phone, so my idea is the following: all real numbers can be divided onto many infinitive sequences, some of them are also infinitive to the left, every number is a value of $g$ at the previous number in the same sequence (if it exists), and some numbers don't occur in any sequence, but their values of $g$ are in one of our sequences. In other words, we draw an arrow from $x$ to $g(x)$ for every $x$. After all I've written and your imagination we can understand our situation. Now $f$ can just divide there sequences by pairs and map them to each other at every pair. Sorry for not too clear explanation
8
On
Here’s a technique for finding the first few terms of a formal power series representing the fractional iterate of a given function like $f(x)=x+x^2$. I repeat that this is a formal solution to the problem, and leaves unaddressed all considerations of convergence of the series answer.
I’m going to find the first six terms of $f^{\circ1/2}(x)$, the “half-th” iterate of $f$, out to the $x^5$-term. Let’s write down the iterates of $f$, starting with the zero-th. \begin{align} f^{\circ0}(x)&=x\\ f^{\circ1}=f&=x&+x^2\\ f^{\circ2}&=x&+2x^2&+2x^3&+x^4\\ f^{\circ3}&\equiv x&+3x^3&+6x^3& + 9x^4& + 10x^5& + 8x^6\\ f^{\circ4}&\equiv x &+ 4x^2& + 12x^3& + 30x^4& + 64x^5& + 118x^6\\ f^{\circ5}&\equiv x& + 5x^2& + 20x^3& + 70x^4& + 220x^5& + 630x^6\\ f^{\circ6}&\equiv x& + 6x^2& + 30x^3& + 135x^4& + 560x^5& + 2170x^6\\ f^{\circ7}&\equiv x& + 7x^2& + 42x^3& + 231x^4& + 1190x^5& + 5810x^6\,, \end{align} where the congruences are modulo all terms of degree $7$ and more.
Now look at the coefficients of the $x$-term: always $1$. Of the $x^2$-term? In $f^{\circ n}$, it’s $C_2(n)=n$. The coefficient of $x^3$ in $f^{\circ n}$ is $C_3(n)=n(n-1)=n^2-n$, as one can see by inspection. Now, a moment’s thought (well, maybe several moments’) tells you that $C_j(n)$, the coefficient of $x^j$ in $f^{\circ n}$, is a polynomial in $n$ of degree $j-1$. And a familiar technique of finite differences shows you that \begin{align} C_4(n)&=\frac{2n^3-5n^2+3n}2\\ C_5(n)&=\frac{3n^4-13n^3+18n^2-8n}3\,, \end{align} I won’t go into the details of that technique. The upshot is that, modulo terms of degree $6$ and higher, you have $f^{\circ n}(x)\equiv x+nx^2+(n^2-n)x^3+\frac12(2n^3-5n^2+3n)x^4+\frac13(3n^4-13n^3+18n^2-8n)x^5$.
Now, you just plug in $n=\frac12$ in this formula to get your desired series. And I’ll leave it to you to go one degree higher, using the iterates I’ve given you.
1
On
Although the question has already been satisfyingly answered I'd like to add some more general framework for questions like this. For polynomials and power series there is the concept of Bell- and Carleman-matrices.
The basic idea is that if you have a Carleman-matrix $F$ associated to a function $f(x)$ then the function $g(x)=f°^{1/2}(x)$ is associated to (matrix-) squareroot $G=F^{1/2}$ (which is again a Carlemanmatrix).
This plays all around only on the coefficients of the formal power series ("FPS") and the typical cases are that of polynomials and/or power series without constant term: $$\begin{eqnarray} f(x)&=&\sum_{k=1}^N a_k \cdot x^k &\qquad \qquad \text{ or }\\ f(x)&=&\sum_{k=1}^\infty a_k \cdot x^k \end{eqnarray}$$ (but can be generalized)
With some relatively simple very general and standard user-procedures using the software Pari/GP I express your problem simply by
n = 16 \\ setting dimension for matrices and vectors
F = carleman(x+x^2,n) \\ making F a carlemanmatrix for f(x)
G = SQRT(F,1,1) \\ matrix-squareroot, flags 1,1 indicate specialized
\\ routine for lower triangular matrices
\\ with units in the diagonal
g = Ser(G[,2]) + O(x^12) \\ extract from G's second column its n coefficients
\\ and write them as FPS with the first 12 terms
getting the leading terms of a formal power series
g(x) = x + 1/2*x^2 - 1/4*x^3 + 1/4*x^4 - 5/16*x^5 + 27/64*x^6 - 9/16*x^7
+ 171/256*x^8 - 69/128*x^9 - 579/2048*x^10 + 10689/4096*x^11 + O(x^12)
in a blink...
To have a function with another fractional iteration height I use like in the logic of fractional powers of scalars:
n = 16 \\ setting dimension for matrices and vectors
h = 3/4 \\ setting some example fractional iteration-height
F = carleman(x+x^2,n) \\ making F a carlemanmatrix for f(x)
L = LOG(F,1,1) \\ matrix-logarithm, flags 1,1 indicate specialized
\\ routine for lower triangular matrices
\\ with units in the diagonal
G = EXP( h * L,1,1) \\ by the Exponential this gives the fractional
\\ h'th power of a triangular matrix,
\\ flags with the equivalent meaning as in LOG
g = Ser(G[,2]) + O(x^12) \\ extract from G's second column its n coefficients
\\ and write them as FPS with the first 12 terms
giving, for this instance with $h=3/4$
g(x) = x + 3/4*x^2 - 3/16*x^3 + 9/64*x^4 - 35/256*x^5 + 35/256*x^6
- 449/4096*x^7 - 1/2048*x^8 + 19041/65536*x^9 - 461901/524288*x^10
+ 1870803/1048576*x^11 + O(x^12)
We can get this even symbolically, where $h$ is left indeterminate, all in he same small framework:
n = 16 \\ setting dimension for matrices and vectors
h = 'h \\ resetting h to its symbolical use
F = carleman(x+x^2,n) \\ making F a carlemanmatrix for f(x)
L = LOG(F,1,1) \\ matrix-logarithm, flags 1,1 indicate specialized
\\ efficient routine for lower triangular matrices
\\ with units in the diagonal (exact rational arithmetic
\\ is possible
G = EXP( h * L,1,1) \\ by the Exponential this gives the fractional h'th
\\ power of a triangular matrix. Same flags
\\ indicate possibility for efficient exact computation
g = Ser(G[,2]) + O(x^n) \\ extract from G's second column its n coefficients
\\ and write them as FPS with the first n terms
giving the two parametric function $g()$ with $x$ , and $h$ for the iteration-height:
g(x,h)= x
+ h*x^2
+ (h^2 - h)*x^3
+ (h^3 - 5/2*h^2 + 3/2*h)*x^4
+ (h^4 - 13/3*h^3 + 6*h^2 - 8/3*h)*x^5
+ O(x^6)
with the coefficients at the powers of $x$ as polynomials in $h$ as already shown in @Lubin's answer.
While the logic of the frational iteration of (-of course- certain classes of) functions is thus reduced to the simple analogue of scalar fractional powers, there are some formal requirements and conditions to make this really helpful.
For instance, while it is as easy as shown to find the leading terms of a formal power series (as far as $f(x)$ has no constant term) it does not mean, that that resulting power series has a finite convergence radius or a nonzero radius of convergence at all, so in classical terms it might be merely useless. (Sometimes one can apply Euler-, Borel- or Noerlund summation for a divergent series to get approximate results anyway, but this is not the focus here. @Will Jagy's answer above shows a related (and might be even better) way to arrive at approximations)
If the basic function $f(x)$ has still no constant term but a coefficient $a_1 \notin \{0,1\}$ the analogue can be done by diagonalization, so for instance for something like $f(x)=1/2 x + 3/4 x^2 - x^3$ we need tell our LOG and EXP-routine, that they shall do things different instead, and shall use diagonalization and apply the iteration height on the exponents of the eigenvalues to create the fractional power of the Carlemanmatrix G for the desired iteration height. But again: while we'll get the correct leading terms of a power series, in most cases that power series for fractional iteration heights shall have radius of convergence being zero, so it cannot be evaluated for any $x$ - except we
- a) introduce methods for divergent summation or
- b) reduce ourselves to the evaluation of the asymptotic formal power series only up to some optimal truncation.
[update Jun 2015]:
Here is a sample-computation, parallel to that of @Will Jagy.
The formal power series for $g(x)$ (=$f(x)$ in the OP and $h(x)$ in Will Jagy's post) is surely only asymptotic (meaning: has convergence radius zero).
But by the functional relation we have $x_{0.5} = ((x_{-20})_{+0.5})_{+20}$ so that we replace $ x' = x_{-20} $ which is closer to the fixpoint zero. If we insert this value in the asymptotic formal power series we have after, say, 64 terms a partial sum for $g(x')$ (locally) converging to some value whose thirty or forty digits stay constant. If we truncate the series we have a good estimate (but, due to the characteristic of the series as asymptotic only) cannot much be improved (except we would replace $x' = x_{-40}$ or an even higher iterate) . Then we use the functional equation again to iterate now 20 times $g(x) = g(x')_{+20}$ . That this procedure is a meaningful approximation shows the following table, where the error is already very small (and can be made smaller as is wished).
The table is computed using the iteration of $x$ towards the fixpoint until $x' <0.01$ (call the required number of inverse iterations "height-offset") , then $64$ terms of the FPS are taken, and the result $g(x')$ was iterated away from the fixpoint zero by the same "height offset".
x g(x) g(g(x)) (x+x^2)- g(g(x))
------------------------------------------------------------------
0.100000000000 0.104772246757 0.110000000000 1.77368712134E-89
0.200000000000 0.218321237354 0.240000000000 8.39980200859E-89
0.300000000000 0.339733963915 0.390000000000 2.06063470023E-88
0.400000000000 0.468317683702 0.560000000000 3.01316590953E-88
0.500000000000 0.603524735182 0.750000000000 5.56377802786E-88
0.600000000000 0.744908688889 0.960000000000 8.30936184814E-88
0.700000000000 0.892096937726 1.19000000000 1.08442417625E-87
0.800000000000 1.04477260629 1.44000000000 1.74985488654E-87
0.900000000000 1.20266208158 1.71000000000 1.94506810567E-87
1.00000000000 1.36552610963 2.00000000000 2.78678464389E-87
1.10000000000 1.53315324992 2.31000000000 3.83614094010E-87
1.20000000000 1.70535494330 2.64000000000 3.76726767336E-87
1.30000000000 1.88196171736 2.99000000000 4.88855868818E-87
1.40000000000 2.06282021339 3.36000000000 6.20513270304E-87
1.50000000000 2.24779082048 3.75000000000 7.73029893345E-87
1.60000000000 2.43674576629 4.16000000000 6.98399441754E-87
1.70000000000 2.62956755754 4.59000000000 8.44163571068E-87
1.80000000000 2.82614769218 5.04000000000 1.00801098914E-86
1.90000000000 3.02638558534 5.51000000000 1.19080669499E-86
2.00000000000 3.23018766572 6.00000000000 1.39339232195E-86
2.10000000000 3.43746660925 6.51000000000 1.61658749748E-86
2.20000000000 3.64814068433 7.04000000000 1.86119108535E-86
2.30000000000 3.86213318872 7.59000000000 2.12798232294E-86
2.40000000000 4.07937196232 8.16000000000 2.41772186544E-86
2.50000000000 4.29978896302 8.75000000000 2.01254742998E-86
2.60000000000 4.52331989580 9.36000000000 2.26122901762E-86
2.70000000000 4.74990388660 9.99000000000 2.52837923787E-86
2.80000000000 4.97948319435 10.6400000000 2.81450750374E-86
2.90000000000 5.21200295561 11.3100000000 3.12011260578E-86
3.00000000000 5.44741095729 12.0000000000 3.44568319447E-86
3.10000000000 5.68565743345 12.7100000000 3.79169822947E-86
3.20000000000 5.92669488324 13.4400000000 4.15862739869E-86
3.30000000000 6.17047790709 14.1900000000 4.54693150987E-86
3.40000000000 6.41696305883 14.9600000000 4.95706285700E-86
3.50000000000 6.66610871201 15.7500000000 5.38946556386E-86
3.60000000000 6.91787493848 16.5600000000 5.84457590632E-86
3.70000000000 7.17222339803 17.3900000000 6.32282261544E-86
3.80000000000 7.42911723769 18.2400000000 6.82462716263E-86
3.90000000000 7.68852099970 19.1100000000 7.35040402841E-86
4.00000000000 7.95040053722 20.0000000000 7.90056095603E-86
4.10000000000 8.21472293692 20.9100000000 8.47549919096E-86
4.20000000000 8.48145644779 21.8400000000 6.68807161944E-86
4.30000000000 8.75057041558 22.7900000000 7.14874364065E-86
4.40000000000 9.02203522214 23.7600000000 7.62849696346E-86
4.50000000000 9.29582222949 24.7500000000 8.12760832582E-86
4.60000000000 9.57190372781 25.7600000000 8.64635043990E-86
4.70000000000 9.85025288737 26.7900000000 9.18499212034E-86
4.80000000000 10.1308437137 27.8400000000 9.74379840606E-86
4.90000000000 10.4136510060 28.9100000000 1.03230306763E-85
5.00000000000 10.6986503181 30.0000000000 1.09229467614E-85
5.10000000000 10.9858179226 31.1100000000 1.15438010481E-85
5.20000000000 11.2751307766 32.2400000000 1.21858445800E-85
5.30000000000 11.5665664903 33.3900000000 1.28493251542E-85
5.40000000000 11.8601032971 34.5600000000 1.35344874129E-85
5.50000000000 12.1557200259 35.7500000000 1.42415729315E-85
5.60000000000 12.4533960747 36.9600000000 1.49708203024E-85
5.70000000000 12.7531113865 38.1900000000 1.57224652159E-85
5.80000000000 13.0548464255 39.4400000000 1.64967405377E-85
5.90000000000 13.3585821561 40.7100000000 1.72938763825E-85
6.00000000000 13.6643000215 42.0000000000 1.81141001853E-85
6.10000000000 13.9719819250 43.3100000000 1.89576367700E-85
6.20000000000 14.2816102113 44.6400000000 1.98247084146E-85
6.30000000000 14.5931676494 45.9900000000 2.07155349145E-85
6.40000000000 14.9066374158 47.3600000000 2.16303336434E-85
6.50000000000 15.2220030795 48.7500000000 2.25693196111E-85
6.60000000000 15.5392485870 50.1600000000 2.35327055204E-85
6.70000000000 15.8583582485 51.5900000000 2.45207018209E-85
6.80000000000 16.1793167241 53.0400000000 2.55335167614E-85
6.90000000000 16.5021090122 54.5100000000 2.65713564400E-85
7.00000000000 16.8267204365 56.0000000000 2.76344248533E-85
7.10000000000 17.1531366352 57.5100000000 2.87229239424E-85
7.20000000000 17.4813435501 59.0400000000 2.98370536390E-85
7.30000000000 17.8113274161 60.5900000000 3.09770119088E-85
7.40000000000 18.1430747514 62.1600000000 3.21429947937E-85
7.50000000000 18.4765723482 63.7500000000 3.33351964530E-85
7.60000000000 18.8118072635 65.3600000000 3.45538092028E-85
7.70000000000 19.1487668109 66.9900000000 3.57990235540E-85
7.80000000000 19.4874385520 68.6400000000 3.70710282498E-85
7.90000000000 19.8278102890 70.3100000000 3.83700103014E-85
8.00000000000 20.1698700566 72.0000000000 3.96961550227E-85
8.10000000000 20.5136061155 73.7100000000 4.10496460640E-85
8.20000000000 20.8590069452 75.4400000000 4.24306654448E-85
8.30000000000 21.2060612373 77.1900000000 4.38393935851E-85
8.40000000000 21.5547578895 78.9600000000 3.33657414671E-85
8.50000000000 21.9050859994 80.7500000000 3.44443732837E-85
8.60000000000 22.2570348583 82.5600000000 3.55437811314E-85
8.70000000000 22.6105939466 84.3900000000 3.66640929124E-85
8.80000000000 22.9657529274 86.2400000000 3.78054354602E-85
8.90000000000 23.3225016420 88.1100000000 3.89679345586E-85
9.00000000000 23.6808301045 90.0000000000 4.01517149613E-85
9.10000000000 24.0407284977 91.9100000000 4.13569004096E-85
9.20000000000 24.4021871676 93.8400000000 4.25836136512E-85
9.30000000000 24.7651966200 95.7900000000 4.38319764567E-85
9.40000000000 25.1297475154 97.7600000000 4.51021096373E-85
9.50000000000 25.4958306654 99.7500000000 4.63941330610E-85
9.60000000000 25.8634370285 101.760000000 4.77081656687E-85
9.70000000000 26.2325577066 103.790000000 4.90443254899E-85
9.80000000000 26.6031839409 105.840000000 5.04027296580E-85
9.90000000000 26.9753071088 107.910000000 5.17834944249E-85
10.0000000000 27.3489187203 110.000000000 5.31867351757E-85
It might be interesting, that we can partially reproduce Will Jagy's computation with the resources which we have computed so far.
The $\lambda()$ - function can be taken from the coefficients of the second column of the $\log$-matrix $L$, and then it needs two more steps to arrive at the Abel-function $\alpha()$:
lam = Ser(-L[,2]) \\ the minus-sign indicates that we want the
\\ log of the inverse f(x): log of sqrt(1+x/4)-1/2
lam_rec = 1/lam \\ Pari/GP allows to compute the formal reciprocal
\\ in the following formal integral the 1/x-term
\\ must be removed as Pari/GP is unable to include
\\ a formal expression for log(x):
alpha = intformal(lam_rec + 1/'x) - logx \\ lx = log(x)
\\ "logx" means, we must further work with that term
Check it out:
lam = Ser( - L[,2])
%995 = -x^2 + x^3 - 3/2*x^4 + 8/3*x^5 - 31/6*x^6 + 157/15*x^7 - 649/30*x^8 + O(x^9)
lam_rec = 1/lam
%996 = -x^-2 - x^-1 + 1/2 - 2/3*x + 13/12*x^2 - 113/60*x^3 + 1187/360*x^4 - 1754/315*x^5 + O(x^6)
alpha = intformal(lam_rec + 1/'x) - logx \\ lx = log(x)
%998 = x^-1 - logx + 1/2*x - 1/3*x^2 + 13/36*x^3 - 113/240*x^4 + O(x^5)
0
On
Well, one more answer, the most q&d trick, giving the formal power series for the half-iterate using the Newton-squareroot-algorithm applied to formal power series. (Thus it is in principle the same logic as the Carleman-ansatz as in my earlier answer but looks stunningly simpler).
In Pari/GP one has the builtin-function "serreverse(f)" finding the reverse of a formal power series (which must not have a constant term as is in your problem).
So we do the following
Z(x) = x + x^2 \\ define the function of which we want the half-iterate
g = x + O(x^32) \\ declare g as formal power series as initial "value"
for(k = 2, 7, g = (Z(serreverse(g))+ g)/2 ) \\ just iterate several times
print(g + O(x^9)) \\ correct to the seventh term:
Result:
x + 1/2*x^2 - 1/4*x^3 + 1/4*x^4 - 5/16*x^5 + 27/64*x^6 - 9/16*x^7
+ 357/512*x^8 + O(x^9)
The eigth term were correct ( 171/256 x^8 ) if we had iterated one more time.
In the same way we would get the formal power series of the half-iterate for the notorious case $g(g(x)) = \exp(x)-1$ just by initializing $Z=exp(x)-1$ in the above code.
0
On
If we assume a function which has a power series expansion, we can just assume f to have a Maclaurin polynomial. Then the coefficients of $f(f(x))$ is a linear combination of iterated convolutions, for which the right hand side is $[0,1,1,0,0,\cdots]$. Also some regularization should likely be employed, punishing too large coefficients for large exponent monomials.
EDDIITT: pretty good approximation ( for $x>4,$ say) with $$ \color{blue}{ h(x) \approx x^{\sqrt 2} + \frac{x^{ \left(\sqrt 2 - 1 \right)}}{\sqrt 2} + (1 - \sqrt 2 ) }$$ The following things are true: as long as you just want $C^1$ with no hope of extending to the complex numbers, you can do it, this is a theorem in the KCG book, I put some pages as a pdf HERE . I have been curious on one technical point for four years now, smoothness of the real restriction of Ecalle's solution at the fixpoint itself, and just wrote to Prof. Ger, maybe he will write back.
Meanwhile, see https://mathoverflow.net/questions/45608/formal-power-series-convergence and the correct answer at https://mathoverflow.net/questions/45608/formal-power-series-convergence/46765#46765
Your function $x^2 + x$ also has derivative $1$ at the fixpoint $0,$ but is not a Mobius transformation. What this means is that there is a solution for $x \geq 0$ that is real analytic for $x > 0$ and can therefore be extended to a holomorphic function in an open set containing the strictly positive real axis. The technique for doing all this is due to Jean Ecalle at Orsay, about 1973. The specific steps are in the KCG book, especially pages 346-347 and 351-352. All steps are done with formal power series, abbreviated FPS in the book.
Anyway, goes like this: define $$ \color{green}{f(x) = \frac{\sqrt{1 + 4 x} - 1}{2} = \frac{2x}{1 + \sqrt{1 + 4 x}} } $$ for $x > -1/4.$ We need to use use rather than the original $x^2 + x$ because we need convergence in the iteration steps.
$$ \color{magenta}{f = x - x^2 + 2 x^3 - 5 x^4 + 14 x^5 - 42 x^6 + 132 x^7 - 429 x^8 + 1430 x^9 - 4862 x^{10} + 16796 x^{11} - 58786 x^{12} + 208012 x^{13} - 742900 x^{14} + 2674440 x^{15} + O(x^{16})} $$
$$ $$
$$ \color{magenta}{ \frac{d f}{dx} = 1 - 2 x + 6 x^2 - 20 x^3 + 70 x^4 - 252 x^5 + 924 x^6 - 3432 x^7 + 12870 x^8 - 48620 x^9 + 184756 x^{10} - 705432 x^{11} + 2704156 x^{12} - 10400600 x^{13} + 40116600 x^{14} - 155117520 x^{15} + O(x^{16}) }$$
Find several terms in the formal power series for $\lambda(x)$ that solves $$ \lambda(f(x)) = f'(x) \lambda(x), $$ or $$ \lambda \left( \frac{\sqrt{1 + 4 x} - 1}{2} \right) = \frac{\lambda(x)}{ \sqrt{1 + 4 x}}, $$ where the power series for $\lambda(x)$ is required to begin with the first term in the power series of $f(x)$ after the initial $x.$ To emphasize, you find the FPS's as above and perform this step with those FPS; I gradually extend the series for $\lambda,$ one coefficient at a time. $$ \color{magenta}{\lambda = - x^2 + x^3 - \frac{3 x^4}{2} + \frac{8 x^5}{3} - \frac{31 x^6}{6} + \frac{157 x^7}{15} - \frac{649 x^8}{30} + \frac{9427 x^9}{210} - \frac{19423 x^{10}}{210} + \frac{6576 x^{11}}{35} - \frac{2627 x^{12}}{7} + \frac{853627 x^{13}}{1155} - \frac{ 2007055 x^{14}}{ 1386} + \frac{3682190 x^{15}}{ 1287} + O(x^{16}))}$$
Next, write several terms for the reciprocal of the series, and use those in $$ \frac{1}{\lambda(x)} = \frac{d \alpha(x)}{dx},$$
$$\color{magenta}{ \frac{d \alpha}{dx} = \frac{-1}{x^2} - \frac{1}{x} + \frac{1}{2} - \frac{2x}{3} + \frac{13x^2}{12} - \frac{113x^3}{60}+ \frac{1187x^4}{360} - \frac{1754x^5}{315} + \frac{14569x^6}{1680} - \frac{176017x^7}{15120} + \frac{ 1745717x^8}{151200} - \frac{ 176434x^9}{51975} - \frac{ 147635381x^{10}}{9979200} + \frac{ 3238110769x^{11}}{129729600} + O(x^{12})}$$
Now, formally integrate to find a short series for $\alpha(x)$ that usually includes a single logarithms term and begins with a few negative powers of $x,$ so it is a logarithm term plus a Laurent expansion. This function $\alpha(x)$ satisfies $$ \alpha(f(x)) = \alpha(x) + 1. $$
$$ \color{magenta}{ \alpha = \frac{1}{x} - \log x + \frac{x}{2} - \frac{x^2}{3} + \frac{13 x^3}{36} - \frac{113 x^4}{240} + \frac{1187x^5}{1800} - \frac{877x^6}{945} + \frac{14569x^7}{11760} - \frac{176017x^8}{120960} + \frac{1745717x^9}{1360800} - \frac{88217x^{10}}{259875} + O(x^{11})}$$
To actually calculate $\alpha(x)$ for some real number $x > 0,$ define $$x_0 = x, x_1 = f(x), \; x_2 = f(x_1), \; \ldots \; x_{n+1} = f(x_n). $$ From the dfining equation for $\alpha,$ we know that $$ \alpha(x_n) - n = \alpha(x). $$ Which is very good, because $x_n$ slowly approaches $0,$ and we can find $\alpha(x)$ to arbitrary accuracy with $$ \color{magenta}{ \alpha(x) = \lim_{n \rightarrow \infty} \alpha(x_n) - n}, $$ where we are using our peculiar Laurent expansion plus logarithm term in the right hand side. We need a second numerical step, which is to have available $\alpha^{-1}(x).$ I did that with ordinary bisection, slow but reliable.
Finally, you were really interested in $$ f^{-1}(x) = x^2 + x. $$ By simple substitution, we have $$ \alpha(f^{-1}(x)) = \alpha(x) - 1. $$
Define $$\color{blue}{ h(x) = \alpha^{-1}\left( \alpha(x) - \frac{1}{2} \right)}, $$ so that $$ \alpha(h(x)) = \alpha(x) - \frac{1}{2}. $$ Then $$ h(h(x)) = \alpha^{-1}\left( \alpha(h(x)) - \frac{1}{2} \right), $$ $$ h(h(x)) = \alpha^{-1}\left( \left( \alpha(x) - \frac{1}{2} \right) - \frac{1}{2} \right) = \alpha^{-1}\left( \alpha(x) - 1 \right) = \alpha^{-1}\left( \alpha(f^{-1}(x)) \right), $$ $$ \color{blue}{ h(h(x)) = f^{-1}(x) = x^2 + x}. $$
This really is the right way to do this. It is just lots of work.
Alright, used gp-pari, the combined Laurent series with log is
If it seems convenient one may include a constant term, in the end it changes nothing.
EDIT, Friday August 29. Quicker than I expected, largely because I still had the C++ program for the sine problem and just had a few changes, all the extra tinkering was in numerical stuff, accuracy demands and so on. The half iterate is called $h(x),$ next column $h(h(x))$ which came out very well, error $h(h(x)) - x - x^2$ in final column.