# Lindstedt Poincare Example - Rayleigh's Equation

For the Lindstedt-Poincare method we introduce a dimensionless time $\tau$

$$
\begin{align*}
    \tau &= \omega t \\
    \dfrac{d}{dt} &= \omega \dfrac{d}{d\tau}\\
    \dfrac{d^2}{dt^2} &= \omega^2 \dfrac{^2}{d\tau^2}
\end{align*}
$$

We will use Rayleigh's Equation as an example

$$
\begin{gather*}
    \ddot{x} + \epsilon \left( \dfrac{1}{3}\dot{x}^3 - \dot{x} \right) + x = 0
\end{gather*}
$$

with initial conditions of 

$$
\begin{gather*}
    x(0) = a\\
    \dot{x}(0) = 0
\end{gather*}
$$

and the system will be investigated up to only the second order.


```python
import sympy as sp
from sympy.simplify.fu import TR0, TR5, TR7, TR8, TR11
from math import factorial
```


```python
N = 3  # first order

# Define the symbolic parameters
epsilon = sp.symbols('epsilon')
omega = sp.symbols('omega')
omega_i = sp.symbols('omega_(0:' + str(N) + ')')
tau = sp.symbols('tau')


x = sp.Function('x')(tau)
xdot = sp.Derivative(x, tau)
xddot = sp.Derivative(xdot, tau)

x0 = sp.Function('x_0')(tau)
x1 = sp.Function('x_1')(tau)
x2 = sp.Function('x_2')(tau)
x3 = sp.Function('x_3')(tau)

# EOM
EOM = omega**2 * xddot + x + epsilon * (omega**3 * xdot**3/3 - omega * xdot)
EOM
```




$\displaystyle \epsilon \left(\frac{\omega^{3} \left(\frac{d}{d \tau} x{\left(\tau \right)}\right)^{3}}{3} - \omega \frac{d}{d \tau} x{\left(\tau \right)}\right) + \omega^{2} \frac{d^{2}}{d \tau^{2}} x{\left(\tau \right)} + x{\left(\tau \right)}$




```python
x_i = (x0, x1, x2, x3)
x_e = sum([epsilon**i * x_i[i] for i in range(N+1)])
x_e
```




$\displaystyle \epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon \operatorname{x_{1}}{\left(\tau \right)} + \operatorname{x_{0}}{\left(\tau \right)}$




```python
_omega = sum([epsilon**i * omega_i[i] for i in range(N)])
_omega
```




$\displaystyle \epsilon^{2} \omega_{2} + \epsilon \omega_{1} + \omega_{0}$




```python
# Substitute these into the EOM
EOM = EOM.subs([
    (x, x_e), (omega, _omega)
])
EOM
```




$\displaystyle \epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon \left(\frac{\left(\epsilon^{2} \omega_{2} + \epsilon \omega_{1} + \omega_{0}\right)^{3} \left(\frac{\partial}{\partial \tau} \left(\epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon \operatorname{x_{1}}{\left(\tau \right)} + \operatorname{x_{0}}{\left(\tau \right)}\right)\right)^{3}}{3} - \left(\epsilon^{2} \omega_{2} + \epsilon \omega_{1} + \omega_{0}\right) \frac{\partial}{\partial \tau} \left(\epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon \operatorname{x_{1}}{\left(\tau \right)} + \operatorname{x_{0}}{\left(\tau \right)}\right)\right) + \epsilon \operatorname{x_{1}}{\left(\tau \right)} + \left(\epsilon^{2} \omega_{2} + \epsilon \omega_{1} + \omega_{0}\right)^{2} \frac{\partial^{2}}{\partial \tau^{2}} \left(\epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon \operatorname{x_{1}}{\left(\tau \right)} + \operatorname{x_{0}}{\left(\tau \right)}\right) + \operatorname{x_{0}}{\left(\tau \right)}$




```python
EOM = sp.expand(EOM).doit()
EOM = sp.expand(EOM)
EOM
```




