Power series approximation for intractable function? (using Solve)

Note: I am using 2*x+cos(x) as an example below. I'm asking this question for any function where Solve[] doesn't yield a closed form result. My specific interest at the moment is getting Mathematica to spit out the (well known) power series expansion solution of Kepler's Equation of the Center, just as a test.

The function 2*x+cos(x) is monotonic increasing (derivative is 2-sin(x)), and thus invertible. There's no simple form for the inverse:


Solve[2*x+Cos[x] == y, x]

Solve::nsmet: This system cannot be solved with the methods available to Solve.


Question: can I get a power series approximation to the solution? Either using Mathematica builtin functions or some sort of voodoo?

• For $y\ll1$ you can expand $\cos(x)=1-x^2/2$ and solve the quadratic equation $1+(2x-y)-x^2/2=0$. For $y\rightarrow\infty$ the solution is $2x\approx y$. Power series can be likewise used to improve it. – yarchik Nov 19 '17 at 14:53
• @yarchik Thanks! Actually, I was asking for functions in general and just using that as an example. I did consider expanding Cos[x] as a power series and solving the polynomial, but that doesn't really give me a power series. – user1722 Nov 19 '17 at 14:57
• Well, formally the solution to your problem is easy: make a series expansion of the l.h.s. and invert the series ('InverseSeries'). However, this will not be a good approach in your case. Much better is to start from a good guess (to guarantee the convergence) and apply (analytically) the Newton's method. – yarchik Nov 19 '17 at 15:04

Use InverseSeries or better yet use InverseFunction

Clear[f, g, gs]

f[x_] = 2 x + Cos[x];

g[x_] := InverseFunction[f][x]

gs[x_] = InverseSeries[Series[f[x], {x, 0, 10}]] // Normal;

Prepend[Table[
{x, f[g[x]], g[f[x]], f[gs[x]], gs[f[x]]},
{x, -2., 2., 0.25}],
Style[#, 14, Bold] & /@
{"x", "f[g[x]]", "g[f[x]]", "f[gs[x]]",
"gs[f[x]]"}] // Grid[#, Frame -> All] &


As with any series expansion,InverseSeries only works well close to the expansion point.

• Didn't know about InverseSeries, thanks! I'm using the results externally, so InverseFunction wouldn't work. – user1722 Nov 19 '17 at 15:27
• @barrycarter - Thanks for the accept. I recommend that you always look at the See Also section of the documentation page (for Series in this case). – Bob Hanlon Nov 19 '17 at 15:33

The new in M12 function AsymptoticSolve can be used to find the perturbation expansions. For example, here's the solution for x when y is expanded about 1 (as in Bob's answer):

AsymptoticSolve[2 x + Cos[x] == y, {x, 0}, {y, 1, 4}]


{{x -> 1/2 (-1 + y) + 1/16 (-1 + y)^2 + 1/64 (-1 + y)^3 + (11 (-1 + y)^4)/ 3072}}

Alternatively, we can solve the equation for when y is 0:

root = x /. First @ Solve[2 x + Cos[x] == 0, x, Reals];
root //InputForm


Root[{Cos[#1] + 2*#1 & , -0.4501836112948735730376641555814704818420.60205991050527}]

Since the above root object is hard to process, I'll introduce a symbol γ that is equal to the root but displays better:

NumericQ[γ] = True;
N[γ, _] := root;


Then, using AsymptoticSolve again:

AsymptoticSolve[2 x + Cos[x] == y, {x, γ}, {y, 0, 3}]
`

{{x -> γ + y/(2 - Sin[γ]) - (y^2 Cos[γ])/( 2 (-2 + Sin[γ])^3) + ( y^3 (-3 Cos[γ]^2 + 2 Sin[γ] - Sin[γ]^2))/( 6 (-2 + Sin[γ])^5)}}