Finding the power series for $y$ where $y + \sin(y) = x$

You should compute the $n$th derivative of the given expression. So, $$ y(x)+\sin(y(x))=x $$ and then $y(0)=0$. The first derivative will give $$ y'(x)(1+\cos(y(x))=1 $$ and so $$ y'(0)=\frac1{2}. $$ Proceeding in a similar way you will get $y''(0)=0$ and $y'''(0)=\frac{1}{16}$ and eventually you will have $$ y(x)=\frac1{2}x+\frac1{16}\frac{x^3}{3!}+\frac1{16}\frac{x^5}{5!}+\frac{43}{256}\frac{x^7}{7!}+\frac{223}{256}\frac{x^9}{9!}+O(x^{11}). $$

This is an instance of the problem of finding a compositional inverse for a power series. Let $F$ in $k[[t]]$ a formal power series, with $F(0) = 0$ and $F'(0) \neq 0$. How to compute $G\in k[[t]]$ such that $F(G) = t$ ? In your example, $F$ is $t + \sin t$ as a formal power series.

A good answer is the Newton's iteration. It both proof the existence of the compositional inverse and gives an efficient way to compute it.

I — Description of the algorithm

Let $n>0$ and assume that you already have $G\in k[[t]]$ such that $F(G) = t \mod t^n$. Define :

• $E = F(G) - t$, the error term ;
• $G_1 = G - E \cdot F'(G)^{-1}$, the new approximation.

Note that $G = G_1 \mod t^n$ since $E = 0 \mod t^n$. What is $F(G_1)$ ? Compute it using taylor expansion $$ \begin{aligned} F(G_1) &= F( G + (G_1 - G) )\\ &= F(G) + (G_1 - G) \cdot F'(G) + \mathcal O((G_1 - G)^2)\\ &= F(G) - E + \mathcal{O}(E^2) \\ &= 0 \mod t^{2n} \end{aligned}$$

That's incredible ! You had $n$ coefficients of the compositional inverse, and now you have $2n$, with a single iteration !

Lagrange inversion theorem is good is you want only one coefficient of the inverse, without having to compute the previous one. If you want the whole expansion, then Newton is way better.

II — Example

Define $$F = t + \sum_n \frac{(-1)^n t^{2n+1}}{(2n+1)!}, $$ and $G_0 = 0$.

The first iteration gives

1/2*t + O(t^2)

The second gives

1/2*t + 1/96*t^3 + O(t^4)

The fifth gives

1/2*t + 1/96*t^3 + 1/1920*t^5 + 43/1290240*t^7 + 223/92897280*t^9 +
60623/326998425600*t^11 + 764783/51011754393600*t^13 +
107351407/85699747381248000*t^15 + 2499928867/23310331287699456000*t^17
+ 596767688063/63777066403145711616000*t^19 +
22200786516383/26786367889321198878720000*t^21 +
64470807442488761/867449737727777704488468480000*t^23 +
3504534741776035061/520469842636666622693081088000000*t^25 +
3597207408242668198973/5845917272495039506088686780416000000*t^27 +
268918457620309807441853/4746884825265972078944013665697792000000*t^29 +
185388032403184965693274807/35316823099978832267343461672791572480000000\
*t^31 + O(t^32)

III — The code

To proceed to the computations, I used Sage, and here is the code :

R.<t> = PowerSeriesRing(QQ)

f = t + sum([ (-1)^n*t^(2*n+1)/factorial(2*n + 1) for n in range(20) ]) + O(t^41)

def newton_it( f, g ) :
  g = parent(g)(g.polynomial()).add_bigoh(2*g.prec())
  return g + (t - f(g))/derivative(f)(g)

newton_it(f, 0 + O(t))
newton_it(f, newton_it(f, 0 + O(t)))
etc.

No more than two lines for the actual newton iteration...

The Bell/Carleman-method, which J.M. mentions is my favorite method for tasks like this, which can be done by paper&pen for a truncation to coefficients at low index . First step is to list sufficiently many formal powers of the powerseries in question. So if $\small f(y)=y + \sin (y) = 2y-y^3/3!+y^5/5!- \ldots + \ldots $

then we list the first few formal powers of f(y): $$\small \begin{array} {lll} f(y)^0&=& 1 + O(x^8)\\ f(y)^1&=& 2*x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + O(x^8)\\ f(y)^2&=& 4*x^2 - 2/3*x^4 + 11/180*x^6 + O(x^8)\\ f(y)^3&=& 8*x^3 - 2*x^5 + 4/15*x^7 + O(x^8)\\ f(y)^4&=& 16*x^4 - 16/3*x^6 + O(x^8) \end{array} $$ then the Carleman-matrix looks like $$\small M_{f(y)}= \begin{bmatrix} 1&.&.&.&.&.&.&.\\ .&2&.&-1/6&.&1/120&.&-1/5040\\ .&.&4&.&-2/3&.&11/180&.\\ .&.&.&8&.&-2&.&4/15\\ .&.&.&.&16&.&-16/3&.\\ .&.&.&.&.&32&.&-40/3\\ .&.&.&.&.&.&64&.\\ .&.&.&.&.&.&.&128 \end{bmatrix} $$ where the coefficients of the above powers of f(y) are filled rowwise into the matrix. This can be done using paper&pen for the first few rows/columns.
Because the matrix is triangular it is then again easy to invert it, by paper&pen at least for some leading terms and this gives: $$\small M_{f^{[-1]}(y)}= \begin{bmatrix} 1&.&.&.&.&.&.&.\\ .&1/2&.&1/96&.&1/1920&.&43/1290240\\ .&.&1/4&.&1/96&.&29/46080&.\\ .&.&.&1/8&.&1/128&.&17/30720 \end{bmatrix} $$ We need do this even only to the second row, because this contains the coefficients for the formal powerseries for the (compositional) inverse
$$\small f(x)=f^{[-1]}(y) = 1/2 x + 1/96 x^3 + 1/1920 x^5 + \ldots $$