$\displaystyle \frac{\epsilon^{16} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3}}{3} + \epsilon^{15} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + \epsilon^{15} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + \epsilon^{14} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + \epsilon^{14} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{14} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + \epsilon^{14} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + \epsilon^{14} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{13} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{13} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + \frac{\epsilon^{13} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3}}{3} + 3 \epsilon^{13} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{13} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{13} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{13} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 2 \epsilon^{13} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \frac{\epsilon^{13} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3}}{3} + \epsilon^{12} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + \epsilon^{12} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + 6 \epsilon^{12} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{12} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{12} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{12} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{12} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{12} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{12} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{12} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{12} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + 2 \epsilon^{12} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{12} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{12} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + \epsilon^{11} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{11} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{11} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{11} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{11} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{11} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{11} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{11} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + \epsilon^{11} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + \epsilon^{11} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{11} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{11} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{11} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + 6 \epsilon^{11} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{11} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{11} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 2 \epsilon^{11} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{11} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + \epsilon^{11} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \frac{\epsilon^{10} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{3}}{3} + 3 \epsilon^{10} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{10} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{10} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{10} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{10} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{10} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 12 \epsilon^{10} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{10} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + 6 \epsilon^{10} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{10} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{10} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + \epsilon^{10} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 2 \epsilon^{10} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \frac{\epsilon^{10} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3}}{3} + 6 \epsilon^{10} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{10} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{10} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{10} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{10} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{10} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{10} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{10} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \frac{\epsilon^{10} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3}}{3} + \epsilon^{9} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{9} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{9} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{9} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{9} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{9} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{9} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{9} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{9} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + 12 \epsilon^{9} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{9} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{9} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{9} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{9} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{9} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 2 \epsilon^{9} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{9} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{9} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{9} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{9} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{9} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{9} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{9} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{9} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + \epsilon^{9} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{9} \omega_{2}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + \epsilon^{8} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + \epsilon^{8} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{8} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{8} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3} + 6 \epsilon^{8} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{8} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 12 \epsilon^{8} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{8} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{8} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{8} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{8} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + 2 \epsilon^{8} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{8} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + \epsilon^{8} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{8} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{8} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{8} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{8} \omega_{1} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + \epsilon^{8} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{7} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)}\right)^{2} + 2 \epsilon^{7} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \frac{\epsilon^{7} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{3}}{3} + 6 \epsilon^{7} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{7} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{7} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{7} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{7} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{7} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 6 \epsilon^{7} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{7} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{7} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 6 \epsilon^{7} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 12 \epsilon^{7} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 2 \epsilon^{7} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{7} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{7} \omega_{0} \omega_{2}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + \epsilon^{7} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{7} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \frac{\epsilon^{7} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3}}{3} + 3 \epsilon^{7} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{7} \omega_{1}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{7} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \frac{\epsilon^{7} \omega_{2}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3}}{3} + \epsilon^{7} \omega_{2}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{6} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{6} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{6} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{6} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 3 \epsilon^{6} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{6} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{6} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{6} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{6} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{6} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{6} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{6} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + 6 \epsilon^{6} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 6 \epsilon^{6} \omega_{0} \omega_{1} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{6} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{6} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{6} \omega_{1}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{6} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{6} \omega_{1} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + 2 \epsilon^{6} \omega_{1} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{6} \omega_{2}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} - \epsilon^{6} \omega_{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{5} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{5} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)}\right)^{2} + \epsilon^{5} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{5} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 6 \epsilon^{5} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{5} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + 3 \epsilon^{5} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{5} \omega_{0}^{2} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{5} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{5} \omega_{0} \omega_{1}^{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 6 \epsilon^{5} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{5} \omega_{0} \omega_{2}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + 2 \epsilon^{5} \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{5} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{5} \omega_{1}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + \epsilon^{5} \omega_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{5} \omega_{1} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} - \epsilon^{5} \omega_{1} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{5} \omega_{2}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} - \epsilon^{5} \omega_{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{4} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{4} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \frac{\epsilon^{4} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{3}}{3} + 3 \epsilon^{4} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + 3 \epsilon^{4} \omega_{0}^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{4} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + 3 \epsilon^{4} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + 2 \epsilon^{4} \omega_{0} \omega_{1} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + 2 \epsilon^{4} \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + 2 \epsilon^{4} \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} - \epsilon^{4} \omega_{0} \frac{d}{d \tau} \operatorname{x_{3}}{\left(\tau \right)} + \frac{\epsilon^{4} \omega_{1}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3}}{3} + \epsilon^{4} \omega_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 2 \epsilon^{4} \omega_{1} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} - \epsilon^{4} \omega_{1} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{4} \omega_{2}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \epsilon^{4} \omega_{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{3} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{3} \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \epsilon^{3} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{3} \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + \epsilon^{3} \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{3} \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + 2 \epsilon^{3} \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 2 \epsilon^{3} \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} - \epsilon^{3} \omega_{0} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon^{3} \omega_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 2 \epsilon^{3} \omega_{1} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \epsilon^{3} \omega_{1} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} - \epsilon^{3} \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} + \epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{2} \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + \epsilon^{2} \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 2 \epsilon^{2} \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 2 \epsilon^{2} \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \epsilon^{2} \omega_{0} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{2} \omega_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \epsilon^{2} \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \frac{\epsilon \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3}}{3} + \epsilon \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 2 \epsilon \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \epsilon \omega_{0} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} + \epsilon \operatorname{x_{1}}{\left(\tau \right)} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} + \operatorname{x_{0}}{\left(\tau \right)}$




