Example 3: Coupled Duffings in 1:1 internal resonance

MMS example on coupled Duffing oscillators in 1:1 internal resonance subject to harmonic forcing.

System description

Nonlinear system.

Illustration of two forced, nonlinearly coupled Duffing oscillators in 1:1 internal resonance.

The system’s equations are

\[\begin{split}\begin{cases} \ddot{x}_{0} + \omega_{0}^{2} x_{0} + c_{0} \dot{x}_{0} + \gamma_{0} x_{0}^{3} + \gamma_{01} x_{0} x_{1}^{2} & = F \cos(\omega t), \\ \ddot{x}_{1} + \omega_{0}^{2} x_{1} + c_{1} \dot{x}_{1} + \gamma_{1} x_{1}^{3} + \gamma_{01} x_{0}^{2} x_{1} & = F \cos(\omega t), \end{cases}\end{split}\]

where

  • \(x_0,\; x_1\) are the oscillators’ coordinates,

  • \(t\) is the time,

  • \(\dot{(\bullet)} = \mathrm{d}(\bullet)/\mathrm{d}t\) is a time derivative,

  • \(c_0,\; c_1\) are the linear viscous damping coefficients,

  • \(\omega_0\) is the oscillators’ natural frequency,

  • \(\gamma_0,\; \gamma_1,\; \gamma_{01}\) are nonlinear coefficients,

  • \(F\) is the forcing amplitude,

  • \(\omega\) is the forcing frequency.

A response around \(\omega_0\) is sought so the frequency is set to

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

where

  • \(\epsilon\) is a small parameter involved in the MMS,

  • \(\sigma\) is the detuning wrt the oscillators’ frequency \(\omega_0\).

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

  • \(c_i = \epsilon \tilde{c}_i\) indicates that damping is weak,

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

  • \(\gamma_i = \epsilon \tilde{\gamma}_i, \; i\in\{0, 1, 01\}\) 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 backbone curve and forced response when only oscillator 1 responds,

  • Evaluate the stability of the forced solution,

  • Compute the coupled-mode backbone curve.

 1# -*- coding: utf-8 -*-
 2
 3#%% Imports and initialisation
 4from sympy import symbols, Function, solve, cos
 5from sympy.physics.vector.printing import init_vprinting
 6init_vprinting(use_latex=True, forecolor='White') # Initialise latex printing 
 7from oscilate import MMS
 8
 9# Parameters and variables
10omega0, F, c0, c1       = symbols(r'\omega_0, F, c_0, c_1', real=True, positive=True)
11gamma0, gamma1, gamma01 = symbols(r'\gamma_0, \gamma_1, \gamma_{01}', real=True)
12t  = symbols('t')
13x0 = Function(r'x_0', real=True)(t)
14x1 = Function(r'x_1', real=True)(t)
15
16# Dynamical system
17Eq0 = x0.diff(t,2) + omega0**2*x0 + gamma0*x0**3 + gamma01*x0*x1**2 + c0*x0.diff(t)
18Eq1 = x1.diff(t,2) + omega0**2*x1 + gamma1*x1**3 + gamma01*x0**2*x1 + c1*x1.diff(t)
19dyn = MMS.Dynamical_system(t, [x0, x1], [Eq0, Eq1], [omega0, omega0], fF=[0, 1], F=F)
20
21# Initialisation of the MMS sytem
22eps            = symbols(r"\epsilon", real=True, positive=True) # Small parameter epsilon
23Ne             = 1      # Order of the expansions
24omega_ref      = omega0 # Reference frequency
25ratio_omegaMMS = 1      # Look for a solution around omega_ref
26
27param_to_scale = (gamma0, gamma1, gamma01, F , c0, c1)
28scaling        = (1     , 1     , 1      , Ne, Ne, Ne)
29param_scaled, sub_scaling = MMS.scale_parameters(param_to_scale, scaling, eps)
30
31mms = MMS.Multiple_scales_system(dyn, eps, Ne, omega_ref, sub_scaling, ratio_omegaMMS=ratio_omegaMMS)
32
33# Application of the MMS
34mms.apply_MMS()
35
36# Evaluation at steady state
37ss = MMS.Steady_state(mms)
38
39# Solve the evolution equations for a given dof
40solve_dof = 1 # dof to solve for
41ss.solve_bbc(solve_dof=solve_dof, c=param_scaled[-2:])
42ss.solve_forced(solve_dof=solve_dof)
43
44# Stability analysis of the 1dof solution
45ss.stability_analysis(coord="cartesian", eigenvalues=False, rewrite_polar=True)
46
47# Computation of the coupled free solution
48print("Manual computation of the coupled free solution")
49Dbeta         = symbols(r"\Delta\beta", real=True) # Phase difference beta0-beta0
50sub_phase     = [(ss.coord.beta[1], ss.coord.beta[0]+Dbeta)] # Substitute the phases by the phase difference
51fa_free       = [fai.subs(ss.sub.sub_free+sub_phase) for fai in ss.sol.fa]
52Dbeta_sol     = solve(fa_free[0], Dbeta) # Solution in terms of phase difference
53chi           = symbols(r"\chi", real=True) # +/- symbol introduced to represent cos(2*Dbeta_sol)
54sub_cos2Dbeta = [(cos(2*Dbeta), chi)]
55fbeta_free    = [fbetai.subs(ss.sub.sub_free + sub_phase + sub_cos2Dbeta) for fbetai in ss.sol.fbeta]
56a02_sol       = solve(fbeta_free[0], ss.coord.a[0]**2)[1]
57sig_bbc       = solve(fbeta_free[1].subs(ss.coord.a[0]**2, a02_sol), mms.sigma)[0].subs(chi**2,1)
58
59
60# %%