# Multiple Scales Stability for forced Duffing System 

Consider 

$$
\begin{gather*}
    \ddot{x} + \omega_0^2 x = -2\epsilon\mu\dot{x} - \epsilon\alpha x^3 + K\cos\Omega t
\end{gather*}
$$

To get the linear terms, introduce epsilon for later Taylor expansion 


```python
import sympy as sp
from math import factorial
```


```python
a0, a1 = sp.symbols('a_0 a_1', real=True)
epsilon = sp.Symbol('epsilon', real=True, positive=True)
a = a0 + epsilon * a1
a
```




$\displaystyle a_{0} + a_{1} \epsilon$




```python
Gamma0, Gamma1 = sp.symbols('Gamma_0 Gamma_1', real=True)
Gamma = Gamma0 + epsilon * Gamma1
Gamma
```




$\displaystyle \Gamma_{0} + \Gamma_{1} \epsilon$




```python
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)

taylor(a**3, epsilon, 4)
```




$\displaystyle a_{1}^{3} \epsilon^{3} + 3 a_{0} a_{1}^{2} \epsilon^{2} + 3 a_{0}^{2} a_{1} \epsilon + a_{0}^{3} + O\left(\epsilon^{4}\right)$




```python
taylor(1/a*sp.cos(Gamma), epsilon, 6)
```




$\displaystyle \frac{\cos{\left(\Gamma_{0} \right)}}{a_{0}} + \epsilon \left(- \frac{\Gamma_{1} \sin{\left(\Gamma_{0} \right)}}{a_{0}} - \frac{a_{1} \cos{\left(\Gamma_{0} \right)}}{a_{0}^{2}}\right) + \epsilon^{2} \left(- \frac{\Gamma_{1}^{2} \cos{\left(\Gamma_{0} \right)}}{2 a_{0}} + \frac{\Gamma_{1} a_{1} \sin{\left(\Gamma_{0} \right)}}{a_{0}^{2}} + \frac{a_{1}^{2} \cos{\left(\Gamma_{0} \right)}}{a_{0}^{3}}\right) + \epsilon^{3} \left(\frac{\Gamma_{1}^{3} \sin{\left(\Gamma_{0} \right)}}{6 a_{0}} + \frac{\Gamma_{1}^{2} a_{1} \cos{\left(\Gamma_{0} \right)}}{2 a_{0}^{2}} - \frac{\Gamma_{1} a_{1}^{2} \sin{\left(\Gamma_{0} \right)}}{a_{0}^{3}} - \frac{a_{1}^{3} \cos{\left(\Gamma_{0} \right)}}{a_{0}^{4}}\right) + \epsilon^{4} \left(\frac{\Gamma_{1}^{4} \cos{\left(\Gamma_{0} \right)}}{24 a_{0}} - \frac{\Gamma_{1}^{3} a_{1} \sin{\left(\Gamma_{0} \right)}}{6 a_{0}^{2}} - \frac{\Gamma_{1}^{2} a_{1}^{2} \cos{\left(\Gamma_{0} \right)}}{2 a_{0}^{3}} + \frac{\Gamma_{1} a_{1}^{3} \sin{\left(\Gamma_{0} \right)}}{a_{0}^{4}} + \frac{a_{1}^{4} \cos{\left(\Gamma_{0} \right)}}{a_{0}^{5}}\right) + \epsilon^{5} \left(- \frac{\Gamma_{1}^{5} \sin{\left(\Gamma_{0} \right)}}{120 a_{0}} - \frac{\Gamma_{1}^{4} a_{1} \cos{\left(\Gamma_{0} \right)}}{24 a_{0}^{2}} + \frac{\Gamma_{1}^{3} a_{1}^{2} \sin{\left(\Gamma_{0} \right)}}{6 a_{0}^{3}} + \frac{\Gamma_{1}^{2} a_{1}^{3} \cos{\left(\Gamma_{0} \right)}}{2 a_{0}^{4}} - \frac{\Gamma_{1} a_{1}^{4} \sin{\left(\Gamma_{0} \right)}}{a_{0}^{5}} - \frac{a_{1}^{5} \cos{\left(\Gamma_{0} \right)}}{a_{0}^{6}}\right) + O\left(\epsilon^{6}\right)$



stable when q is greater than zero


```python
mu, sigma, alpha, omega0 = sp.symbols('mu sigma alpha omega_0', real=True)
q = mu**2 + (sigma - 3 * alpha * a0**2 / 8 / omega0) * (sigma - 9 * alpha * a0**2 / 8 / omega0)
q
```




$\displaystyle \mu^{2} + \left(- \frac{9 a_{0}^{2} \alpha}{8 \omega_{0}} + \sigma\right) \left(- \frac{3 a_{0}^{2} \alpha}{8 \omega_{0}} + \sigma\right)$



Substitute in $\sigma$ as a function of $a_0$ into $q$ and find when $q= 0$ (see inverted frequency response diagram)

This is the top branch


```python
k = sp.Symbol('k')
topBranchCriterium = sp.simplify(q.subs([
    (sigma, 3/8*alpha/omega0*a0**2 + sp.sqrt(k**2 / 4 / omega0**2 / a0**2 - mu**2))
]))
topBranchCriterium = sp.Eq(topBranchCriterium, 0)
topBranchCriterium
```




$\displaystyle \frac{0.25 \left(- 1.5 a_{0}^{2} \alpha \sqrt{- 4 \mu^{2} + \frac{k^{2}}{a_{0}^{2} \omega_{0}^{2}}} + \frac{1.0 k^{2}}{a_{0}^{2} \omega_{0}}\right)}{\omega_{0}} = 0$



