What do you do to find the power series for an inverse relationship such as for $y$ in $y + \sin(y) = x$? I'm not sure where to begin.

(Similarly, the Lambert $W$ function has such a power series for $y$ in $y = W(x)$ ($y e^y = x$) which allows computation.)

  • $\begingroup$ Lagrange's inversion theorem may help. $\endgroup$ – lhf Apr 23 '12 at 14:33
  • $\begingroup$ The smallest zeros of the derivative of $f(x)=x+\sin x$ have absolute value $\pi$. Therefore the radius of convergence of the resulting power series cannot exceed $|f(\pi)|=\pi$. But is it equal to $\pi$? $\endgroup$ – Jyrki Lahtonen Apr 23 '12 at 16:27

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}). $$

  • $\begingroup$ Can I define this power series using summation notation? $\endgroup$ – jnm2 Apr 23 '12 at 14:04
  • $\begingroup$ @DavidMitra: It is the analytical solution and the simplest one to explain the technique. Of course, I could take another solution of the initial equation $y_0$ and work it out similarly. $\endgroup$ – Jon Apr 23 '12 at 14:04
  • $\begingroup$ @Jon oops, sorry. I wasn't thinking clearly.. $\endgroup$ – David Mitra Apr 23 '12 at 14:06
  • 1
    $\begingroup$ @jnm2: This particular case seems to have the simple expression $$y(x)=\sum_{n=0}^\infty\frac{x^{2n+1}}{2^{2n+1}}.$$ $\endgroup$ – Jon Apr 23 '12 at 14:07
  • 3
    $\begingroup$ Also sprach Mathematica: $$y(x)=\frac x2+\frac{x^3}{96}+\frac{x^5}{1920}+\frac{43x^7}{1290240}+\frac{223x^9}{92897280}+O(x^{11}).$$ But Jon described the process. He had my +1 already. Gotta go home, now. $\endgroup$ – Jyrki Lahtonen Apr 23 '12 at 15: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 +
*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)))

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

  • $\begingroup$ Shouldn't that be 'you had $n$ coefficients of the inverse, and now you have $2n$'? Convergence is very good, but it's not quite tower-of-exponentials good. :-) $\endgroup$ – Steven Stadnicki Apr 23 '12 at 16:16
  • $\begingroup$ You're right, of course ! $\endgroup$ – Lierre Apr 23 '12 at 16:19
  • $\begingroup$ +1 Very impressive convergence, indeed! But why did your example solve for the inverse of $t+\sinh t$ as opposed to $t+\sin t$? It doesn't make much difference, of course. And I really needed to spend quite some time reading this before the alternating signs began to bug me :-) $\endgroup$ – Jyrki Lahtonen Apr 23 '12 at 17:02
  • $\begingroup$ Oh... You're right too... (In fact, I admit I don't know, in sage, how to get a taylor expansion in a power series ring.) $\endgroup$ – Lierre Apr 23 '12 at 17:12

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 $$


lhf already brought up Lagrangian inversion (a special case of the more general Lagrange-Bürmann series):

$$f^{(-1)}(x)=\sum_{k=0}^\infty \frac{x^{k+1}}{(k+1)!} \left(\left.\frac{\mathrm d^k}{\mathrm dt^k}\left(\frac{t}{f(t)}\right)^{k+1}\right|_{t=0}\right)$$

There are a number of nice series reversion algorithms that make use of coefficients from the series to be inverted. Here is a Mathematica implementation of an algorithm due to Henry Thacher (also used here):

a = Rest[CoefficientList[Series[(x + Sin[x])/2, {x, 0, 20}], x]];
n = Length[a];
      c[i, j + 1] = Sum[c[k, 1]c[i - k, j], {k, 1, i - j}];
      , {j, i - 1, 1, -1}];
    c[i, 1] = Boole[i == 1] - Sum[a[[j]] c[i, j], {j, 2, i}]
    , {i, n}];
Table[c[i, 1]/2^i, {i, n}]

{1/2, 0, 1/96, 0, 1/1920, 0, 43/1290240, 0, 223/92897280, 0, 
60623/326998425600, 0, 764783/51011754393600, 0, 
107351407/85699747381248000, 0, 2499928867/23310331287699456000, 0, 

Thacher's method expects the series to be inverted to take the form $x+\cdots$ (and the series for $x+\sin\,x$ starts out $2x+\cdots$), so some rescaling is necessary before feeding the series coefficients to the algorithm; the division by powers of $2$ at the end recovers the coefficients of the original function to be inverted.

There a lot more methods (e.g. Carleman matrices, Bell polynomials); search around for more information on them.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.