```python
# Collect the coefficients for the epsilons 
epsilon_Eq = sp.collect(EOM, epsilon, evaluate=False)
epsilon_0_Eq = sp.Eq(epsilon_Eq[epsilon**0], 0)
epsilon_0_Eq
```




$\displaystyle \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} + \operatorname{x_{0}}{\left(\tau \right)} = 0$




```python
epsilon_1_Eq = sp.Eq(epsilon_Eq[epsilon], 0)
epsilon_1_Eq
```




$\displaystyle \frac{\omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3}}{3} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 2 \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \omega_{0} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} + \operatorname{x_{1}}{\left(\tau \right)} = 0$




```python
epsilon_2_Eq = sp.Eq(epsilon_Eq[epsilon**2], 0)
epsilon_2_Eq
```




$\displaystyle \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 2 \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 2 \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \omega_{0} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \omega_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \omega_{1} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} + \operatorname{x_{2}}{\left(\tau \right)} = 0$




```python
epsilon_3_Eq = sp.Eq(epsilon_Eq[epsilon**3], 0)
epsilon_3_Eq
```




$\displaystyle \omega_{0}^{3} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \omega_{0}^{3} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} \left(\frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)}\right)^{2} + 3 \omega_{0}^{2} \omega_{1} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{2} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} + \omega_{0}^{2} \omega_{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + \omega_{0} \omega_{1}^{2} \left(\frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)}\right)^{3} + 2 \omega_{0} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 2 \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} - \omega_{0} \frac{d}{d \tau} \operatorname{x_{2}}{\left(\tau \right)} + \omega_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 2 \omega_{1} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{0}}{\left(\tau \right)} - \omega_{1} \frac{d}{d \tau} \operatorname{x_{1}}{\left(\tau \right)} - \omega_{2} \frac{d}{d \tau} \operatorname{x_{0}}{\left(\tau \right)} + \operatorname{x_{3}}{\left(\tau \right)} = 0$




```python
sp.dsolve(epsilon_0_Eq, x0) 
```




$\displaystyle \operatorname{x_{0}}{\left(\tau \right)} = C_{1} e^{- \frac{i \tau}{\omega_{0}}} + C_{2} e^{\frac{i \tau}{\omega_{0}}}$



Not as convenient as working with the polar form


```python
a = sp.symbols('a')
beta = sp.symbols('beta')
phi = tau + beta
x0_polar = a * sp.cos(tau + beta)
x0_polar
```




$\displaystyle a \cos{\left(\beta + \tau \right)}$



Update $\epsilon^1$ equation 


```python
epsilon_1_Eq = sp.expand(epsilon_1_Eq.subs(x0, x0_polar)).doit()
epsilon_1_Eq = sp.expand(TR8(sp.expand(TR7(TR5((epsilon_1_Eq))))))
epsilon_1_Eq
```




$\displaystyle - \frac{a^{3} \omega_{0}^{3} \sin{\left(\beta + \tau \right)}}{4} + \frac{a^{3} \omega_{0}^{3} \sin{\left(3 \beta + 3 \tau \right)}}{12} - 2 a \omega_{0} \omega_{1} \cos{\left(\beta + \tau \right)} + a \omega_{0} \sin{\left(\beta + \tau \right)} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + \operatorname{x_{1}}{\left(\tau \right)} = 0$



In the equation above we can see that the $\cos\phi$ terms are secular which causes the approximation fail due to resonance. Therefore, we must set those terms to go to zero.

$\omega_1$ = 0

$a = \pm \dfrac{2}{\omega_0}$

and let us pick $a = 2/\omega_0$


```python
omega_1_new = 0
```


```python
epsilon_1_Eq = epsilon_1_Eq.subs([(omega_i[1], omega_1_new), (a, 2/omega_i[0])])
epsilon_1_Eq = sp.simplify(epsilon_1_Eq)
epsilon_1_Eq
```




$\displaystyle \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + \operatorname{x_{1}}{\left(\tau \right)} + \frac{2 \sin{\left(3 \beta + 3 \tau \right)}}{3} = 0$




```python
temp = sp.dsolve(epsilon_1_Eq)
temp
```




$\displaystyle \operatorname{x_{1}}{\left(\tau \right)} = C_{1} e^{- \frac{i \tau}{\omega_{0}}} + C_{2} e^{\frac{i \tau}{\omega_{0}}} + \frac{2 \sin{\left(3 \beta + 3 \tau \right)}}{3 \left(9 \omega_{0}^{2} - 1\right)}$



We only care about the particular solution


```python
x1_p = 2 * sp.sin(3 * beta + 3 * tau) / 3 / (9 * omega_i[0]**2 - 1)
x1_p
```




$\displaystyle \frac{2 \sin{\left(3 \beta + 3 \tau \right)}}{3 \left(9 \omega_{0}^{2} - 1\right)}$




```python
epsilon_2_Eq = epsilon_2_Eq.subs([
    (omega_i[1], omega_1_new), (x0, x0_polar), (x1, x1_p)
]).doit()
epsilon_2_Eq = sp.expand(TR8(sp.expand(TR7(TR5((epsilon_2_Eq))))))
epsilon_2_Eq
```




$\displaystyle - \frac{a^{2} \omega_{0}^{3} \cos{\left(\beta + \tau \right)}}{18 \omega_{0}^{2} - 2} - \frac{a^{2} \omega_{0}^{3} \cos{\left(5 \beta + 5 \tau \right)}}{18 \omega_{0}^{2} - 2} + \frac{a^{2} \omega_{0}^{3} \cos{\left(3 \beta + 3 \tau \right)}}{9 \omega_{0}^{2} - 1} - 2 a \omega_{0} \omega_{2} \cos{\left(\beta + \tau \right)} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} - \frac{2 \omega_{0} \cos{\left(3 \beta + 3 \tau \right)}}{9 \omega_{0}^{2} - 1} + \operatorname{x_{2}}{\left(\tau \right)} = 0$



Remove the secular term by setting

$$
\omega_2 = -\dfrac{a\omega_0^2}{36\omega_0^2 - 4}
$$


```python
omega_2_new = -a*omega_i[0]**2 / (36*omega_i[0]**2 - 4)
```


```python
epsilon_2_Eq = epsilon_2_Eq.subs([(omega_i[2], omega_2_new), (a, 2/omega_i[0])])
epsilon_2_Eq = sp.simplify(epsilon_2_Eq)
epsilon_2_Eq
```




$\displaystyle \frac{9 \omega_{0}^{4} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 9 \omega_{0}^{2} \operatorname{x_{2}}{\left(\tau \right)} - \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{2}}{\left(\tau \right)} + 2 \omega_{0} \cos{\left(3 \beta + 3 \tau \right)} - 2 \omega_{0} \cos{\left(5 \beta + 5 \tau \right)} - \operatorname{x_{2}}{\left(\tau \right)}}{9 \omega_{0}^{2} - 1} = 0$




```python
temp = sp.dsolve(epsilon_2_Eq)
temp
```




$\displaystyle \operatorname{x_{2}}{\left(\tau \right)} = C_{1} e^{- \frac{i \tau}{\omega_{0}}} + C_{2} e^{\frac{i \tau}{\omega_{0}}} - \frac{2 \omega_{0} \cos{\left(5 \beta + 5 \tau \right)}}{225 \omega_{0}^{4} - 34 \omega_{0}^{2} + 1} + \frac{2 \omega_{0} \cos{\left(3 \beta + 3 \tau \right)}}{81 \omega_{0}^{4} - 18 \omega_{0}^{2} + 1}$




```python
x2_p = (-2 * omega_i[0] * sp.cos(5*beta + 5*tau) / (225 * omega_i[0]**4 - 34*omega_i[0]**2 + 1) 
        + 2 *omega_i[0]*sp.cos(3*beta + 3*tau) / (81*omega_i[0]**4 - 18*omega_i[0]**2 + 1))
x2_p
```