This is the bottom branch


```python
bottomBranchCriterium = sp.simplify(q.subs([
    (sigma, 3/8*alpha/omega0*a0**2 - sp.sqrt(k**2 / 4 / omega0**2 / a0**2 - mu**2))
]))
bottomBranchCriterium = sp.Eq(bottomBranchCriterium, 0)
bottomBranchCriterium
```




$\displaystyle \frac{0.25 \left(1.5 a_{0}^{2} \alpha \sqrt{- 4 \mu^{2} + \frac{k^{2}}{a_{0}^{2} \omega_{0}^{2}}} + \frac{1.0 k^{2}}{a_{0}^{2} \omega_{0}}\right)}{\omega_{0}} = 0$




```python
temp = sp.Eq(
    (1.5*a0**2*alpha * sp.sqrt(-4*mu + k**2 / a0**2 / omega0**2))**2,
    (-k**2 / a0**2 / omega0)**2
)
sp.solve(temp, a0)[0]
```




$\displaystyle - \begin{cases} 0.467327632592034 \sqrt{\frac{0.286178560638333 k^{2}}{\mu \omega_{0}^{2}} - \left(\frac{0.081898168569028 k^{4}}{\mu^{2} \omega_{0}^{4}} - \left(- \frac{k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333}\right)^{0.5} - \left(- \frac{0.046875 k^{6}}{\mu^{3} \omega_{0}^{6} \left(\frac{0.081898168569028 k^{4}}{\mu^{2} \omega_{0}^{4}} - \left(- \frac{k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333}\right)^{0.5}} + \frac{0.163796337138056 k^{4}}{\mu^{2} \omega_{0}^{4}} + \left(- \frac{k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333}\right)^{0.5}} & \text{for}\: \frac{0.111111111111111 k^{4}}{\alpha^{2} \mu \omega_{0}^{2}} = 0.0 \\0.556957466689064 \sqrt{\frac{0.201481862216914 k^{2}}{\mu \omega_{0}^{2}} - 1.0 \left(\frac{0.0405949408023956 k^{4}}{\mu^{2} \omega_{0}^{4}} + 1.0 \left(\left(\frac{0.0037078857421875 k^{16}}{\alpha^{4} \mu^{6} \omega_{0}^{12}} - \frac{k^{12}}{\alpha^{6} \mu^{3} \omega_{0}^{6}}\right)^{0.5} + \frac{0.0608924112035933 k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333} + \frac{k^{4}}{\alpha^{2} \mu \omega_{0}^{2} \left(\left(\frac{0.0037078857421875 k^{16}}{\alpha^{4} \mu^{6} \omega_{0}^{12}} - \frac{k^{12}}{\alpha^{6} \mu^{3} \omega_{0}^{6}}\right)^{0.5} + \frac{0.0608924112035933 k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333}}\right)^{0.5} - \left(- \frac{0.0163582885389041 k^{6}}{\mu^{3} \omega_{0}^{6} \left(\frac{0.0405949408023956 k^{4}}{\mu^{2} \omega_{0}^{4}} + 1.0 \left(\left(\frac{0.0037078857421875 k^{16}}{\alpha^{4} \mu^{6} \omega_{0}^{12}} - \frac{k^{12}}{\alpha^{6} \mu^{3} \omega_{0}^{6}}\right)^{0.5} + \frac{0.0608924112035933 k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333} + \frac{k^{4}}{\alpha^{2} \mu \omega_{0}^{2} \left(\left(\frac{0.0037078857421875 k^{16}}{\alpha^{4} \mu^{6} \omega_{0}^{12}} - \frac{k^{12}}{\alpha^{6} \mu^{3} \omega_{0}^{6}}\right)^{0.5} + \frac{0.0608924112035933 k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333}}\right)^{0.5}} + \frac{0.0811898816047911 k^{4}}{\mu^{2} \omega_{0}^{4}} - 1.0 \left(\left(\frac{0.0037078857421875 k^{16}}{\alpha^{4} \mu^{6} \omega_{0}^{12}} - \frac{k^{12}}{\alpha^{6} \mu^{3} \omega_{0}^{6}}\right)^{0.5} + \frac{0.0608924112035933 k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333} - \frac{k^{4}}{\alpha^{2} \mu \omega_{0}^{2} \left(\left(\frac{0.0037078857421875 k^{16}}{\alpha^{4} \mu^{6} \omega_{0}^{12}} - \frac{k^{12}}{\alpha^{6} \mu^{3} \omega_{0}^{6}}\right)^{0.5} + \frac{0.0608924112035933 k^{8}}{\alpha^{2} \mu^{3} \omega_{0}^{6}}\right)^{0.333333333333333}}\right)^{0.5}} & \text{otherwise} \end{cases}$



Equation for sigma - using plus sign


```python
sigma_func = 3*alpha/omega0*a0**2/8 + sp.sqrt(k**2 / 4 / omega0**2 / a0**2 - mu**2)
sigma_func
```




$\displaystyle \frac{3 a_{0}^{2} \alpha}{8 \omega_{0}} + \sqrt{- \mu^{2} + \frac{k^{2}}{4 a_{0}^{2} \omega_{0}^{2}}}$



Search for $a_0$ values such that the slope is zero - notice this matches the solutions above 


```python
sp.solve(sp.diff(sigma_func, a0), a0)

# not solvable with sympy
# but let's assume that it is the same as the huge answer shown above
```

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