# Conservative SDOF - Lindstedt Poincare

The straight-forward expansion fails in part because the frequency is fixed in the expansion. Thus, secular terms arise which we cannot eliminate.

$$
\begin{gather*}
    \ddot{x} + \omega_0^2 x + \sum_{n=2}^{\infty} \alpha_n x^n.
\end{gather*}
$$

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


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


```python
N = 3

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


x = sp.Function('x')(tau)
xdot = sp.Derivative(x, tau)
xddot = sp.Derivative(xdot, tau)
            
x1 = sp.Function('x_1')(tau)
x2 = sp.Function('x_2')(tau)
x3 = sp.Function('x_3')(tau)

# EOM
EOM = omega**2 * xddot + omega_i[0]**2 * x + sum([alpha_i[i-2] * x**i for i in range(2,N+1)])
EOM
```




$\displaystyle \alpha_{2} x^{2}{\left(\tau \right)} + \alpha_{3} x^{3}{\left(\tau \right)} + \omega^{2} \frac{d^{2}}{d \tau^{2}} x{\left(\tau \right)} + \omega_{0}^{2} x{\left(\tau \right)}$




```python
x_i = (x1, x2, x3)
x_e = sum([epsilon**i * x_i[i-1] for i in range(1,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)}$




```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 \alpha_{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)}\right)^{2} + \alpha_{3} \left(\epsilon^{3} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{2} \operatorname{x_{2}}{\left(\tau \right)} + \epsilon \operatorname{x_{1}}{\left(\tau \right)}\right)^{3} + \omega_{0}^{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)}\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)}\right)$




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




$\displaystyle \alpha_{2} \epsilon^{6} \operatorname{x_{3}}^{2}{\left(\tau \right)} + 2 \alpha_{2} \epsilon^{5} \operatorname{x_{2}}{\left(\tau \right)} \operatorname{x_{3}}{\left(\tau \right)} + 2 \alpha_{2} \epsilon^{4} \operatorname{x_{1}}{\left(\tau \right)} \operatorname{x_{3}}{\left(\tau \right)} + \alpha_{2} \epsilon^{4} \operatorname{x_{2}}^{2}{\left(\tau \right)} + 2 \alpha_{2} \epsilon^{3} \operatorname{x_{1}}{\left(\tau \right)} \operatorname{x_{2}}{\left(\tau \right)} + \alpha_{2} \epsilon^{2} \operatorname{x_{1}}^{2}{\left(\tau \right)} + \alpha_{3} \epsilon^{9} \operatorname{x_{3}}^{3}{\left(\tau \right)} + 3 \alpha_{3} \epsilon^{8} \operatorname{x_{2}}{\left(\tau \right)} \operatorname{x_{3}}^{2}{\left(\tau \right)} + 3 \alpha_{3} \epsilon^{7} \operatorname{x_{1}}{\left(\tau \right)} \operatorname{x_{3}}^{2}{\left(\tau \right)} + 3 \alpha_{3} \epsilon^{7} \operatorname{x_{2}}^{2}{\left(\tau \right)} \operatorname{x_{3}}{\left(\tau \right)} + 6 \alpha_{3} \epsilon^{6} \operatorname{x_{1}}{\left(\tau \right)} \operatorname{x_{2}}{\left(\tau \right)} \operatorname{x_{3}}{\left(\tau \right)} + \alpha_{3} \epsilon^{6} \operatorname{x_{2}}^{3}{\left(\tau \right)} + 3 \alpha_{3} \epsilon^{5} \operatorname{x_{1}}^{2}{\left(\tau \right)} \operatorname{x_{3}}{\left(\tau \right)} + 3 \alpha_{3} \epsilon^{5} \operatorname{x_{1}}{\left(\tau \right)} \operatorname{x_{2}}^{2}{\left(\tau \right)} + 3 \alpha_{3} \epsilon^{4} \operatorname{x_{1}}^{2}{\left(\tau \right)} \operatorname{x_{2}}{\left(\tau \right)} + \alpha_{3} \epsilon^{3} \operatorname{x_{1}}^{3}{\left(\tau \right)} + \epsilon^{7} \omega_{2}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + 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)} + 2 \epsilon^{5} \omega_{0} \omega_{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + \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_{2}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + 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_{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^{3} \omega_{0}^{2} \operatorname{x_{3}}{\left(\tau \right)} + \epsilon^{3} \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + 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_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon^{2} \omega_{0}^{2} \operatorname{x_{2}}{\left(\tau \right)} + \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)} + \epsilon \omega_{0}^{2} \operatorname{x_{1}}{\left(\tau \right)} + \epsilon \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)}$




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




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




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




