# Conservative SDOF - Harmonic Balance (2)

- Good for strongly and weakly NL systems
- Requires significant knowledge about system studied 
- Similar to Galerkin method, but in time

## Overview

Assume an harmonic expansion for $x(t)$

$$
\begin{gather*}
     x(t) = \sum^N_{m = 0} A_m \cos\left( m\omega t + m \beta_0\right)
\end{gather*}
$$

substitute this into the EOM and enforce zero projected residual, i.e. taking the inner product with orthogonal assumption.

$$
\begin{gather*}
    \left< EOM_{LHS},\phantom{-} \cos\left(m\omega t + m\beta_0\right) \right> = 0\\
    m = 0, 1, 2, \cdots, N
\end{gather*}
$$

which is 

$$
\begin{gather*}
    \int_0^{\frac{2\pi}{\omega} } \left( EOM_{LHS}\bullet \cos\left(m\omega t + m\beta_0\right) \right)dt
\end{gather*}
$$

Then we solve for $A_0, A_2, A_3, \cdots, A_n$ and $\omega$ in terms of $A_1$ and other system parameters. In this example we solve the following equation

$$
\begin{gather*}
    \ddot{x} + \omega_0^2 x + \alpha_2 x^2 + \alpha_3 x^3 = 0.
\end{gather*}
$$

## Assumed Harmonic Expansion - up to and including $\cos\left(2\left(\omega t + \beta \right)\right)$ and $\alpha_2 = 0$ no spurious solutions

Assume a harmonic expansion for $x(t)$ - will use this to satisfy the cubic equation, so include harmonic terms at twice the linear natrual frequency 

Note: phi = omega * t + beta0


```python
import sympy as sp
from sympy.simplify.fu import TR4, TR7, TR8, TR11
import numpy as np
from itertools import combinations_with_replacement
from math import factorial
```


```python
M = 2  # order of harmonic expansion
Ai = sp.symbols('A_(0:' + str(M+1) + ')', real=True)  # coefficient A
t, w, w0 = sp.symbols('t omega omega_0', real=True, positive=True)
b0 = sp.Symbol('beta_0')

# harmonic expansion of x(t)
phi = sp.Symbol('phi')
xh = sp.Function('x_h')(phi)
xh = sum([Ai[i] * sp.cos(i * phi) for i in range(M+1)])
xh
```




$\displaystyle A_{0} + A_{1} \cos{\left(\phi \right)} + A_{2} \cos{\left(2 \phi \right)}$




```python
# Defining the EOM
N = 3  # order of powered terms of alpha 
ai = sp.symbols('alpha_(2:' + str(N+1) + ')')
x = sp.Function('x')(phi)

# EOM
EOM = w**2 * sp.diff(x, phi, 2) + w0**2 * x + sum([ai[i-2] * x**i for i in range(2, N+1)])
EOM = sp.Eq(EOM, 0)
EOM = EOM.subs(ai[0], 0)  #! set alpha_2 = 0
EOM
```




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




```python
# Substitute x(t) for the harmonic expansion of x(t) into EOM
err = EOM.lhs.subs(x, xh).doit().expand()
err = TR7(err).expand()
err = TR8(err).expand()
err
```




