Example 1: Duffing oscillator

MMS example on the Duffing oscillator subject to harmonic forcing. This configuration was studied by Nayfeh and Mook in Nonlinear Oscillations (1995), sections 4.1 and 4.1.1.

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.

A direct response is sought so the frequency (either backbone curve frequency or forcing frequency) is close from the oscilltor’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.

  • Apply the MMS to the system up to order \(N_e = 3\). The MMS results can be visualized in LaTeX in an IPython-based interactive environment (e.g., VS Code’s Python Interactive Window or Jupyter Notebook). For example, to visualize the modulation function \(f_\beta(a, \beta)\), use

    mms.sol.fbeta  # Display the symbolic expression for :math:`f_\beta(a, \beta)`

  • Evaluate the MMS results at steady state. The solutions are stored in ss.sol (e.g., ss.sol.fa, ss.sol.fbeta).

  • Compute the forced response and the backbone curve. These results are also stored in ss.sol (e.g., ss.sol.sigma, ss.sol.F, ss.sol.omega_bbc).

  • Evaluate the stability of the computed forced solution. The stability results are stored in ss.stab (e.g., ss.stab.eigvals, ss.stab.bif_a).

  • Evaluate the steady-state results for given numerical values of the parameters and plot the results using methods like ss.plot_FRC() or ss.plot_ARC().

 1# -*- coding: utf-8 -*-
 2
 3#%% Imports and initialisation
 4from sympy import symbols, Function
 5from sympy.physics.vector.printing import init_vprinting, vlatex
 6init_vprinting(use_latex=True, forecolor='White') # Initialise latex printing 
 7from oscilate import MMS
 8
 9# Parameters and variables
10omega0, F, c = symbols(r'\omega_0, F, c', real=True,positive=True)
11gamma        = symbols(r'\gamma',real=True)
12t            = symbols('t')
13x            = Function(r'x', real=True)(t)
14
15# Dynamical system
16Eq = x.diff(t,2) + omega0**2*x + gamma*x**3 + c*x.diff(t)
17dyn = MMS.Dynamical_system(t, x, Eq, omega0, F=F)
18
19# Initialisation of the MMS sytem
20eps            = symbols(r"\epsilon", real=True, positive=True) # Small parameter epsilon
21Ne             = 3      # Order of the expansions
22omega_ref      = omega0 # Reference frequency
23ratio_omegaMMS = 1      # Look for responses around the reference frequency (Default behaviour)
24
25param_to_scale = (gamma, F , c )
26scaling        = (1    , Ne, Ne)
27param_scaled, sub_scaling = MMS.scale_parameters(param_to_scale, scaling, eps)
28
29mms = MMS.Multiple_scales_system(dyn, eps, Ne, omega_ref, sub_scaling)
30
31# Application of the MMS
32mms.apply_MMS()
33
34# Evaluation at steady state
35ss = MMS.Steady_state(mms)
36
37# Solve the evolution equations for a given dof
38solve_dof = 0 # dof to solve for
39ss.solve_forced(solve_dof=solve_dof)
40ss.solve_bbc(solve_dof=solve_dof, c=param_scaled[-1])
41
42# Stability analysis
43ss.stability_analysis(coord="polar", eigenvalues=True, bifurcation_curves=True)
44
45# Plot the steady state results
46# -----------------------------
47import numpy as np
48
49# Set parameters' numerical values
50a0 = np.linspace(1e-10, 1.2, 1000)
51
52dic_numpy = dict(
53    omega0 = (omega0, 1),
54    c      = (c, 1e-2),
55    gamma  = (gamma, 0.2),
56    a      = (ss.coord.a[0], a0),
57    )
58
59F_val     = 1e-2
60omega_val = 1.05
61
62# Compute and plot the frequency-response curves (FRC)
63dic_FRC = dic_numpy | dict(F=(dyn.forcing.F, F_val)) # Parameters for the FRC
64FRC     = MMS.visualisation.numpise_FRC(mms, ss, dyn, dic_FRC)
65kwargs  = dict(phase_name=vlatex(ss.sol.cos_phase[0].args[0]),  # Plot parameters
66               amp_name=vlatex(ss.coord.a[0]))
67MMS.visualisation.plot_FRC(FRC, **kwargs)
68
69# Compute and plot the amplitude-response curves (ARC)
70dic_ARC = dic_numpy | dict(omega=(mms.omega, omega_val)) # Parameters for the ARC
71ARC     = MMS.visualisation.numpise_ARC(mms, ss, dyn, dic_ARC)
72kwargs  = dict(phase_name=vlatex(ss.sol.cos_phase[0].args[0]), # Plot parameters
73               amp_name=vlatex(ss.coord.a[0]))
74MMS.visualisation.plot_ARC(ARC, **kwargs)
75
76# %%