$\displaystyle \alpha_{2} \operatorname{x_{1}}^{2}{\left(\tau \right)} + \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} \omega_{1} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} = 0$




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




$\displaystyle 2 \alpha_{2} \operatorname{x_{1}}{\left(\tau \right)} \operatorname{x_{2}}{\left(\tau \right)} + \alpha_{3} \operatorname{x_{1}}^{3}{\left(\tau \right)} + \omega_{0}^{2} \operatorname{x_{3}}{\left(\tau \right)} + \omega_{0}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{3}}{\left(\tau \right)} + 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_{1}^{2} \frac{d^{2}}{d \tau^{2}} \operatorname{x_{1}}{\left(\tau \right)} = 0$




```python
sp.dsolve(epsilon_1_Eq, x1) 
```




$\displaystyle \operatorname{x_{1}}{\left(\tau \right)} = C_{1} \sin{\left(\tau \right)} + C_{2} \cos{\left(\tau \right)}$



Not as convenient as working with the polar form


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




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



Update $\epsilon^2$ equation 


```python
epsilon_2_Eq = sp.expand(TR7(epsilon_2_Eq.subs(x1, x1_polar)))
epsilon_2_Eq
```




$\displaystyle \frac{a^{2} \alpha_{2} \cos{\left(2 \beta + 2 \tau \right)}}{2} + \frac{a^{2} \alpha_{2}}{2} + \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} \omega_{1} \frac{\partial^{2}}{\partial \tau^{2}} a \cos{\left(\beta + \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


```python
epsilon_2_Eq = epsilon_2_Eq.subs(omega_i[1], 0)
epsilon_2_Eq
```




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




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




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




```python
x2_p = a**2 * alpha_i[0] * sp.cos(2*phi) / 6 / omega_i[0]**2 - a**2 * alpha_i[0] / 2 / omega_i[0]**2
x2_p
```




$\displaystyle \frac{a^{2} \alpha_{2} \cos{\left(2 \beta + 2 \tau \right)}}{6 \omega_{0}^{2}} - \frac{a^{2} \alpha_{2}}{2 \omega_{0}^{2}}$



We only take the particular solution into consideration, and update $\epsilon^3$ equation


```python
epsilon_3_Eq = epsilon_3_Eq.subs([
    (omega_i[1], 0), (x1, x1_polar), (x2, x2_p)
])
```


```python
epsilon_3_Eq = TR8(epsilon_3_Eq)
epsilon_3_Eq
```




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



We remove the secular term by setting

$$
\begin{gather*}
    \omega_2 = \dfrac{(9 \alpha_3 \omega^2_0 - 10 \alpha_2^2) a^2}{24\omega_0^3}
\end{gather*}
$$


```python
omega_2_new = (9 * alpha_i[1] * omega_i[0]**2 - 10 * alpha_i[0]**2) * a**2 / 24 / omega_i[0]**3
omega = omega_i[0] + epsilon**2 *  omega_2_new
omega
```




$\displaystyle \frac{a^{2} \epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right)}{24 \omega_{0}^{3}} + \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^{2} \epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right)}{24 \omega_{0}^{3}} + \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 &= \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*}
$$


```python
# Imposing provided initial conditions 
a0 = sp.symbols('a_0', real=True, positive=True)
beta0 = sp.symbols('beta_0', real=True, positive=True)

x0 = epsilon * a0 * sp.cos(beta0)
x0
```