$\displaystyle A_{0}^{3} \alpha_{3} + 3 A_{0}^{2} A_{1} \alpha_{3} \cos{\left(\phi \right)} + 3 A_{0}^{2} A_{2} \alpha_{3} \cos{\left(2 \phi \right)} + \frac{3 A_{0} A_{1}^{2} \alpha_{3} \cos{\left(2 \phi \right)}}{2} + \frac{3 A_{0} A_{1}^{2} \alpha_{3}}{2} + 3 A_{0} A_{1} A_{2} \alpha_{3} \cos{\left(\phi \right)} + 3 A_{0} A_{1} A_{2} \alpha_{3} \cos{\left(3 \phi \right)} + \frac{3 A_{0} A_{2}^{2} \alpha_{3} \cos{\left(4 \phi \right)}}{2} + \frac{3 A_{0} A_{2}^{2} \alpha_{3}}{2} + A_{0} \omega_{0}^{2} + \frac{3 A_{1}^{3} \alpha_{3} \cos{\left(\phi \right)}}{4} + \frac{A_{1}^{3} \alpha_{3} \cos{\left(3 \phi \right)}}{4} + \frac{3 A_{1}^{2} A_{2} \alpha_{3} \cos{\left(2 \phi \right)}}{2} + \frac{3 A_{1}^{2} A_{2} \alpha_{3} \cos{\left(4 \phi \right)}}{4} + \frac{3 A_{1}^{2} A_{2} \alpha_{3}}{4} + \frac{3 A_{1} A_{2}^{2} \alpha_{3} \cos{\left(\phi \right)}}{2} + \frac{3 A_{1} A_{2}^{2} \alpha_{3} \cos{\left(3 \phi \right)}}{4} + \frac{3 A_{1} A_{2}^{2} \alpha_{3} \cos{\left(5 \phi \right)}}{4} - A_{1} \omega^{2} \cos{\left(\phi \right)} + A_{1} \omega_{0}^{2} \cos{\left(\phi \right)} + \frac{3 A_{2}^{3} \alpha_{3} \cos{\left(2 \phi \right)}}{4} + \frac{A_{2}^{3} \alpha_{3} \cos{\left(6 \phi \right)}}{4} - 4 A_{2} \omega^{2} \cos{\left(2 \phi \right)} + A_{2} \omega_{0}^{2} \cos{\left(2 \phi \right)}$



Note there are terms with frequency content as high as $6\phi$

To find $A_0$, $\omega$, $A_2$, we project the error back on the assumed expansion and set the projections to zero

Note: same as setting the coefficients in front of each $\cos(n\phi)$ to zero for $n = 0, 1, 2$


```python
proj_eqn = {}
for i in range(M+1):
    proj_eqn['eq'+str(i)] = sp.integrate(sp.expand_trig(err * sp.cos(i * phi)), (phi, 0, 2*sp.pi))
```


```python
proj_eqn['eq0']
```




$\displaystyle 2 \pi A_{0}^{3} \alpha_{3} + 3 \pi A_{0} A_{1}^{2} \alpha_{3} + 3 \pi A_{0} A_{2}^{2} \alpha_{3} + 2 \pi A_{0} \omega_{0}^{2} + \frac{3 \pi A_{1}^{2} A_{2} \alpha_{3}}{2}$




```python
proj_eqn['eq1']
```




$\displaystyle 3 \pi A_{0}^{2} A_{1} \alpha_{3} + 3 \pi A_{0} A_{1} A_{2} \alpha_{3} + \frac{3 \pi A_{1}^{3} \alpha_{3}}{4} + \frac{3 \pi A_{1} A_{2}^{2} \alpha_{3}}{2} - \pi A_{1} \omega^{2} + \pi A_{1} \omega_{0}^{2}$




```python
proj_eqn['eq2']
```




$\displaystyle 3 \pi A_{0}^{2} A_{2} \alpha_{3} + \frac{3 \pi A_{0} A_{1}^{2} \alpha_{3}}{2} + \frac{3 \pi A_{1}^{2} A_{2} \alpha_{3}}{2} + \frac{3 \pi A_{2}^{3} \alpha_{3}}{4} - 4 \pi A_{2} \omega^{2} + \pi A_{2} \omega_{0}^{2}$



if $A_1$ is small neglect terms with $A_0^2, A_2^2, A_0A_2, A_1^2A_0, A_0^3, A_1^2A_2, \cdots$

This is because there is the assumption of 
$$
\begin{gather*}
    A_1 > A_i \phantom{---} i = 0, 2, 3, \cdots
\end{gather*}
$$

