12
$\begingroup$

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

  • 0
    [Lagrange's inversion theorem](http://en.wikipedia.org/wiki/Lagrange_inversion_theorem) may help.2012-04-23
  • 0
    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$?2012-04-23

4 Answers 4

10

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

  • 0
    Can I define this power series using summation notation?2012-04-23
  • 0
    @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.2012-04-23
  • 0
    @Jon oops, sorry. I wasn't thinking clearly..2012-04-23
  • 1
    @jnm2: This particular case seems to have the simple expression $$y(x)=\sum_{n=0}^\infty\frac{x^{2n+1}}{2^{2n+1}}.$$2012-04-23
  • 0
    That isn't right though. It converges to something way off.2012-04-23
  • 0
    @Jon, a nice exposition, but does the series really come to that?? That's a geometric series, so the sum is a rational function, which is too good to be true, surely??2012-04-23
  • 0
    Hmm. May be you forget to divide the $x^n$ term with $n!$ ??2012-04-23
  • 0
    @JyrkiLahtonen: I think you are right. I will fix this. The final sum should be $$y(x)=\sum_{n=0}^\infty\frac{x^{2n+1}}{2^{2n+1}(2n+1)!}.$$ Thank you very much!2012-04-23
  • 0
    That's converging to the wrong function too... I was kind of wondering how to find the summation, not just what this particular summation is.2012-04-23
  • 0
    @jnm2: What are you asking for? Please, explain better your thoughts. Do you have anything known that we don't?2012-04-23
  • 0
    No, but I have the graph of $y + \sin(y) = x$ and these power series just aren't matching up.2012-04-23
  • 0
    I'm trying something similar (substituting $y(x)$ to the series of sine), and I get a better match with $y=x/2+x^3/96+O(x^5)$. Can you check the calculation for $y'''(0)$, please?2012-04-23
  • 0
    I'm getting $x^3/16$. See [W|A](http://www.wolframalpha.com/input/?i=third+derivative+of+y%28x%29+%2B+sin%28y%28x%29%29+%3D+x) for third derivative and substitute $y(x) = 0$, $y'(x) = \frac1{2}$ and $y''(x) = 0$. @Jon, I'll upvote if you make sure your work is correct.2012-04-23
  • 0
    If $y=x/2+x^3/16$, I get $$y+(y-y^3/6)=x+5x^3/48+ O(x^5).$$ If $y=x/2+x^3/96$, I get $$y+(y-y^3/6)=x+O(x^5).$$ I think we should get something like that latter one??2012-04-23
  • 0
    @jnm2, Mathematica gave me the equation $$-y'(0)^3+2y'''(0)=0$$ after substituting $y(0)=0$. With $y'(0)=1/2$, this gives $y'''(0)=1/16$, so the series should begin like $$y(x)=\frac x2+\frac1{3!}\frac{x^3}{16}+O(x^5),$$ right?2012-04-23
  • 0
    @JyrkiLahtonen, Oh, you're right. I forgot the factorial. Well hey, if you have Mathematica, do this: http://reference.wolfram.com/mathematica/ref/InverseSeries.html2012-04-23
  • 0
    @jnm2, I'll be! I totally didn't know about that. I just calculated that $y^{(5)}(0)=1/16$, too, but I'll check this out.2012-04-23
  • 3
    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.2012-04-23
  • 0
    Thank you a lot, people! Now it is all fine.2012-04-23
9

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. = 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...

  • 0
    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. :-)2012-04-23
  • 0
    You're right, of course !2012-04-23
  • 0
    +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 :-)2012-04-23
  • 0
    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.)2012-04-23
6

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

5

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]; Do[     Do[       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,  596767688063/63777066403145711616000} 

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.