Example 5: Parametrically excited Duffing oscillator

MMS example on the Duffing oscillator subject to parametric forcing. This configuration was studied by Nayfeh and Mook in Nonlinear Oscillations (1995), section 5.7.3.

System description

Illustration of a parametrically forced Duffing oscillator through the time-varying stiffness \(\omega_0^2 + 2 F \cos(\omega t)\).

The system’s equation is

\[\ddot{x} + c \dot{x} + \omega_0^2 x + \gamma x^3 = -2 x 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.

A parametric response at twice the oscillator’s frequency is sought so the frequency is set to

\[\omega = 2\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 \tilde{c}\) indicates that damping is weak,

• \(F = \epsilon \tilde{F}\) indicates that forcing is weak,

• \(\gamma = \epsilon \tilde{\gamma}\) indicates that nonlinearities are weak.

Code description

The script below allows to

• Construct the dynamical system.

• Apply the MMS to the system,

• Evaluate the MMS results at steady state,

• Compute the forced response and the backbone curve,

• Evaluate the stability of the computed forced solution.

# -*- 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)
fF = -2*x # Parametric forcing
dyn = MMS.Dynamical_system(t, x, Eq, omega0, fF=fF, F=F)

# Initialisation of the MMS sytem
eps            = symbols(r"\epsilon", real=True, positive=True) # Small parameter epsilon
Ne             = 1      # Order of the expansions
omega_ref      = omega0 # Reference frequency
ratio_omegaMMS = 2      # Look for a solution around n*omega_ref

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

mms = MMS.Multiple_scales_system(dyn, eps, Ne, omega_ref, sub_scaling, ratio_omegaMMS=ratio_omegaMMS)

# Application of the MMS
mms.apply_MMS(rewrite_polar="all")

# Evaluation at steady state
ss = MMS.Steady_state(mms)

# Solve the evolution equations for a given dof
solve_dof = 0 # dof to solve for
ss.solve_bbc(solve_dof=solve_dof, c=param_scaled[-1])
ss.solve_forced(solve_dof=solve_dof)

# Stability analysis
ss.stability_analysis(coord="polar", eigenvalues=True)

# %%

Validation

The notations used in the above code are related to those in Nonlinear Oscillations as follows:

Link between notations from Nonlinear Oscillations and this document.

Current Notation	Nonlinear Oscillations	Description
\(u\)	\(x\)	Oscillator’s dof
\(\omega_0\)	\(\omega\)	Natural frequency
\(\omega\)	1	Forcing frequency
\(2\omega_0+\epsilon\sigma\)	\(\omega+\epsilon\sigma\)	Forcing frequency definition
\(F\)	\(1\)	Forcing amplitude
\(\tilde{c}\)	\(2\mu\)	Damping coefficient (scaled)
\(\tilde{\gamma}\)	\(\alpha\)	Nonlinear coefficient (scaled)
\(x_0^{\textrm{h}}\)	[no name]	Leading order homogeneous solution
\(a_0\)	\(a\)	Oscillator’s amplitude
\(2\beta_0\)	\(\psi\)	Oscillator’s autonomous phase
\(\textrm{D}_1(\bullet)\)	\((\bullet)'\)	Slow time derivative
\(\lambda_i\)	\(\lambda_i\)	Eigenvalues of the Jacobian matrix

Note that the above table implies that the detuning \(\sigma\) in Nonlinear Oscillations’s notations is half that of the current notation. Section 5.7.3 from Nonlinear Oscillations gives the following results

\[\begin{split}\begin{cases} u(t) & = a \cos(t - \frac{1}{2}\psi) + \mathcal{O}(\epsilon), \\ a' & = - \dfrac{a}{2 \omega} \sin \psi - \mu a, \\ a\psi' & = 2\sigma a - \dfrac{a}{\omega} \cos \psi - \dfrac{3 \alpha}{4 \omega} a^3, \\ \lambda_{1,2} & = - \mu \pm \sqrt{\mu^2 + \frac{3}{4} \alpha a^2 \cos \psi}, \quad \cos \psi = 2 \sigma \omega - \frac{3}{4} \alpha a^2. \end{cases}\end{split}\]

In the current notations, this is equivalent to

\[\begin{split}\begin{cases} x_0^{\textrm{h}}(t) & = a_0 \cos(\frac{\omega}{2} t - \beta_0), \\ \textrm{D}_1 a_0 & = - \dfrac{\tilde{c} a_{0}}{2} - \dfrac{\tilde{F} a_{0} \sin{\left(2 \beta_{0} \right)}}{2 \omega_{0}}, \\ a_0 \textrm{D}_1 \beta_0 & = \dfrac{\sigma a_{0}}{2} - \dfrac{\tilde{F} a_{0} \cos{\left(2 \beta_{0} \right)}}{2 \omega_{0}} - \dfrac{3 \tilde{\gamma} a_{0}^{3}}{8 \omega_{0}}, \\ \lambda_{1,2} & = \epsilon \dfrac{1}{2} \dfrac{ \left(- \omega_{0} \tilde{c} \pm \sqrt{\omega_{0}^{2} \tilde{c}^{2} + 3 \omega_{0} \sigma \tilde{\gamma} a_{0}^{2} - \frac{9}{4} \tilde{\gamma}^{2} a_{0}^{4}}\right)}{ \omega_{0}}, \end{cases}\end{split}\]

which are the outputs for ss.sol.x[0][0].simplify(), ss.sol.faO[0][1], ss.sol.fbetaO[0][1] and ss.stab.eigvals, respectively. Note the \(\epsilon\) factor for the eigenvalues. It is related to the eigenvalues being computed from the reconstituted modulation equations rather than just the \(1^{\textrm{st}\) order slow time modulation equations.