$\displaystyle a_{0} \epsilon \cos{\left(\beta_{0} \right)}$




```python
xdot0 = -epsilon * a0 * sp.sin(beta0)
xdot0
```




$\displaystyle - a_{0} \epsilon \sin{\left(\beta_{0} \right)}$



Expand the $x_1$ homogeneous coefficients


```python
A_i = sp.symbols('A_(1:' + str(N+1) + ')')
B_i = sp.symbols('B_(1:' + str(N+1) + ')')

a_i = sum([epsilon ** i * A_i[i] for i in range(N)])
a_i
```




$\displaystyle A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}$




```python
beta_i = sum([epsilon ** i * B_i[i] for i in range(N)])
beta_i
```




$\displaystyle B_{1} + B_{2} \epsilon + B_{3} \epsilon^{2}$



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


```python
omega = omega.subs(a, a_i)
omega
```




$\displaystyle \frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{24 \omega_{0}^{3}} + \omega_{0}$



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


```python
x1_new = a_i * sp.cos(omega * t + beta_i) 
x2_new = x2_p.subs([   
    (tau, omega * t), (a, a_i), (beta, beta_i)
])
x_new = epsilon * x1_new + epsilon ** 2 * x2_new
```


```python
x_new
```




$\displaystyle \epsilon^{2} \left(\frac{\alpha_{2} \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2} \cos{\left(2 B_{1} + 2 B_{2} \epsilon + 2 B_{3} \epsilon^{2} + 2 t \left(\frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{24 \omega_{0}^{3}} + \omega_{0}\right) \right)}}{6 \omega_{0}^{2}} - \frac{\alpha_{2} \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{2 \omega_{0}^{2}}\right) + \epsilon \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right) \cos{\left(B_{1} + B_{2} \epsilon + B_{3} \epsilon^{2} + t \left(\frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{24 \omega_{0}^{3}} + \omega_{0}\right) \right)}$




```python
xdot_new = sp.diff(x_new, t)
xdot_new
```




$\displaystyle - \frac{\alpha_{2} \epsilon^{2} \left(\frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{12 \omega_{0}^{3}} + 2 \omega_{0}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2} \sin{\left(2 B_{1} + 2 B_{2} \epsilon + 2 B_{3} \epsilon^{2} + 2 t \left(\frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{24 \omega_{0}^{3}} + \omega_{0}\right) \right)}}{6 \omega_{0}^{2}} - \epsilon \left(\frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{24 \omega_{0}^{3}} + \omega_{0}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right) \sin{\left(B_{1} + B_{2} \epsilon + B_{3} \epsilon^{2} + t \left(\frac{\epsilon^{2} \left(- 10 \alpha_{2}^{2} + 9 \alpha_{3} \omega_{0}^{2}\right) \left(A_{1} + A_{2} \epsilon + A_{3} \epsilon^{2}\right)^{2}}{24 \omega_{0}^{3}} + \omega_{0}\right) \right)}$



For IC analysis, taylor $x$ and $\dot{x}$ up to order $\epsilon^3$


```python
# The built in sympy series command takes forever (bad algorithm) so we just implement our original Taylor series function
def taylor(expr, var, upto):
    res = expr.subs(var, 0)
    dexpr = sp.diff(expr, var)
    for n in range(1,upto+1):
        res += dexpr.subs(var, 0) / factorial(n) * var**n
        dexpr = sp.diff(dexpr, var)
    return res + sp.O(var**upto)

xTaylor = taylor(x_new, epsilon, 3)
xTaylor
```




$\displaystyle \epsilon^{2} \left(\frac{A_{1}^{2} \alpha_{2} \cos{\left(2 B_{1} + 2 \omega_{0} t \right)}}{6 \omega_{0}^{2}} - \frac{A_{1}^{2} \alpha_{2}}{2 \omega_{0}^{2}} - A_{1} B_{2} \sin{\left(B_{1} + \omega_{0} t \right)} + A_{2} \cos{\left(B_{1} + \omega_{0} t \right)}\right) + A_{1} \epsilon \cos{\left(B_{1} + \omega_{0} t \right)} + O\left(\epsilon^{3}\right)$




