# Conservative SDOF - Multiple Scales 

- Introduces multiple time scales (Homogenation)
- Treate damped systems easier then L-P
- Built-in stability

Introduce new independent time variables

$$
\begin{gather*}
    T_n = \epsilon^n t
\end{gather*}
$$

and 

$$
\begin{align*}
    \frac{d}{dt} &= \frac{\partial}{\partial T_0}\frac{dT_0}{dt} + \frac{\partial}{\partial T_1}\frac{dT_1}{dt} +  \frac{\partial}{\partial T_2}\frac{dT_2}{dt}  + \cdots\\
    &= \frac{\partial}{\partial T_0} + \epsilon \frac{\partial}{\partial T_1} + \epsilon^2 \frac{\partial}{\partial T_2} + \cdots\\
    &= D_0 + \epsilon D_1 + \epsilon^2 D_2 + \cdots
\end{align*}
$$

$$
\begin{align*}
    \frac{d^2}{dt^2} &= \left( D_0 + \epsilon D_1 + \epsilon^2 D_2 + \cdots \right)^2
\end{align*}
$$

Introducing the Expansion for $x(t)$

$$
\begin{align*}
    x(t) &= x_0(T_0,T_1,T_2,\cdots) + \epsilon x_1(T_0,T_1,T_2,\cdots) + \epsilon^2(T_0,T_1,T_2,\cdots) + \epsilon^3(T_0,T_1,T_2,\cdots) + O(\epsilon^4) 
\end{align*}
$$


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


```python
# Functions for multiple scales 

# Function for Time operator 
def Dt(f, n, Ts, e=sp.Symbol('epsilon')):
    if n==1:
        return sp.expand(sum([e**i * sp.diff(f, T_i) for i, T_i in enumerate(Ts)]))
    return Dt(Dt(f, 1, Ts, e), n-1, Ts, e)

def collect_epsilon(f, e=sp.Symbol('epsilon')):
    N = sp.degree(f, e)
    f_temp = f
    collected_dict = {}
    for i in range(N, 0, -1):
        collected_term = f_temp.coeff(e**i)
        collected_dict[e**i] = collected_term
        delete_terms = sp.expand(e**i * collected_term)
        f_temp -= delete_terms
    collected_dict[e**0] = f_temp
    return collected_dict
```


```python
N = 3

f = sp.Function('f')
t = sp.Symbol('t', real=True)

# Define the symbolic parameters
epsilon = sp.symbols('epsilon')

T_i = sp.symbols('T_(0:' + str(N) + ')', real=True)
alpha_i = sp.symbols('alpha_(2:' + str(N+1) + ')', real=True)
omega0 = sp.Symbol('omega_0', real=True)

# x0 = sp.Function('x_0')(*T_i)
x1 = sp.Function('x_1')(*T_i)
x2 = sp.Function('x_2')(*T_i)
x3 = sp.Function('x_3')(*T_i)

# Expansion for x(t)
x_e = epsilon*x1 + epsilon**2 * x2 + epsilon**3 * x3
x_e
```




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




```python
# Derivatives with time operators
xd = Dt(x_e, 1, T_i, epsilon)
xdd = Dt(x_e, 2, T_i, epsilon)

# EOM
EOM = xdd + sp.expand(omega0**2 * x_e) + sp.expand(sum([alpha_i[i-2] * x_e**i for i in range(2,N+1)]))
EOM
```