equation 0 tells us that $A_0$ is proportional to $A_1$ squared - We have a constant term multiplying $A_0$ (in addition to others) and no $A_1$ terms but only $A_1^2$ terms.


```python
deglist = list(map(int, sp.degree_list(proj_eqn['eq0'], gens=Ai)))
A1_deg = deglist[1]
removing_combos = []
Ai_ = tuple(filter(lambda A: A != Ai[1], Ai))

# A1 ** n -> n can be up to O(A1**4) which is square of A1**2
keep_terms = []
for i in range(A1_deg+1, A1_deg**2+1):
    temp = []
    for j in range(i):
        temp.append(Ai[1])
    keep_terms.append(tuple(temp))

temp = list(combinations_with_replacement(Ai_, A1_deg))
for t in temp:
    if t not in keep_terms:
        removing_combos.append(t)
    
for i in range(A1_deg+1, max(deglist)+1):
    temp = list(combinations_with_replacement(Ai, i))
    for t in temp:
        if t not in keep_terms:
            removing_combos.append(t)
        
eq_copy = proj_eqn['eq0']
for term in removing_combos:
    eq_copy = eq_copy.subs(np.prod(term), 0)

eq_copy = sp.Eq(eq_copy, 0)

eq0_new = eq_copy
eq0_ordered = eq0_new
epsilon = sp.Symbol('epsilon')
for A in Ai_:
    eq0_ordered = eq0_ordered.subs(A, epsilon**A1_deg * A)
eq0_ordered = sp.collect(eq0_ordered.subs(Ai[1], epsilon * Ai[1]).lhs, epsilon)
eq0_ordered = sp.Eq(eq0_ordered, 0)

eq_copy = sp.Eq(Ai[0], sp.solve(eq_copy, Ai[0])[0])
eq_copy
```




$\displaystyle A_{0} = 0$



equation 2 tells $A_2$ is also proportional to $A_1$ squared 


```python
deglist = list(map(int, sp.degree_list(proj_eqn['eq2'], gens=Ai)))
A1_deg = deglist[1]
Ai_ = tuple(filter(lambda A: A != Ai[1], Ai))

for i in range(A1_deg+1, max(deglist)+1):
    temp = list(combinations_with_replacement(Ai, i))
    for t in temp:
        if t not in removing_combos:
            if t not in keep_terms:
                removing_combos.append(t)

eq_copy = proj_eqn['eq2']
for term in removing_combos:
    eq_copy = eq_copy.subs(np.prod(term), 0)

eq_copy = sp.Eq(eq_copy, 0)

eq2_new = eq_copy
eq2_ordered = eq2_new
for A in Ai_:
    eq2_ordered = eq2_ordered.subs(A, epsilon**A1_deg * A)
eq2_ordered = sp.collect(eq2_ordered.subs(Ai[1], epsilon * Ai[1]).lhs, epsilon)
eq2_ordered = sp.Eq(eq2_ordered, 0)

eq_copy = sp.Eq(Ai[2], sp.solve(eq_copy, Ai[2])[0])
eq_copy
```




$\displaystyle A_{2} = 0$




```python
eq_copy = proj_eqn['eq1']
for term in removing_combos:
    eq_copy = eq_copy.subs(np.prod(term), 0)

eq_copy = sp.Eq(eq_copy, 0)
eq1_new = eq_copy
eq1_ordered = eq1_new
for A in Ai_:
    eq1_ordered = eq1_ordered.subs(A, epsilon**A1_deg * A)
eq1_ordered = sp.collect(eq1_ordered.subs(Ai[1], epsilon * Ai[1]).lhs, epsilon)
eq1_ordered = sp.Eq(eq1_ordered, 0)
```

To facilitate the solution, we order the coefficients using epsilon and then ignore all terms quadratic and higher in $A_0$ and $A_2$ - this means we keep terms up to tand including $O(\epsilon^3)$ -> small amplitude assumption

$A_1$ is ordered at $\epsilon$ and $A_0$ and $A_2$ are order at $\epsilon^2$


```python
# Eq0
eq0_ordered
```




$\displaystyle 2 \pi A_{0} \epsilon^{2} \omega_{0}^{2} = 0$




```python
# Eq1 
eq1_ordered
```




$\displaystyle \frac{3 \pi A_{1}^{3} \alpha_{3} \epsilon^{3}}{4} + \epsilon \left(- \pi A_{1} \omega^{2} + \pi A_{1} \omega_{0}^{2}\right) = 0$




```python
# Eq2
eq2_ordered
```




$\displaystyle \epsilon^{2} \left(- 4 \pi A_{2} \omega^{2} + \pi A_{2} \omega_{0}^{2}\right) = 0$




```python
A0_sol = sp.solve(eq0_ordered, Ai[0])[0]
A0_sol
```




$\displaystyle 0$




```python
A2_sol = sp.solve(eq2_ordered, Ai[2])[0]
A2_sol
```




$\displaystyle 0$




```python
eq1_ordered_update = eq1_ordered.subs([
    (Ai[0], A0_sol), (Ai[2], A2_sol)
])
omegaSol = sp.solve(eq1_ordered_update, w)
```


```python
omegaSol[0]
```




$\displaystyle - \frac{\sqrt{3 A_{1}^{2} \alpha_{3} \epsilon^{2} + 4 \omega_{0}^{2}}}{2}$




```python
omegaSol[1]
```




$\displaystyle \frac{\sqrt{3 A_{1}^{2} \alpha_{3} \epsilon^{2} + 4 \omega_{0}^{2}}}{2}$



We are keeping up to $O(\epsilon^4)$ only - use this to get our final expression for $\omega$


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


```python
omega1 = taylor(omegaSol[0], epsilon, 4)
omega1
```




$\displaystyle - \omega_{0} - \frac{3 A_{1}^{2} \alpha_{3} \epsilon^{2}}{8 \omega_{0}} + O\left(\epsilon^{4}\right)$




```python
omega2 = taylor(omegaSol[1], epsilon, 4)
omega2
```




$\displaystyle \omega_{0} + \frac{3 A_{1}^{2} \alpha_{3} \epsilon^{2}}{8 \omega_{0}} + O\left(\epsilon^{4}\right)$




```python
omega = sp.simplify(omega2)
omega
```




$\displaystyle \omega_{0} + \frac{3 A_{1}^{2} \alpha_{3} \epsilon^{2}}{8 \omega_{0}} + O\left(\epsilon^{4}\right)$




```python
t = sp.Symbol('t')
phi_temp = omega.removeO() * t
phi_temp
```




$\displaystyle t \left(\frac{3 A_{1}^{2} \alpha_{3} \epsilon^{2}}{8 \omega_{0}} + \omega_{0}\right)$




```python
x = xh.subs([
    (Ai[0], A0_sol), (Ai[2], A2_sol), (phi, omega.removeO()*t)
])
x = x.subs(w, omega.removeO())
x
```




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




```python
initial_amplitude = sp.simplify(x.subs([
    (t, 0), (epsilon, 1)
]))
initial_amplitude
```




$\displaystyle A_{1}$




```python
ic = float(sp.simplify(initial_amplitude.subs([
    (Ai[1], 1)
])))
ic
```




    1.0




```python
# We plot the phase plane of this system 

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


```python
# System
def sys(t, x):
    return [x[1], -x[0] - x[0]**3]

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

Using the initial condition found, we get response with period $2\pi$ not $4\pi$


```python
x0 = [[0.1, 0.0]]

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

DEplot(sys, (0, 2*np.pi, 1000), x0, x1, x2, {})
```

An attempt to use a much bigger initial condition


```python
x0 = [[3.0, 0.0]]

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

DEplot(sys, (0, 2*np.pi, 1000), x0, x1, x2, {})
```

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