$\displaystyle - \frac{2 \omega_{0} \cos{\left(5 \beta + 5 \tau \right)}}{225 \omega_{0}^{4} - 34 \omega_{0}^{2} + 1} + \frac{2 \omega_{0} \cos{\left(3 \beta + 3 \tau \right)}}{81 \omega_{0}^{4} - 18 \omega_{0}^{2} + 1}$




```python
omega = omega_i[0] + epsilon * omega_1_new + epsilon**2 *  omega_2_new
omega
```




$\displaystyle - \frac{a \epsilon^{2} \omega_{0}^{2}}{36 \omega_{0}^{2} - 4} + \omega_{0}$




```python
t = sp.symbols('t')
x_new = epsilon * a * sp.cos(omega * t + beta)
x_new
```




$\displaystyle a \epsilon \cos{\left(\beta + t \left(- \frac{a \epsilon^{2} \omega_{0}^{2}}{36 \omega_{0}^{2} - 4} + \omega_{0}\right) \right)}$



We can find $x_3(t)$ after this, but we do not have to. This is because we will "reconstitute" the solution to $O(\epsilon^3)$

$$
\begin{align*}
    \begin{cases}
        x &= x_0 + \epsilon x_1 + \epsilon^2 x_2 + O(\epsilon^3)\\
        \omega &= \omega_0 + \epsilon \omega_1 + \epsilon^2 \omega_2 + O(\epsilon^3)\\
        \tau &= \omega t 
    \end{cases}
\end{align*}
$$

Establish the expansion for $x$ (found during the Lindstedt-Poincare analysis)


```python
x0_new = a * sp.cos(omega * t + beta) 
x1_new = x1_p.subs(tau, omega * t )
x2_new = x2_p.subs(tau, omega * t)
x_final = x0_new + epsilon * x1_new + epsilon ** 2 * x2_new
```


```python
x_final
```




$\displaystyle a \cos{\left(\beta + t \left(- \frac{a \epsilon^{2} \omega_{0}^{2}}{36 \omega_{0}^{2} - 4} + \omega_{0}\right) \right)} + \epsilon^{2} \left(- \frac{2 \omega_{0} \cos{\left(5 \beta + 5 t \left(- \frac{a \epsilon^{2} \omega_{0}^{2}}{36 \omega_{0}^{2} - 4} + \omega_{0}\right) \right)}}{225 \omega_{0}^{4} - 34 \omega_{0}^{2} + 1} + \frac{2 \omega_{0} \cos{\left(3 \beta + 3 t \left(- \frac{a \epsilon^{2} \omega_{0}^{2}}{36 \omega_{0}^{2} - 4} + \omega_{0}\right) \right)}}{81 \omega_{0}^{4} - 18 \omega_{0}^{2} + 1}\right) + \frac{2 \epsilon \sin{\left(3 \beta + 3 t \left(- \frac{a \epsilon^{2} \omega_{0}^{2}}{36 \omega_{0}^{2} - 4} + \omega_{0}\right) \right)}}{3 \left(9 \omega_{0}^{2} - 1\right)}$




```python
# From the initial conditions we know that omega0 = 1, a = 2, and beta = 0
x_final = x_final.subs([
    (omega_i[0], 1), (a, 2), (beta, 0)])
x_final
```




$\displaystyle \epsilon^{2} \left(\frac{\cos{\left(3 t \left(1 - \frac{\epsilon^{2}}{16}\right) \right)}}{32} - \frac{\cos{\left(5 t \left(1 - \frac{\epsilon^{2}}{16}\right) \right)}}{96}\right) + \frac{\epsilon \sin{\left(3 t \left(1 - \frac{\epsilon^{2}}{16}\right) \right)}}{12} + 2 \cos{\left(t \left(1 - \frac{\epsilon^{2}}{16}\right) \right)}$




```python
omega.subs([(a,2), (omega_i[0],1)])
```




$\displaystyle 1 - \frac{\epsilon^{2}}{16}$




```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, DOP853
from typing import List
```