```python
xdotTaylor = taylor(xdot_new, epsilon, 3)
xdotTaylor
```




$\displaystyle \epsilon^{2} \left(- \frac{A_{1}^{2} \alpha_{2} \sin{\left(2 B_{1} + 2 \omega_{0} t \right)}}{3 \omega_{0}} - A_{1} B_{2} \omega_{0} \cos{\left(B_{1} + \omega_{0} t \right)} - A_{2} \omega_{0} \sin{\left(B_{1} + \omega_{0} t \right)}\right) - A_{1} \epsilon \omega_{0} \sin{\left(B_{1} + \omega_{0} t \right)} + O\left(\epsilon^{3}\right)$



Order $\epsilon$ initial condition for $x$



```python
temp = sp.collect(xTaylor.removeO(), epsilon, evaluate=False)
ICeq1 = sp.Eq(temp[epsilon].subs(t, 0), x0 / epsilon)
ICeq1
```




$\displaystyle A_{1} \cos{\left(B_{1} \right)} = a_{0} \cos{\left(\beta_{0} \right)}$



Order $\epsilon^2$ initial condition equation for $x$


```python
ICeq2 = sp.Eq(temp[epsilon**2].subs(t, 0), 0)
ICeq2
```




$\displaystyle \frac{A_{1}^{2} \alpha_{2} \cos{\left(2 B_{1} \right)}}{6 \omega_{0}^{2}} - \frac{A_{1}^{2} \alpha_{2}}{2 \omega_{0}^{2}} - A_{1} B_{2} \sin{\left(B_{1} \right)} + A_{2} \cos{\left(B_{1} \right)} = 0$



Order $\epsilon$ initial condition for $\dot{x}$


```python
temp = sp.collect(xdotTaylor.removeO(), epsilon, evaluate=False)
ICeq3 = sp.Eq(temp[epsilon].subs(t, 0), -a0 * omega_i[0] * sp.sin(beta0))
ICeq3
```




$\displaystyle - A_{1} \omega_{0} \sin{\left(B_{1} \right)} = - a_{0} \omega_{0} \sin{\left(\beta_{0} \right)}$



Order $\epsilon^2$ initial condition equation for $\dot{x}$


```python
ICeq4 = sp.Eq(temp[epsilon**2].subs(t, 0), 0)
ICeq4
```




$\displaystyle - \frac{A_{1}^{2} \alpha_{2} \sin{\left(2 B_{1} \right)}}{3 \omega_{0}} - A_{1} B_{2} \omega_{0} \cos{\left(B_{1} \right)} - A_{2} \omega_{0} \sin{\left(B_{1} \right)} = 0$




```python
res = sp.solve([ICeq1, ICeq2, ICeq3, ICeq4], [A_i[0], A_i[1], B_i[0], B_i[1]], force=True)
```


```python
res[0][0].subs(sp.atan(sp.tan(beta0/2)), beta0/2)
```




$\displaystyle a_{0}$




```python
sp.simplify(res[0][1].subs(sp.atan(sp.tan(beta0/2)), beta0/2))
```




$\displaystyle \frac{a_{0}^{2} \alpha_{2} \left(2 \sin{\left(2 \beta_{0} \right)} + \sin{\left(4 \beta_{0} \right)}\right)}{24 \omega_{0}^{2} \sin{\left(\beta_{0} \right)}}$




```python
res[0][2].subs(sp.atan(sp.tan(beta0/2)), beta0/2)
```




$\displaystyle \beta_{0}$




```python
sp.simplify(TR8(res[0][3].subs(sp.atan(sp.tan(beta0/2)), beta0/2)))
```




$\displaystyle - \frac{a_{0} \alpha_{2} \left(9 \sin{\left(\beta_{0} \right)} + \sin{\left(3 \beta_{0} \right)}\right)}{12 \omega_{0}^{2}}$



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