# Solution of equation with power series (perturbation)

So I need to use Mathematica to find the solution of $y=x- \epsilon \sin(2y)$ as a power series in terms of $\epsilon$. I'd assume I'd need to create an equation $f=x-y- \epsilon \sin(2y)$, then express $y(\epsilon)=\sum_n a_n \epsilon^n$, then input into series, but I can't seem to get it to work. Some help would be appreciated.

• Hello Lenny! Welcome to the Mathematica StackExchange community. Have you tried to do anything by yourself? I mean, is there any piece of code which you've tried to use?
– Rod
Jun 13, 2013 at 10:13
• I have. What I have tried so far: eq[eps_,y_]:=x-eps Sin[2y]; y[eps_]:=Sum[Subscript[a,n] eps^n,{n,0,20}]; Series[eq[eps,y[eps]],{eps,0,2}]; Jun 13, 2013 at 10:16
• Is eps any Mathematica function? Or is it just an argument of a function?
– Rod
Jun 13, 2013 at 10:19
• Hint: let $x=3\eta/2-\xi, y=\eta/2$, then your equation will reduce to Kepler's Equation. Jun 13, 2013 at 11:03
• It's just I think this type of equations are not as innocent as they look like, so instead of trying it with trivial methods, adopting some sophisticated results might be more meaningful. Anyway, you can try this method. Jun 13, 2013 at 11:33

Let define the equation to solve $f=x-y\epsilon\sin(2 x)\equiv 0$ and series expansion of $y$ in powers of $\varepsilon$.

f = x - y - \[Epsilon] Sin[2 y];
ord = 3;
y = x + Sum[a[n] \[Epsilon]^n, {n, ord}] + O[\[Epsilon]]^(ord + 1);


Then expand the equation and solve for any value of $\varepsilon$ parameter

eqs = Normal[Series[f, {\[Epsilon], 0, ord}]];
SolveAlways[eqs == 0, \[Epsilon]]
(* => {{a -> 2 (a - a^3 + a Cos[4 x]), a -> -2 a Cos[2 x], Sin[2 x] -> -a}} *)


If we are interested in the expansion coefficients $a_n$ in terms of $x$ then following code will refine the solution

sol = First@Solve[Equal@@@First@SolveAlways[eqs == 0, \[Epsilon]], a/@Range[ord]]
(* => {a -> -Sin[2 x], a -> 2 Cos[2 x] Sin[2 x],  a -> 2 Sin[2 x] (-1 - Cos[4 x] + Sin[2 x]^2)} *)


One should also do check

f/.sol//Simplify
(* => O[\[Epsilon]]^4 *)


The new in M12 function AsymptoticSolve can be used to find the perturbation expansions:

AsymptoticSolve[y == x - ϵ Sin[2 y], {y, x}, {ϵ, 0, 4}]


{{y -> x - ϵ Sin[2 x] + 2 ϵ^2 Cos[2 x] Sin[2 x] - 2 ϵ^3 (2 Cos[2 x]^2 Sin[2 x] - Sin[2 x]^3) + 8/3 ϵ^4 (3 Cos[2 x]^3 Sin[2 x] - 5 Cos[2 x] Sin[2 x]^3)}}

First, one can express $\varepsilon$ as a function of $y$:

Clear[eps]
eps[y_] = ε /. Solve[x == y - ε Sin[2 y], ε][]


$$-(x-y) \csc (2 y).$$

Then, expand into series around $y=x$, since it is a solution for $\varepsilon=0$. Applying InverseSeries after that gives:

res = InverseSeries[Series[eps[y], {y, x, 3}]] /. y -> ε // FullSimplify


$$x+\varepsilon \sin (2 x)+\varepsilon ^2 \sin (4 x)+\varepsilon ^3 \sin (2 x) (3 \cos (4 x)+1)+O\left(\varepsilon ^4\right).$$

Check:

Series[y - ε Sin[2 y] /. y -> res, {ε, 0, 3}] // FullSimplify


$$x+O\left(\varepsilon ^4\right).$$