```python
# Generate phase plane of Rayleigh's equation to confirm limit cycle 

def rayleigh(t, x, e):
    return [x[1], -x[0] - e*(x[1]**3 / 3 - x[1])]

def solve_diffeq(func, t, tspan, ic, parameters={}, algorithm='DOP853', stepsize=np.inf):
    return solve_ivp(fun=func, t_span=tspan, t_eval=t, y0=ic, method=algorithm, 
                     args=tuple(parameters.values()), atol=1e-8, rtol=1e-5, max_step=stepsize)

def phasePlane(x1, x2, func, params):
    X1, X2 = np.meshgrid(x1, x2)  # create grid
    u, v = np.zeros(X1.shape), np.zeros(X2.shape)
    NI, NJ = X1.shape
    for i in range(NI):
        for j in range(NJ):
            x = X1[i, j]
            y = X2[i, j]
            dx = func(0, (x, y), *params.values())  # compute values on grid
            u[i, j] = dx[0]
            v[i, j] = dx[1]
    M = np.hypot(u, v)
    u /= M
    v /= M
    return X1, X2, u, v, M

def DEplot(sys: object, tspan: tuple, x0: List[List[float]], 
           x: np.ndarray, y: np.ndarray, params: dict):
    if len(tspan) != 3:
        raise Exception('tspan should be tuple of size 3: (min, max, number of points).')
    # Set up the figure the way we want it to look
    plt.figure(figsize=(12, 9))

    X1, X2, dx1, dx2, M = phasePlane(
        x, y, sys, params
    )

    # Quiver plot
    plt.quiver(X1, X2, dx1, dx2, M, scale=None, pivot='mid')
    plt.grid()
    
    if tspan[0] < 0:
        t1 = np.linspace(0, tspan[0], tspan[2])
        t2 = np.linspace(0, tspan[1], tspan[2])
        if min(tspan) < 0:
            t_span1 = (np.max(t1), np.min(t1))
        else:
            t_span1 = (np.min(t1), np.max(t1))
        t_span2 = (np.min(t2), np.max(t2))
        for x0i in x0:
            sol1 = solve_diffeq(sys, t1, t_span1, x0i, params)
            plt.plot(sol1.y[0, :], sol1.y[1, :], '-r')
            sol2 = solve_diffeq(sys, t2, t_span2, x0i, params)
            plt.plot(sol2.y[0, :], sol2.y[1, :], '-r')
    else:
        t = np.linspace(tspan[0], tspan[1], tspan[2])
        t_span = (np.min(t), np.max(t))
        for x0i in x0:
            sol = solve_diffeq(sys, t, t_span, x0i, params)
            plt.plot(sol.y[0, :], sol.y[1, :], '-r')
        

    plt.xlim([np.min(x), np.max(x)])
    plt.ylim([np.min(y), np.max(y)])
    plt.show()
```


```python
x10 = np.arange(0, 10, 1)
x20 = np.arange(0, 10, 1)
x0 = np.stack((x10, x20), axis=-1)

p = {'e': 0.01}

x1 = np.linspace(-5, 5, 20)
x2 = np.linspace(-5, 5, 20)

DEplot(rayleigh, (-8, 8, 1000), x0, x1, x2, p)
```


```python
# Compare the approximation to the actual solution 
# let a = 2

tmax = 2
tmin = 30
tspan = np.linspace(tmin, tmax, 1000)
# ODE solver solution 
sol =  solve_diffeq(rayleigh, tspan, (tmin, tmax), [2, 0], p)

# Approximation 
def rayleigh_LP2(t, e):
    omega = 1 - e**2/16
    A = 2 * np.cos(t * omega)
    B = e * np.sin(3 * t * omega)/12
    C = e**2 * np.cos(3 * t * omega)/32
    D = -e**2 * np.cos(5 * t * omega) / 96
    return A + B + C + D

approx = rayleigh_LP2(tspan, 0.01)

plt.figure(figsize=(12, 9))
plt.plot(tspan, sol.y[0, :])
plt.plot(tspan, approx)
plt.grid(True)
plt.xlabel('$t$')
plt.ylabel('$x$')
plt.show()
```

### Better than first order Multiple Scale method

1ca242857ddccdedcdec111784c1aafe7fd8c53d
ipynb
Jupyter Notebook
notebooks/lindstedt_poincare_example.ipynb
smallpondtom/pNLsys
ebaea9c2945391fccd076f8006baee7066951097
MIT
2021-12-17T16:44:14.000Z
2022-02-22T05:40:43.000Z
notebooks/lindstedt_poincare_example.ipynb
smallpondtom/pNLsys
ebaea9c2945391fccd076f8006baee7066951097
MIT
notebooks/lindstedt_poincare_example.ipynb
smallpondtom/pNLsys
ebaea9c2945391fccd076f8006baee7066951097
MIT
Qwen/Qwen-72B
1. YES
2. YES