$\displaystyle \alpha_{2} \epsilon^{6} \operatorname{x_{3}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + 2 \alpha_{2} \epsilon^{5} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \alpha_{2} \epsilon^{4} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \alpha_{2} \epsilon^{4} \operatorname{x_{2}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + 2 \alpha_{2} \epsilon^{3} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \alpha_{2} \epsilon^{2} \operatorname{x_{1}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + \alpha_{3} \epsilon^{9} \operatorname{x_{3}}^{3}{\left(T_{0},T_{1},T_{2} \right)} + 3 \alpha_{3} \epsilon^{8} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + 3 \alpha_{3} \epsilon^{7} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + 3 \alpha_{3} \epsilon^{7} \operatorname{x_{2}}^{2}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 6 \alpha_{3} \epsilon^{6} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \alpha_{3} \epsilon^{6} \operatorname{x_{2}}^{3}{\left(T_{0},T_{1},T_{2} \right)} + 3 \alpha_{3} \epsilon^{5} \operatorname{x_{1}}^{2}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 3 \alpha_{3} \epsilon^{5} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{2}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + 3 \alpha_{3} \epsilon^{4} \operatorname{x_{1}}^{2}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \alpha_{3} \epsilon^{3} \operatorname{x_{1}}^{3}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{7} \frac{\partial^{2}}{\partial T_{2}^{2}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{6} \frac{\partial^{2}}{\partial T_{2}^{2}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{6} \frac{\partial^{2}}{\partial T_{2}\partial T_{1}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{5} \frac{\partial^{2}}{\partial T_{2}^{2}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{5} \frac{\partial^{2}}{\partial T_{1}^{2}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{5} \frac{\partial^{2}}{\partial T_{2}\partial T_{1}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{5} \frac{\partial^{2}}{\partial T_{2}\partial T_{0}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{4} \frac{\partial^{2}}{\partial T_{1}^{2}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{4} \frac{\partial^{2}}{\partial T_{2}\partial T_{1}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{4} \frac{\partial^{2}}{\partial T_{2}\partial T_{0}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{4} \frac{\partial^{2}}{\partial T_{1}\partial T_{0}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{3} \omega_{0}^{2} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{3} \frac{\partial^{2}}{\partial T_{1}^{2}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{3} \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{3} \frac{\partial^{2}}{\partial T_{2}\partial T_{0}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{3} \frac{\partial^{2}}{\partial T_{1}\partial T_{0}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{2} \omega_{0}^{2} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon^{2} \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \epsilon^{2} \frac{\partial^{2}}{\partial T_{1}\partial T_{0}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon \omega_{0}^{2} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + \epsilon \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)}$




```python
# Ordered Equations by epsilon 
epsilon_Eq = collect_epsilon(EOM)
epsilon0_Eq = sp.Eq(epsilon_Eq[epsilon**0], 0)
epsilon0_Eq
```




$\displaystyle \text{True}$




```python
epsilon1_Eq = sp.Eq(epsilon_Eq[epsilon**1], 0)
epsilon1_Eq
```




$\displaystyle \omega_{0}^{2} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} = 0$




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




$\displaystyle \alpha_{2} \operatorname{x_{1}}^{2}{\left(T_{0},T_{1},T_{2} \right)} + \omega_{0}^{2} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \frac{\partial^{2}}{\partial T_{1}\partial T_{0}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} = 0$




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




$\displaystyle 2 \alpha_{2} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \alpha_{3} \operatorname{x_{1}}^{3}{\left(T_{0},T_{1},T_{2} \right)} + \omega_{0}^{2} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + \frac{\partial^{2}}{\partial T_{1}^{2}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \frac{\partial^{2}}{\partial T_{2}\partial T_{0}} \operatorname{x_{1}}{\left(T_{0},T_{1},T_{2} \right)} + 2 \frac{\partial^{2}}{\partial T_{1}\partial T_{0}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} = 0$




```python
# Find the solution for epsilon-1
A = sp.Function('A')(*T_i[1::])
x1_sol = A * sp.exp(sp.I * omega0 * T_i[0]) + sp.conjugate(A) * sp.exp(-sp.I * omega0 * T_i[0])
x1_sol
```




$\displaystyle A{\left(T_{1},T_{2} \right)} e^{i T_{0} \omega_{0}} + e^{- i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}$




```python
# Update the epsilon-2 equation
epsilon2_Eq = epsilon2_Eq.subs(x1, x1_sol).doit()
epsilon2_Eq = sp.expand(epsilon2_Eq)
epsilon2_Eq
```




$\displaystyle \alpha_{2} A^{2}{\left(T_{1},T_{2} \right)} e^{2 i T_{0} \omega_{0}} + 2 \alpha_{2} A{\left(T_{1},T_{2} \right)} \overline{A{\left(T_{1},T_{2} \right)}} + \alpha_{2} e^{- 2 i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}^{2} + \omega_{0}^{2} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + 2 i \omega_{0} e^{i T_{0} \omega_{0}} \frac{\partial}{\partial T_{1}} A{\left(T_{1},T_{2} \right)} - 2 i \omega_{0} e^{- i T_{0} \omega_{0}} \overline{\frac{\partial}{\partial T_{1}} A{\left(T_{1},T_{2} \right)}} + \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} = 0$



The secular terms will be cancelled out by 

$$
\begin{gather*}
    D_1 A = 0
\end{gather*}
$$


```python
epsilon2_Eq = epsilon2_Eq.subs(sp.diff(A, T_i[1]), 0)
epsilon2_Eq
```




$\displaystyle \alpha_{2} A^{2}{\left(T_{1},T_{2} \right)} e^{2 i T_{0} \omega_{0}} + 2 \alpha_{2} A{\left(T_{1},T_{2} \right)} \overline{A{\left(T_{1},T_{2} \right)}} + \alpha_{2} e^{- 2 i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}^{2} + \omega_{0}^{2} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} + \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{2}}{\left(T_{0},T_{1},T_{2} \right)} = 0$



The particular solution of $x_2$ is 

$$
\begin{gather*}
    x_2 = \frac{\alpha_2 A^2}{3 \omega_0^2} e^{2i\omega_0 T_0} - \frac{\alpha_2 }{\omega^2_0}A\overline{A} + cc
\end{gather*}
$$


```python
x2_p = alpha_i[0] * A**2 / 3/omega0**2 * sp.exp(2*sp.I*omega0*T_i[0]) - alpha_i[0]/omega0**2 * A * sp.conjugate(A) 
x2_p
```




$\displaystyle \frac{\alpha_{2} A^{2}{\left(T_{1},T_{2} \right)} e^{2 i T_{0} \omega_{0}}}{3 \omega_{0}^{2}} - \frac{\alpha_{2} A{\left(T_{1},T_{2} \right)} \overline{A{\left(T_{1},T_{2} \right)}}}{\omega_{0}^{2}}$




```python
epsilon3_Eq = epsilon3_Eq.subs([
    (sp.diff(A, T_i[1]), 0), (x1, x1_sol), (x2, x2_p)
]).doit()
epsilon3_Eq = sp.expand(epsilon3_Eq)
epsilon3_Eq = epsilon3_Eq.subs(sp.diff(A, T_i[1], 2), 0)
epsilon3_Eq
```




$\displaystyle \frac{2 \alpha_{2}^{2} A^{3}{\left(T_{1},T_{2} \right)} e^{3 i T_{0} \omega_{0}}}{3 \omega_{0}^{2}} - \frac{4 \alpha_{2}^{2} A^{2}{\left(T_{1},T_{2} \right)} e^{i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}}{3 \omega_{0}^{2}} - \frac{2 \alpha_{2}^{2} A{\left(T_{1},T_{2} \right)} e^{- i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}^{2}}{\omega_{0}^{2}} + \frac{8 i \alpha_{2} A{\left(T_{1},T_{2} \right)} e^{2 i T_{0} \omega_{0}} \frac{\partial}{\partial T_{1}} A{\left(T_{1},T_{2} \right)}}{3 \omega_{0}} + \alpha_{3} A^{3}{\left(T_{1},T_{2} \right)} e^{3 i T_{0} \omega_{0}} + 3 \alpha_{3} A^{2}{\left(T_{1},T_{2} \right)} e^{i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}} + 3 \alpha_{3} A{\left(T_{1},T_{2} \right)} e^{- i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}^{2} + \alpha_{3} e^{- 3 i T_{0} \omega_{0}} \overline{A{\left(T_{1},T_{2} \right)}}^{3} + \omega_{0}^{2} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} + 2 i \omega_{0} e^{i T_{0} \omega_{0}} \frac{\partial}{\partial T_{2}} A{\left(T_{1},T_{2} \right)} - 2 i \omega_{0} e^{- i T_{0} \omega_{0}} \overline{\frac{\partial}{\partial T_{2}} A{\left(T_{1},T_{2} \right)}} + \frac{\partial^{2}}{\partial T_{0}^{2}} \operatorname{x_{3}}{\left(T_{0},T_{1},T_{2} \right)} = 0$



The to cancel out the secular term we let 

$$
\begin{gather*}
    2i\omega_0 D_2 A + \dfrac{9\alpha_3 \omega_0^2 - 10\alpha_2^2 }{3\omega_0^2}A^2\overline{A} = 0
\end{gather*}
$$

Question: What if the secular terms arising from $i\omega_0$ and $-i \omega_0$ are handled together - do we get a single real equation?

Substitute the polar $A$

$$
\begin{gather*}
    A = \dfrac{1}{2}a e^{i\beta}
\end{gather*}
$$



```python
x3_sec = sp.Eq(2*sp.I*omega0*sp.diff(A, T_i[2]) + (9*alpha_i[1]*omega0**2 - 10*alpha_i[0]**2)/3/omega0**2 * A**2 * sp.conjugate(A), 0)

a = sp.Symbol('a', real=True)
beta = sp.Symbol('beta', real=True)

temp = x3_sec.subs(A, a*sp.exp(sp.I * beta)/2)
temp
```




$\displaystyle \frac{a^{3} \left(- \frac{10 \alpha_{2}^{2}}{3} + 3 \alpha_{3} \omega_{0}^{2}\right) e^{i \beta}}{8 \omega_{0}^{2}} + 2 i \omega_{0} \frac{\partial}{\partial T_{2}} \frac{a e^{i \beta}}{2} = 0$




```python
temp = sp.expand(temp)
temp_im = sp.im(temp.lhs)
temp_im
```




$\displaystyle - \frac{5 a^{3} \alpha_{2}^{2} \sin{\left(\beta \right)}}{12 \omega_{0}^{2}} + \frac{3 a^{3} \alpha_{3} \sin{\left(\beta \right)}}{8} + 2 \omega_{0} \operatorname{re}{\left(\frac{\partial}{\partial T_{2}} \frac{a e^{i \beta}}{2}\right)}$




```python
temp_re = sp.re(temp.lhs)
temp_re
```




$\displaystyle - \frac{5 a^{3} \alpha_{2}^{2} \cos{\left(\beta \right)}}{12 \omega_{0}^{2}} + \frac{3 a^{3} \alpha_{3} \cos{\left(\beta \right)}}{8} - 2 \omega_{0} \operatorname{im}{\left(\frac{\partial}{\partial T_{2}} \frac{a e^{i \beta}}{2}\right)}$



Thus separating into real and imaginary parts we obtain

$$
\begin{align*}
    \omega_0 D_2 a &=0\\
    omega_0 a D_2 \beta + \dfrac{10\alpha_2^2 - 9\alpha_3\omega_0^2}{24\omega_0^2}a^3 &= 0
\end{align*}
$$

$a$ is a constant and 

$$
\begin{gather*}
    D_2\beta = - \dfrac{10\alpha_2^2 - 9\alpha_3\omega_0^3a}{24\omega_0^2}a^3
    \beta = \dfrac{9\alpha_3\omega_0^2 - 10\alpha_2^2  }{24\omega_0^3}a^2 T_2 + \beta_0
\end{gather*}
$$

Here $\beta_0$ is a constant. Now using $T_2 = \epsilon^2 t$ we find that 

$$
\begin{gather*}
    A = \dfrac{1}{2}a \exp\left[ i\dfrac{9\alpha_3\omega_0^2 - 10\alpha_2^2  }{24\omega_0^3}a^3 \epsilon^2 t + i\beta_0 \right]
\end{gather*}
$$

and substituting in the expressions for $x_1$ and $x_2$ into the equations we have, we obtain the following final results

$$
\begin{gather*}
    x = \epsilon a \cos(\omega t + \beta_0) - \dfrac{\epsilon^2 a^2\alpha_2}{2\omega_0^2}\left[ 1 - \dfrac{1}{3}\cos(2\omega t + 2\beta_0) \right] + O(\epsilon^3)
\end{gather*}
$$

where 

$$
\begin{gather*}
    \omega = \omega_0 \left[ 1 + \dfrac{9\alpha_3 \omega_0^2 - 10\alpha_2^2}{24\omega_0^4}\epsilon^2 a^2 \right] + O(\epsilon^3)
\end{gather*}
$$

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