Duffing oscillator - complex form

MMS example on the Duffing oscillator subject to a direct harmonic forcing. This configuration was studied by Nayfeh and Mook [2], sections 4.1 and 4.1.1, but here it is solved using a complex form of the dynamical system, as in [5].

System description

Nonlinear system. — Illustration of a forced Duffing oscillator.

The system’s equation is

\[\ddot{x} + c \dot{x} + \omega_0^2 x + \gamma x^3 = F \cos(\omega t),\]

where

\(x\) is the oscillator’s coordinate,
\(t\) is the time,
\(\dot{(\bullet)} = \mathrm{d}(\bullet)/\mathrm{d}t\) is a time derivative,
\(c\) is the linear viscous damping coefficient,
\(\omega_0\) is the oscillator’s natural frequency,
\(\gamma\) is the nonlinear coefficient,
\(F\) is the forcing amplitude,
\(\omega\) is the forcing frequency.

This system of equations is transformed through the introduction of the complex coordinates

\[\begin{split}\begin{cases} x & = z + \bar{z}, \\ \dot{x} & = \textrm{j} \omega_0 (z - \bar{z}), \end{cases}\end{split}\]

or equivalently

\[\begin{split}\begin{cases} z & = \frac{1}{2} (x - \frac{\textrm{j}}{\omega_0} \dot{x}), \\ \bar{z} & = \frac{1}{2} (x + \frac{\textrm{j}}{\omega_0} \dot{x}), \end{cases}\end{split}\]

resulting in the complex first order dynamical system (on \(z(t)\))

\[\dot{z} = \textrm{j} \omega z + \frac{\textrm{j} \gamma \left(z + \overline{z}\right)^{3}}{2 \omega} - \frac{c z}{2} + \frac{c \overline{z}}{2}.\]

A direct response is sought so the frequency (either backbone curve frequency or forcing frequency) is close from the oscillator’s, such that

\[\omega = \omega_0 + \epsilon \sigma,\]

where

\(\epsilon\) is a small parameter involved in the MMS,
\(\sigma\) is the detuning wrt the reference frequency \(\omega_0\).

The parameters are then scaled to indicate how weak they are:

\(c = \epsilon^{N_e} \tilde{c}\) indicates that damping is weak (of order \(N_e\)),
\(F = \epsilon^{N_e} \tilde{F}\) indicates that forcing is weak (of order \(N_e\)),
\(\gamma = \epsilon \tilde{\gamma}\) indicates that nonlinearities are weak (of order 1).

Note that \(N_e\) is the order up to which the solutions are sought and the time scales are constructed. It is therefore chosen here to scale the damping and forcing at maximum order.

Code description

The script below allows to

Construct the dynamical system in oscillator and complex forms
Apply the MMS to the system up to order \(N_e=5\) using the oscillator and complex forms
Evaluate the MMS results at steady state
Compute the backbone curve and the associated mapping from the normal coordinate’s amplitude to the displacement (initial) coordinate. These results are in accordance with those from the normal form theory [6].

# -*- coding: utf-8 -*-

#%% Imports and initialisation
from sympy import symbols, Function
from sympy.physics.vector.printing import init_vprinting, vlatex
init_vprinting(use_latex=True, forecolor='White') # Initialise latex printing 
from oscilate import MMS

# Parameters and variables
omega0, F, c = symbols(r'\omega_0, F, c', real=True,positive=True)
gamma        = symbols(r'\gamma',real=True)
t            = symbols('t')
x            = Function(r'x', real=True)(t)

# Dynamical system
Eq = x.diff(t,2) + omega0**2*x + gamma*x**3 + c*x.diff(t)
dyn = MMS.Dynamical_system(t, x, Eq, omega0, F=F)
dyn.complex_form() # Construct the complex form of the system

# Initialisation of the MMS sytem
eps            = symbols(r"\epsilon", real=True, positive=True) # Small parameter epsilon
Ne             = 5      # Order of the expansions
omega_ref      = omega0 # Reference frequency
ratio_omegaMMS = 1      # Look for responses around the reference frequency (Default behaviour)

param_to_scale = (gamma, F , c )
scaling        = (1    , Ne, Ne)
param_scaled, sub_scaling = MMS.scale_parameters(param_to_scale, scaling, eps)

# Complex form of the MMS
mms = MMS.Multiple_scales_complex(dyn, eps, Ne, omega_ref, sub_scaling)
mms.apply_MMS(orders_polar="all") # Application of the MMS
ss = MMS.Steady_state(mms) # Steady state analysis
ss.solve_bbc(solve_dof=0, c=param_scaled[-1]) # Backbone curve computation

# Oscillator form of the MMS
mms_o = MMS.Multiple_scales_oscillator(dyn, eps, Ne, omega_ref, sub_scaling)
mms_o.apply_MMS(orders_polar="all") # Application of the MMS
ss_o = MMS.Steady_state(mms_o) # Steady state analysis
ss_o.solve_bbc(solve_dof=0, c=param_scaled[-1]) # Backbone curve computation

# Oscillator form of the MMS - order 1 
mms_o1 = MMS.Multiple_scales_oscillator(dyn, eps, 1, omega_ref, sub_scaling)
mms_o1.apply_MMS(orders_polar="all") # Application of the MMS
ss_o1 = MMS.Steady_state(mms_o1) # Steady state analysis
ss_o1.solve_bbc(solve_dof=0, c=param_scaled[-1]) # Backbone curve computation

# Plot the steady state results
# -----------------------------

# Set parameters' numerical values
import numpy as np
param = [(omega0, 1),
         (gamma, 0.5),
         (ss.coord.a[0], np.linspace(1e-10, 1.6, 1000))]

# Frequency response
param_FRC = param + [(dyn.forcing.F, 1e-2)]
BBC    = MMS.visualisation.Backbone_curve(mms, ss, dyn, param_FRC)
BBC_o  = MMS.visualisation.Backbone_curve(mms_o, ss_o, dyn, param_FRC)
BBC_o1 = MMS.visualisation.Backbone_curve(mms_o1, ss_o1, dyn, param_FRC)

# Plot
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(BBC.omega, BBC.xmax, label=f"complex form - Ne={Ne}")
ax.plot(BBC_o.omega, BBC_o.xmax, label=f"oscillator form - Ne={Ne}")
ax.plot(BBC_o1.omega, BBC_o1.xmax, label=f"oscillator form - Ne=1")
ax.axvline(BBC.omegaMMS, c="k")
ax.set_xlabel(r"$\omega_{\text{nl}}$")
ax.set_ylabel(r"$max(x(t))$")
ax.legend()
ax.set_ymargin(0)

# %%

Plot outputs

The plot outputs shown below are generated from the code above.

The figure below displays the backbone curve as the peak oscillator amplitude as a function of the oscillation frequency. Three curves are shown, associated to the results obtained with the complex (blue) and oscillator (orange) forms, and the later if pushing the expansions only up to \(N_e=1\) (green). One can see how these results deviate at large amplitudes.