A wide array of contemporary problems can be represented by nonlinear ordinary differential equations with solutions that can be represented by Fourier Series:
These can all be observed in a quick literature search on 'Harmonic Balance'.
Mousai is designed to solve the cases where $\mathbf{u}(t)$ can be represented by
\begin{equation} \mathbf{u}(t)=A_0 + \sum_{n=0}^{\infty} A_{n}\sin(n \times \omega t) + B_n\cos(n \times \omega t) \end{equation}In practicality this is represented in the form:
\begin{equation} \mathbf{u}(t)=\sum_{n=-\infty}^{\infty} \mathbf{U}_{n}e^{j(n \times \omega t)} \end{equation}where $U_n = \frac{A_n}{2}-j\frac{B_n}{2}$ because the math is far simpler, if less intuitive for the analyst. This also works better computationally and is closely related to the form of a discrete Fourier transform- a process much more suitable to computation.
The solution is presumed to be of the form:
\begin{equation} \mathbf{z}(t)=\sum_{n=-\infty}^{\infty} \mathbf{Z}_{n}e^{j(n \times \omega t)} \end{equation}In practice this holds true even if the governing equation is nonlinear.
When the model is linear, we can use superposition to solve for one term at a time. This solution is very well known (See , and ).
For a linear system, the state equation reduces to
\begin{equation} \dot{\mathbf{z}}(t)=A \mathbf{z}(t) + B \mathbf{u}(t) \end{equation}where
\begin{equation} A = \frac{\partial\mathbf{f}(\mathbf{z}(t),\mathbf{u}(t))}{\partial\mathbf{z}(t)},\qquad B = \frac{\partial\mathbf{f}(\mathbf{z}(t),\mathbf{u}(t))}{\partial\mathbf{u}(t)}\end{equation}Taking the Fourier transform, this is
\begin{equation}j\omega\mathbf{Z}(\omega)=A\mathbf{Z}(\omega)+B\mathbf{U}(\omega)\end{equation}The solution is:
\begin{equation} \label{eq:linsoln} \mathbf{Z}(\omega) = \left(Ij\omega-A\right)^{-1}B\mathbf{U}(\omega) \end{equation}where the magnitudes and phases of the elements of $\mathbf{Z}$ provide the amplitudes and phases of the harmonic response of each state at the frequency $\omega$.
Thus, each value $\mathbf{Z}_n(n\omega)$ is
\begin{equation} \mathbf{Z}_n(n\omega) = \left(Ijn\omega-A\right)^{-1}B\mathbf{U}_n(n\omega) \end{equation}===========
When a model is nonlinear closed-form solutions do not typically exist– numerical methods must then be used to estimate the solution.
Numerically finding the oscillatory response, after dissipation of the transient response, requires long time marching.
- Without energy dissipation term in the model this is not feasible because transients do not attenuate with time. When appears possible it is because numerical dissipation enables it. However, the numerical dissipation phenomenon isn't real- the solution represents that of a system with an equivalent amount of dissipation, not the non-dissipative system being modeled.
- With dissipation, simulations require tens, hundreds, or thousands of cycles, therefore tens of thousands of times steps may be necessary to simulate until a time at which the transient isn't noticeable.
- Time marching can be inaccurate or unstable. Numerical energy generation or loss can substantially impact the solution.For a linear system in the frequency domain this is
\begin{equation}j\omega\mathbf{Z}(\omega)=\mathbf{f}(\mathbf{Z}(\omega),\mathbf{U}(\omega))\end{equation}\begin{equation}j\omega\mathbf{Z}(\omega)=A\mathbf{Z}(\omega)+B\mathbf{U}(\omega)\end{equation}where
\begin{equation}A = \frac{\partial \mathbf{f}(\mathbf{Z}(\omega),\mathbf{U}(\omega))}{\partial\mathbf{Z}(\omega)},\qquad B = \frac{\partial \mathbf{f}(\mathbf{Z}(\omega),\mathbf{U}(\omega))}{\partial\mathbf{U}(\omega)}\end{equation}are constant matrices.
The solution is:
\begin{equation}\mathbf{Z}(\omega) = \left(Ij\omega-A\right)^{-1}B\mathbf{U}(\omega)\end{equation}where the magnitudes and phases of the elements of $\mathbf{Z}$ provide the amplitudes and phases of the harmonic response of each state at the frequency $\omega$.
This is actually a function call to a Finite Element Package, CFD, Matlab function, - whatever your solver uses to get derivatives
We can also find $\dot{\mathbf{z}}(t)$ from the derivative of the Fourier Series:
In [4]:
# Define our function (Python)
def duff_osc_ss(x, params):
omega = params['omega']
t = params['cur_time']
xd = np.array([[x[1]],
[-x[0] - 0.1 * x[0]**3 - 0.1 * x[1] + 1 * sin(omega * t)]])
return xd
In [5]:
# Arguments are name of derivative function, number of states, driving frequency,
# form of the equation, and number of harmonics
t, x, e, amps, phases = ms.hb_time(duff_osc_ss, num_variables=2, omega=.1,
eqform='first_order', num_harmonics=5)
print('Displacement amplitude is ', amps[0])
print('Velocity amplitude is ', amps[1])
In [6]:
def pltcont():
time, xc = ms.time_history(t, x)
disp_plot, _ = plt.plot(time, xc.T[:, 0], t,
x.T[:, 0], '*b', label='Displacement')
vel_plot, _ = plt.plot(time, xc.T[:, 1], 'r',
t, x.T[:, 1], '*r', label='Velocity')
plt.legend(handles=[disp_plot, vel_plot])
plt.xlabel('Time (sec)')
plt.title('Response of Duffing Oscillator at 0.0159 rad/sec')
plt.ylabel('Response')
plt.legend
plt.grid(True)
In [7]:
fig=plt.figure()
ax=fig.add_subplot(111)
time, xc = ms.time_history(t, x)
disp_plot, _ = ax.plot(time, xc.T[:, 0], t,
x.T[:, 0], '*b', label='Displacement')
vel_plot, _ = ax.plot(time, xc.T[:, 1], 'r',
t, x.T[:, 1], '*r', label='Velocity')
ax.legend(handles=[disp_plot, vel_plot])
ax.set_xlabel('Time (sec)')
ax.set_title('Response of Duffing Oscillator at 0.0159 rad/sec')
ax.set_ylabel('Response')
ax.legend
ax.grid(True)
In [19]:
pltcont()# abbreviated plotting function
In [20]:
time, xc = ms.time_history(t, x)
disp_plot, _ = plt.plot(time, xc.T[:, 0], t,
x.T[:, 0], '*b', label='Displacement')
vel_plot, _ = plt.plot(time, xc.T[:, 1], 'r',
t, x.T[:, 1], '*r', label='Velocity')
plt.legend(handles=[disp_plot, vel_plot])
plt.xlabel('Time (sec)')
plt.title('Response of Duffing Oscillator at 0.0159 rad/sec')
plt.ylabel('Response')
plt.legend
plt.grid(True)
In [21]:
omega = np.arange(0, 3, 1 / 200) + 1 / 200
amp = sp.zeros_like(omega)
amp[:] = np.nan
t, x, e, amps, phases = ms.hb_time(duff_osc_ss, num_variables=2,
omega=1 / 200, eqform='first_order', num_harmonics=1)
for i, freq in enumerate(omega):
# Here we try to obtain solutions, but if they don't work,
# we ignore them by inserting `np.nan` values.
x = x - sp.average(x)
try:
t, x, e, amps, phases =
ms.hb_time(duff_osc_ss, x0=x,
omega=freq, eqform='first_order', num_harmonics=1)
amp[i] = amps[0]
except:
amp[i] = np.nan
if np.isnan(amp[i]):
break
plt.plot(omega, amp)
In [10]:
omegal = np.arange(3, .03, -1 / 200) + 1 / 200
ampl = sp.zeros_like(omegal)
ampl[:] = np.nan
t, x, e, amps, phases = ms.hb_time(duff_osc_ss, num_variables=2,
omega=3, eqform='first_order', num_harmonics=1)
for i, freq in enumerate(omegal):
# Here we try to obtain solutions, but if they don't work,
# we ignore them by inserting `np.nan` values.
x = x - np.average(x)
try:
t, x, e, amps, phases = ms.hb_time(duff_osc_ss, x0=x,
omega=freq, eqform='first_order', num_harmonics=1)
ampl[i] = amps[0]
except:
ampl[i] = np.nan
if np.isnan(ampl[i]):
break
In [8]:
plt.plot(omega,amp, label='Up sweep')
plt.plot(omegal,ampl, label='Down sweep')
plt.legend()
plt.title('Amplitude versus frequency for Duffing Oscillator')
plt.xlabel('Driving frequency $\\omega$')
plt.ylabel('Amplitude')
plt.grid()
In [11]:
def two_dof_demo(x, params):
omega = params['omega']
t = params['cur_time']
force_amplitude = params['force_amplitude']
alpha = params['alpha']
# The following could call an external code to obtain the state derivatives
xd = np.array([[x[1]],
[-2 * x[0] - alpha * x[0]**3 + x[2]],
[x[3]],
[-2 * x[2] + x[0]]] + force_amplitude * np.sin(omega * t))
return xd
In [12]:
parameters = {'force_amplitude': 0.2}
parameters['alpha'] = 0.4
t, x, e, amps, phases = ms.hb_time(two_dof_demo, num_variables=4,
omega=1.2, eqform='first_order', params=parameters)
amps
Out[12]:
In [13]:
alpha = np.linspace(-1, .45, 2000)
amp = np.zeros_like(alpha)
for i, alphai in enumerate(alpha):
parameters['alpha'] = alphai
t, x, e, amps, phases = ms.hb_time(two_dof_demo, num_variables=4, omega=1.2,
eqform='first_order', params=parameters)
amp[i] = amps[0]
In [12]:
plt.plot(alpha,amp)
plt.title('Amplitude of $x_1$ versus $\\alpha$')
plt.ylabel('Amplitude of $x_1$')
plt.xlabel('$\\alpha$')
plt.grid()
In [16]:
def two_dof_coulomb(x, params):
omega = params['omega']
t = params['cur_time']
force_amplitude = params['force_amplitude']
mu = params['mu']
# The following could call an external code to obtain the state derivatives
xd = np.array([[x[1]],
[-2 * x[0] - mu * np.abs(x[1]) + x[2]],
[x[3]],
[-2 * x[2] + x[0]]] + force_amplitude * np.sin(omega * t))
return xd
In [17]:
parameters = {'force_amplitude': 0.2}
parameters['mu'] = 0.1
t, x, e, amps, phases = ms.hb_time(two_dof_coulomb, num_variables=4,
omega=1.2, eqform='first_order', params=parameters)
amps
Out[17]:
In [18]:
mu = np.linspace(0, 1.0, 200)
amp = np.zeros_like(mu)
for i, mui in enumerate(mu):
parameters['mu'] = mui
t, x, e, amps, phases = ms.hb_time(two_dof_coulomb, num_variables=4, omega=1.2,
eqform='first_order', num_harmonics=3, params=parameters)
amp[i] = amps[0]
In [32]:
plt.plot(mu,amp)
plt.title('Amplitude of $x_1$ versus $\\mu$')
plt.ylabel('Amplitude of $x_1$')
plt.xlabel('$\\mu$')
plt.grid()
In [25]:
out = ms.hb_time(lambda x, v,
params: np.array([[-x - .1 * x**3 - .1 * v + 1 *
sin(params['omega'] * params['cur_time'])]]),
num_variables=1, omega=.7, num_harmonics=1)
out[3][0]
Out[25]:
OK - that's a bit obtuse. I wouldn't do that normally, but Mousai can.
pip install mousai
In [ ]: