# 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 '13 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}]; – Lenny Jun 13 '13 at 10:16
• Is eps any Mathematica function? Or is it just an argument of a function? – Rod Jun 13 '13 at 10:19
• Hint: let $x=3\eta/2-\xi, y=\eta/2$, then your equation will reduce to Kepler's Equation. – Silvia Jun 13 '13 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. – Silvia Jun 13 '13 at 11:33

## 3 Answers

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

in agreement with mmal's answer.

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