Modules

MMS Module

dyn_sys

Started on Tue Feb 15 17:25:59 2022

@author: Vincent MAHE

Analyse systems of coupled nonlinear equations using the Method of Multiple Scales (MMS). This sub-module defines the dynamical system.

class oscilate.MMS.dyn_sys.Dynamical_system(t, x, Eq, omegas, **kwargs)

Bases: object

The dynamical system studied. See Dynamical systems of interest for a detailed description of the dynamical system.

Parameters:

• t (sympy.Symbol) – time \(t\).

• x (sympy.Function or list of sympy.Function) – Unknown(s) of the problem.

• Eq (sympy.Expr or list of sympy.Expr) – System’s equations without forcing, which can be defined separately (see parameters F and fF). Eq is the unforced system of equations describing the system’s dynamics.

• omegas (sympy.Symbol or list of sympy.Symbol) – The natural frequency of each oscillator.

• F (sympy.Symbol or 0, optional) – Forcing amplitude \(F\). Default is 0.

• fF (sympy.Expr or list of sympy.Expr, optional) – For each oscillator, specify the coefficient multiplying the forcing terms in the equation. It can be used to define parametric forcing. Typically, if the forcing is \(x F \cos(\omega t)\), then fF = x. Default is a list of 1, so the forcing is direct.

class oscilate.MMS.dyn_sys.Forcing(F, fF)

Bases: object

Define the forcing on the system as

• A forcing amplitude F,

• Forcing coefficients fF, used to introduce parametric forcing or simply weight the harmonic forcing.

For the \(i^\textrm{th}\) oscillator, denoting fF[i] as \(f_{F,i}(\boldsymbol{x}(t), \dot{\boldsymbol{x}}(t), \ddot{\boldsymbol{x}}(t))\), the forcing term on that oscillator is \(f_{F,i} F \cos(\omega t)\).

mms

Started on Tue Feb 15 17:25:59 2022

@author: Vincent MAHE

Analyse systems of coupled nonlinear equations using the Method of Multiple Scales (MMS). This sub-module defines the multiple scales system from the dynamical one, and the application of the MMS.

oscilate.MMS.mms.Chain_rule_d2fdt2(f, tS, eps)

Apply the chain rule to express second order time derivatives in terms of the time scales’ derivatives.

Parameters:

• f (sympy.Function) – Time scales-dependent function \(f(\boldsymbol{t})\).

• tS (list) – Time scales \(\boldsymbol{t}^\intercal = [t_0, \cdots, t_{N_e}]\).

• eps (sympy.Symbol) – Small parameter \(\epsilon\).

Returns:

d2fdt2 – \(\mathrm{d}^2 f/ \mathrm{d}t^2\) expressed in terms of the time scales.

Return type:

sympy.Function

Notes

Consider a time scales-dependent function \(f(t_0, t_1, ...)\), where \(t_0\) is the fast time and \(t_1, ...\) are the slow times. The Chain Rule is applied to give the expression of \(\mathrm{d}^2 f/ \mathrm{d}t^2\) in terms of the time scales.

oscilate.MMS.mms.Chain_rule_dfdt(f, tS, eps)

Apply the chain rule to express first order time derivatives in terms of the time scales’ derivatives.

Parameters:

• f (sympy.Function) – Time scales-dependent function \(f(\boldsymbol{t})\).

• tS (list) – Time scales \(\boldsymbol{t}^\intercal = [t_0, \cdots, t_{N_e}]\).

• eps (sympy.Symbol) – Small parameter \(\epsilon\).

Returns:

dfdt – \(\mathrm{d} f/ \mathrm{d}t\) expressed in terms of the time scales.

Return type:

sympy.Function

Notes

Consider a time scales-dependent function \(f(t_0, t_1, ...)\), where \(t_0\) is the fast time and \(t_1, ...\) are the slow times. The Chain Rule is applied to give the expression of \(\mathrm{d} f/ \mathrm{d}t\) in terms of the time scales.

class oscilate.MMS.mms.Coord_MMS(mms)

Bases: object

The coordinates used in the MMS.

class oscilate.MMS.mms.Forcing_MMS(F, f_order, fF, forcing_term)

Bases: object

Define the forcing on the system as

• A forcing amplitude F

• A scaling order f_order for the forcing

• Forcing coefficients fF

• Forcing terms (direct or parametric) forcing_term

class oscilate.MMS.mms.Multiple_scales_system(dynamical_system, eps, Ne, omega_ref, sub_scaling, ratio_omegaMMS=1, eps_pow_0=0, **kwargs)

Bases: object

The multiple scales system. See Application of the MMS for a detailed description of the dynamical system.

Parameters:

• dynamical_system (Dynamical_system) – The dynamical system.

• eps (sympy.Symbol) – Small perturbation parameter \(\epsilon\).

• Ne (int) – Truncation order of the asymptotic series and order of the slowest time scale.

• omega_ref (sympy.Symbol) – Reference frequency \(\omega_{\textrm{ref}}\) of the MMS. Not necessarily the frequency around which the MMS is going to be applied, see ratio_omegaMMS.

• sub_scaling (list of tuples) – Substitutions to do to scale the equations. Links small parameters to their scaled counterpart through \(\epsilon\).

• ratio_omegaMMS (int or sympy.Rational, optional) –

  Specify the frequency omegaMMS around which the MMS is going to be applied in terms of \(\omega_{\textrm{ref}}\). Denoting ratio_omegaMMS as \(r_{\textrm{MMS}}\), this means that

  \[\omega_{\textrm{MMS}} = r_{\textrm{MMS}} \omega_{\textrm{ref}}.\]

  Use ratio_omegaMMS=Rational(p,q) for

  \[q \omega_{\textrm{MMS}} = p \omega_{\textrm{ref}}\]

  to get better-looking results than the float \(p/q\). Default is 1.

• eps_pow_0 (int, optional) –

  Order of the leading order term in the asymptotic series of each oscillators’ response. For the \(i^{\textrm{th}}\) oscillator and denoting eps_pow_0 as \(\lambda_0\), this means that

  \[x_i = \epsilon^{\lambda_0} x_{i,0} + \epsilon^{\lambda_0+1} x_{i,1} + \cdots.\]

  Default is 0.

• ratio_omega_osc (list of int or sympy.Rational or None, optional) –

  Specify the natural frequencies of the oscillators \(\omega_i\) in terms of the reference frequency \(\omega_{\textrm{ref}}\). Denoting ratio_omega_osc[i] as \(r_i\), this means that

  \[\omega_i \approx r_i \omega_{\textrm{ref}}.\]

  Use ratio_omega_osc[i]=Rational(p,q) for

  \[q \omega_{i} \approx p \omega_{\textrm{ref}}\]

  to get better-looking results than the float \(p/q\). Default is None for each oscillator, so the \(\omega_i\) are arbitrary and there are no internal resonances. Detuning can be introduced through the detunings keyword argument.

• detunings (list of sympy.Symbol or int, optional) –

  The detuning of each oscillator. Denoting detunings[i] as \(\delta_i\), this means that

  \[\omega_i = r_i \omega_{\textrm{ref}} + \delta_i.\]

  Default is 0 for each oscillator.

Notes

Description of the method of multiple scales.

apply_MMS(rewrite_polar=0)

Apply the MMS.

Parameters:

rewrite_polar (str, int or list of int, optional) – The orders at which the solutions will be rewritten in polar form. See sol_x_polar().

Notes

The application of the MMS is operated by successively calling the following methods.

1. system_t0(): An equivalent system is written in terms of the fast time scale \(t_0\). This introduces the temporary unknowns \(\tilde{x}_{i,j}(t_0)\) and allows the use of dsolve().

2. sol_order_0(): Leading order solutions \(x_{i,0}(\boldsymbol{t})\) are defined.

3. secular_analysis(): The leading order solutions are introduced in the equations and the secular terms at each order are identified. Cancelling those secular terms is a condition for bounded solutions. It leads to a system of modulation equations governing the slow time evolution of the complex amplitude of the homogeneous leading order solutions. Each equation takes the form

   \[\textrm{D}_{j} A_i(\boldsymbol{t}_s) = f_{A_i}^{(j)}(\boldsymbol{A}, t_1).\]

   After cancelling the secular terms the higher order equations are solved successively to express the higher order solutions \(x_{i,j}(\boldsymbol{t}),\; j>0\) in terms of the leading order ones.

4. autonomous_phases(): The phase coordinates are changed from \(\phi_i(\boldsymbol{t}_s)\) to \(\beta_i(\boldsymbol{t}_s)\) to cancel the slow time \(t_1\) in the secular terms. This will be used afterwards to obtain an autonomous system.

5. modulation_equations(): The secular conditions are split into real and imaginary parts, polar coordinates are used and the autonomous phases are introduced, resulting in an autonomous system of modulation equations on polar coordinates. Equations come by two, one representing the amplitude modulation while the other represents the phase’s, such that

   \[\begin{split}\begin{cases} \textrm{D}_{j} a_i(\boldsymbol{t}_s) & = f_{a_i}^{(j)}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_i \textrm{D}_{j} \beta_i(\boldsymbol{t}_s) & = f_{\beta_i}^{(j)}(\boldsymbol{a}, \boldsymbol{\beta}). \end{cases}\end{split}\]

   This is the key result of the MMS. The modulations on each time scale can be combined to reintroduce the physical time, resulting in a system of the form

   \[\begin{split}\begin{cases} \dfrac{\textrm{d}}{dt} a_i(t) & = f_{a_i}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_i \dfrac{\textrm{d}}{dt} \beta_i(t) & = f_{\beta_i}(\boldsymbol{a}, \boldsymbol{\beta}). \end{cases}\end{split}\]

   This is known as the reconstitution method.

6. sol_x_polar(): The leading and higher order solutions are rewritten in terms of polar coordinates using \(\cos\) and \(\sin\) functions.

asymptotic_series(dynamical_system, eps_pow_0=0)

Define the asymptotic series.

Notes

The series expansion for oscillator \(i\) (and for a leading order term \(\epsilon^0 = 1\)) takes the form (see Application of the MMS)

\[x_i(t) = x_{i,0}(t) + \epsilon x_{i,1}(t) + \epsilon^2 x_{i,2}(t) + \cdots + \epsilon^{N_e} x_{i,N_e}(t) + \mathcal{O}(\epsilon^{N_e+1}).\]

On top of introducing the terms of the asymptotic series, this function prepares substitutions from

1. dof \(x_i(t)\) to temporary \(t\)-dependent asymptotic terms \(x_{i,j}(t)\), such that

   \[x_i(t) = x_{i,0}(t) + \epsilon x_{i,1}(t) + \epsilon^2 x_{i,2}(t) + \cdots + \epsilon^{N_e} x_{i,N_e}(t),\]

2. Temporary \(x_{i,j}(t)\) to the time scales-dependent terms \(x_{i,j}(\boldsymbol{t})\), such that

   \[x_i(\boldsymbol{t}) = x_{i,0}(\boldsymbol{t}) + \epsilon x_{i,1}(\boldsymbol{t}) + \epsilon^2 x_{i,2}(\boldsymbol{t}) + \cdots + \epsilon^{N_e} x_{i,N_e}(\boldsymbol{t}).\]

autonomous_phases()

Define phase coordinates that render an autonomous system.

Notes

Define new phase coordinates \(\beta_i\) to transform nonautonomous equations into autonomous ones. The \(\beta_i\) are defined as

\[\beta_i = \frac{r_i}{r_{\textrm{MMS}}} \sigma t_1 - \phi_i,\]

where we recall that \(\omega = r_{\textrm{MMS}} \omega_{\textrm{ref}} + \epsilon \sigma\) and \(\omega_{i,0} = r_i \omega_{\textrm{ref}}\). See details on this choice in Application of the MMS.

compute_EqO(dynamical_system)

Compute the system of equations for each oscillator at each order of \(\epsilon\). This system is described in Application of the MMS.

The output EqO is a list of lists:

• The \(1^{\text{st}}\) level lists are associated to the equations for each oscillator,

• The \(2^{\text{nd}}\) level lists are associated to the orders of \(\epsilon\) from the lowest to the highest order.

Parameters:

dynamical_system (Dynamical_system) – The dynamical system.

find_harmonics()

Determine the harmonics contained in the MMS solutions.

Returns:

harmonics – list of the harmonics appearing in the MMS solutions.

Return type:

list

forcing_MMS(dynamical_system)

Rewrite the forcing terms.

Parameters:

dynamical_system (Dynamical_system) – The dynamical system.

Notes

The initial forcing terms are

\[f_{F,i}(\boldsymbol{x}(t), \dot{\boldsymbol{x}}(t), \ddot{\boldsymbol{x}}(t)) F \cos(\omega t), \quad i=1,...,N.\]

Rewritting them involves

1. Replacing the \(x_i(t)\) by their series expansions written in terms of time scales,

2. Scaling the forcing and the parameters in \(f_{F,i}\) if any,

3. Truncating terms whose order is larger than the largest order retained in the MMS,

4. Rewrite \(\cos(\omega t)\) as

   \[\cos(\omega t) = \frac{1}{2} e^{\mathrm{j}(\omega_{\textrm{MMS}} + \epsilon \sigma)t} + cc = \frac{1}{2} e^{\mathrm{j}(\omega_{\textrm{MMS}}t_0 + \sigma t_1)} + cc.\]

modulation_equations()

Derive the modulation equations of the polar coordinates system.

Notes

Derive the modulation equations of the polar coordinates system (defined in polar_coordinates() and autonomous_phases()) from the secular conditions. For oscillator \(i\), these are defined as

\[\begin{split}\begin{cases} \dfrac{\textrm{d} a_i}{\textrm{d} t} & = f_{a_i}(\boldsymbol{a}, \boldsymbol{\beta}), \\ a_i \dfrac{\textrm{d} \beta_i}{\textrm{d} t} & = f_{\beta_i}(\boldsymbol{a}, \boldsymbol{\beta}), \end{cases}\end{split}\]

where \(\boldsymbol{a}\) and \(\boldsymbol{\beta}\) are vectors containing the polar amplitudes and phases.

The aim here is to compute all the \(f_{a_i}\) and \(f_{\beta_i}\). This is done by:

1. Introducing polar coordinates in the secular terms

2. Splitting the real and imaginary parts of the (complex) secular terms

3. Using the autonomous phase coordinates

4. Collecting the terms governing the slow amplitude and phase dynamics.

oscillators_frequencies()

Gives the expression of every oscillator frequency in terms of the reference frequency, possibly with a detuning.

Notes

For the \(i^\textrm{th}\) oscillator, this leads to

\[\omega_i = r_i \omega_{\textrm{ref}} + \delta_i .\]

An associated first-order natural frequency \(\omega_{i,0}\) is defined by neglecting the detuning \(\delta_i\), which is at least of order \(\epsilon\), resulting in

\[\omega_{i,0} = r_i \omega_{\textrm{ref}}.\]

polar_coordinates()

Define the polar coordinates.

Notes

Define polar coordinates such that, for oscillator \(i\), the complex amplitude of the homogeneous leading order solution is defined as

\[A_i(\boldsymbol{t}_s) = \frac{1}{2} a_i(\boldsymbol{t}_s) e^{\textrm{j} \phi_i(\boldsymbol{t}_s)},\]

where \(a_i(\boldsymbol{t}_s)\) and \(\phi_i(\boldsymbol{t}_s)\) are the solution’s amplitude and phase, respectively.

reconstitution()

Use the reconstitution method to combine the modulation equations at each order. This reconstitution is based on the chain rule relation

\[\dfrac{\textrm{d}(\bullet)}{\textrm{d}t} = \sum_{i=0}^{N_e} \epsilon^{i} \mathrm{D}_i (\bullet) + \mathcal{O}(\epsilon^{N_e+1}).\]

Note that some MMS approaches do not apply this reconstitution step.

secular_analysis()

Identify the secular terms in the MMS equations. This allows to:

1. Compute the modulation equations of the complex amplitudes \(A_i(\boldsymbol{t}_s)\), coming from the elimination of the secular terms,

2. Derive nonsecular MMS equations, i.e. MMS equations with the secular terms cancelled,

3. Use the nonsecular equations to express the higher order solutions \(x_{i,j}(\boldsymbol{t}),\; j>0\) in terms of the \(\boldsymbol{A}(\boldsymbol{t}_s)\).

sol_higher_order(EqO_t0, xO_t0, io, ix)

Compute the higher order solutions \(x_{i,j}(\boldsymbol{t}_s),\; j>0\).

Parameters:

• EqO_t0 (list of list of sympy.Expr) – The MMS equations at each order and for each oscillator written with \(t_0\) as the only independent variable.

• xO_t0 (list of list of sympy.Function) – Oscillators’ response at each order written in terms of \(t_0\) only, \(\tilde{x}_{i,j}(t_0)\).

• io (int) – The order of \(\epsilon\).

• ix (int) – The dof number.

sol_order_0()

Compute the leading-order solutions for each oscillator.

Notes

For oscillator \(i\), the homogeneous solution takes the general form

\[x_{i,0}^{\textrm{h}}(\boldsymbol{t}) = A_i(\boldsymbol{t}_s) e^{\textrm{j} \omega_{i,0} t_0} + cc,\]

where \(A_i(\boldsymbol{t}_s)\) is an unknown complex amplitude.

If the oscillator is subject to hard forcing (i.e. forcing appears at leading order), then the particular solution

\[x_{i,0}^{\textrm{p}}(t_0, t_1) = B_i e^{\textrm{j} \omega t} + cc = B_i e^{\textrm{j} (\omega_{\textrm{MMS}} t_0 + \sigma t_1)} + cc\]

is also taken into account. \(B_i\) is a time-independent function of the forcing parameters.

sol_x_polar(rewrite_polar=0)

Write the solutions using the polar coordinates and \(\cos\) and \(\sin\) functions.

Parameters:

rewrite_polar (str or int or list of int, optional) – The orders at which the solutions will be rewritten in polar form. If "all", then all solution orders will be rewritten. If int, then only a single order will be rewritten. If list of int, then the listed orders will be rewritten. Default is 0, so only the leading order solution will be rewritten.

system_t0()

Rewrite the system with the fast time scale as the only independent variable.

Notes

This is a trick to use dsolve(), which only accepts functions of 1 variable, to solve higher order equations. The higher order equations are rewritten in terms of temporary coordinates \(\tilde{x}_{i,j}(t_0)\) in place of \(x_{i,j}(\boldsymbol{t})\), with \(i,j\) denoting the oscillator number and \(\epsilon\) order, respectively. This is equivalent to temporary considering that \(\boldsymbol{A}(\boldsymbol{t}_s)\) does not depend on the slow times, which is of no consequence as there are no slow time derivatives appearing in the higher order equations at this stage. Indeed, they were either substituted using the complex modulation equations, or they disappeared when eliminating the secular terms.

time_scales()

Define the time scales.

Notes

The time scales are defined as (see Application of the MMS)

\[t_0 = t, \quad t_1 = \epsilon t, \quad \cdots, \quad t_{N_e} = \epsilon^{N_e} t.\]

Substitutions from the physical time \(t\) to the time scales \(t_i, \; i=0, ..., N_e\) are also prepared. Note that \(t_0\) is refered-to as the fast time as it captures the oscillations. Other time scales are refered-to as slow time as they capture the modulations.

class oscilate.MMS.mms.Sol_MMS

Bases: object

Solutions obtained when applying the MMS.

class oscilate.MMS.mms.Substitutions_MMS(sub_t, sub_xO_t, sub_x, sub_scaling, sub_omega, sub_sigma)

Bases: object

Substitutions used in the MMS.

oscilate.MMS.mms.cartesian_to_polar(y, sub_polar, sub_phase=None)

Rewrites an expression or a matrix y from cartesian to polar coordinates.

Parameters:

• y (sympy.Expr or sympy.Matrix) – A sympy expression or matrix written in cartesian coordinates.

• sub_polar (list) – A list of substitutions to perform to go from cartesian to polar coordinates.

• sub_phase (list, optional) – Additional substitutions to try and get rid of phases, so that only the amplitude remains in the expression. Default is None.

Returns:

yp – The initial expression or matrix written in polar coordinates.

Return type:

sympy.Expr or sympy.Matrix

oscilate.MMS.mms.rescale(expr, mms)

Rescales a scaled expression.

Parameters:

• expr (sympy.Expr) – An unscaled expression, i.e. an expression appearing at some order of \(\epsilon\).

• mms (Multiple_scales_system) – The mms system, containing substitutions to scale an expression.

Returns:

expr_scaled – The scaled expression.

Return type:

sympy.Expr

oscilate.MMS.mms.scale_parameters(param, scaling, eps)

Scale parameters with the scaling parameter \(\epsilon\).

Parameters:

• param (list of sympy.Symbol and/or sympy.Function) – Unscaled parameters.

• scaling (list of int or float) – The scaling for each parameter.

• eps (sympy.Symbol) – Small parameter \(\epsilon\).

Returns:

• param_scaled (list of sympy.Symbol and/or sympy.Function) – Scaled parameters.

• sub_scaling (list of 2 lists of tuple) – Substitutions from scaled to unscaled parameters and vice-versa.

  • \(1^{\text{st}}\) list: The substitutions to do to introduce the scaled parameters in an expression.

  • \(2^{\text{nd}}\) list: The substitutions to do to reintroduce the unscaled parameters in a scaled expression.

Notes

For a given parameter \(p\) and a scaling order \(\lambda\), the associated scaled parameter \(\tilde{p}\) is

\[p = \epsilon^{\lambda} \tilde{p} .\]

steady_state

Started on Tue Feb 15 17:25:59 2022

@author: Vincent MAHE

Analyse systems of coupled nonlinear equations using the Method of Multiple Scales (MMS). This sub-module defines the steady state system from the multiple scales one, and the functions for steady state evaluation.

class oscilate.MMS.steady_state.Coord_SS

Bases: object

The coordinates used in the steady state analysis.

class oscilate.MMS.steady_state.Forcing_SS(mms)

Bases: object

Define the forcing on the system.

class oscilate.MMS.steady_state.Sol_SS(ss, mms)

Bases: object

Solutions obtained when evaluating at steady state.

class oscilate.MMS.steady_state.Stab_SS

Bases: object

Stability analysis parameters and outputs.

class oscilate.MMS.steady_state.Steady_state(mms)

Bases: object

Steady state analysis of the multiple scales system. See Evaluation at steady state for a detailed description of the dynamical system.

Parameters:

mms (Multiple_scales_system) – The multiple scales system.

Notes

Description of the steady state analysis.

Jacobian_cartesian()

Compute the Jacobian of the modulation equations systems expressed in cartesian coordinates (see stability_analysis()).

Returns:

J – Jacobian of the cartesian system.

Return type:

sympy.Matrix

Jacobian_polar()

Compute the Jacobian of the modulation equations systems expressed in polar coordinates (see stability_analysis()).

Returns:

J – Jacobian of the polar system.

Return type:

sympy.Matrix

additional_cartesian_substitutions(fp, fq)

Reformulate the already-existing substitutions from polar to cartesian to try and substitute leftover polar terms.

Parameters:

• fp (list of sympy.Expr) – Evolution functions for the cartesian coordinates \(p_i\). There are polar coordinates remaining.

• fq (list of sympy.Expr) – Evolution functions for the cartesian coordinates \(q_i\). There are polar coordinates remaining.

Returns:

• fp (list of sympy.Expr) – Evolution functions for the cartesian coordinates \(p_i\). Additional substitutions were performed to get rid of polar coordinates.

• fq (list of sympy.Expr) – Evolution functions for the cartesian coordinates \(q_i\). Additional substitutions were performed to get rid of polar coordinates.

bifurcation_curves(detJ, var_a=False, var_sig=True, solver=<function solve_poly2>)

Compute bifurcation curves, defined by the simple bifurcation points of the slow time system obtained for any forcing frequency and amplitude.

Parameters:

• detJ (sympy.Expr) – The determinant of the matrix.

• var_a (bool, optional) – Consider the \(i^{\textrm{th}}\) oscillator’s amplitude \(a_i\) as the variable and find the bifurcation curve as an expression for \(a_i\). detJ is rarely a quadratic polynomial in \(a_i\), so this can rarely be computed easily. Default is False.

• var_sig (bool, optional) – Consider the detuning \(\sigma\) as the variable and find the bifurcation curve as an expression for \(\sigma\). detJ is often a quadratic polynomial in \(\sigma\), so this can often be computed. Default is True.

• solver (function, optional) – The solver to use to compute the bifurcation curves. Available are solver called as solve(expr, x), which solve expr=0 for x. solve() can be used but is sometimes slow. Default is solve_poly2().

Returns:

• bif_a (list) – The bifurcation curves for \(a_i^2\).

• bif_sig (list) – The bifurcation curves for \(\sigma\).

Notes

Simple bifurcations are associated to real eigenvalues of the Jacobian matrix crossing the imaginary axis, hence passing through 0. As the determinant is the product of the eigenvalues, the zeros of the determinant give the bifurcation points of the system.

cartesian_coordinates()

Define cartesian coordinates from the polar ones.

Notes

The homogeneous leading order solution for oscillator \(i\) expressed in polar coordinates takes the form

\[\begin{split}\begin{split} x_{i,0}^{\textrm{h}}(t) & = a_i \cos\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t - \beta_i \right), \\ & = a_i \cos(\beta_i) \cos\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right) + a_i \sin(\beta_i) \sin\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right) \end{split}\end{split}\]

The polar coordinates are defined as

\[\begin{split}\begin{cases} p_i = a_i \cos (\beta_i), \\ q_i = a_i \sin (\beta_i), \end{cases}\end{split}\]

such that the leading order solution can be written as

\[x_{i,0}^{\textrm{h}}(t) = p_i \cos\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right) + q_i \sin\left(\frac{r_i}{r_{\textrm{MMS}}}\omega t\right).\]

check_cartesian_substitutions(a, beta, fp, fq)

Check if substitutions from polar to cartesian coordinates are complete.

Parameters:

• a (list of sympy.Symbol) – Amplitudes of the leading order solutions.

• beta (list of sympy.Symbol) – Phases of the leading order solutions.

• fp (list of sympy.Expr) – Evolution functions for the cartesian coordinates \(p_i\).

• fq (list of sympy.Expr) – Evolution functions for the cartesian coordinates \(q_i\).

Returns:

substitution_OK – True if substitutions are complete. False otherwise.

Return type:

bool

eigenvalues(A, detA=None, trA=None)

Computes the eigenvalues of a matrix \(\textrm{A}\).

Parameters:

• A (sympy.Matrix) – The matrix whose eigenvalues are to be computed.

• detA (sympy.Expr) – Determinant of A. Default is None.

• trA (sympy.Expr) – Trace of A. Default is None.

Returns:

eigvals – The eigenvalues of \(\textrm{A}\).

Return type:

list

modulation_equations_SS(mms)

Evaluate the modulation equations at steady state (polar system).

modulation_equations_cartesian()

Compute the modulation equations of the cartesian coordinates system.

Notes

Write the modulation equations using the cartesian coordinates (defined in cartesian_coordinates()). For oscillator \(i\), this results in

\[\begin{split}\begin{cases} \dfrac{\textrm{d} p_i}{\textrm{d} t} & = f_{p_i}(\boldsymbol{p}, \boldsymbol{q}), \\ \dfrac{\textrm{d} q_i}{\textrm{d} t} & = f_{q_i}(\boldsymbol{p}, \boldsymbol{q}), \end{cases}\end{split}\]

where \(\boldsymbol{p}(t)\) and \(\boldsymbol{q}(t)\) are vectors containing the cartesian coordinates.

polar_coordinates_SS(mms)

Introduce time-independent amplitudes and phases (polar coordinates).

solve_F()

Solve the forced response in terms of the forcing amplitude \(F\). Returns \(F(a_i)\)

solve_a()

Solve the forced response in terms of the oscillator’s amplitude. For readability, the output actually returned is \(a_i^2(\sigma, F)\).

solve_bbc(c=[], solve_dof=None)

Find the backbone curve (bbc) of a given oscillator with the other oscillators’ amplitude set to 0.

Parameters:

• c (list, optional) – Damping terms. They will be set to 0 to compute the backbone curve. Note that these are the scaled damping terms. Default is [].

• solve_dof (None or int, oprtional) – The oscillator number to solve for. If None, no oscillator is solved for. Default is None.

Notes

The backbone curve describes the frequency of free oscillations as a function of the oscillator’s amplitude. In the presence of small damping (as in the case in the MMS) the frequency of free oscillations is close from the resonance frequency. The backbone curve can therefore be interpreted as the backbone of the forced response.

The backbone curve of oscillator \(i\) typically takes the form

\[\omega_{\textrm{bbc}}^{(i)} = k\omega_{i} + f_{\textrm{bbc}}^{(i)} (a_i),\]

where \(k=1,\; k<1,\; k>1\) are associated to direct, superharmonic and subharmonic responses, respectively.

solve_forced(solve_dof=None)

Solve the forced response of an oscillator.

Parameters:

solve_dof (None or int, optional) – The oscillator to solve for. If None, no oscillator is solved for. Default is None.

Notes

Find the steady state solution for a given oscillator with the other oscillators’ amplitude set to 0.

To do so, one must choose an oscillator to chose for, say oscillator \(i\). Then, the following methods are called:

1. substitution_solve_dof(): Set the other oscillators’ amplitude to 0, i.e. \(a_j = 0 \; \forall j \neq i\).

2. solve_phase(): express the oscillator’s phase \(\beta_i\) as a function of its amplitude \(a_i\).

3. solve_sigma(): find the expression of \(\sigma(a_i)\).

4. solve_a(): find the expression of \(a_i (\sigma, F)\).

5. solve_F(): find the expression of \(F(a_i)\).

solve_phase()

Find solutions for the oscillator’s phase \(\beta_i\). The solutions actually returned are \(\sin(\beta_i)\) and \(\cos(\beta_i)\).

solve_sigma()

Solve the forced response in terms of the detuning \(\sigma\). Returns \(\sigma(a_i)\). It is recalled that \(\omega = \omega_{\textrm{MMS}} + \epsilon \sigma\).

stability_analysis(coord='cartesian', rewrite_polar=False, eigenvalues=False, bifurcation_curves=False, trace_curves=False, analyse_blocks=False, kwargs_bif={})

Evaluate the stability of a steady state solution. See Stability analysis for a detailed description of the dynamical system.

Parameters:

• coord (str, optional) – Either "cartesian" or "polar". Specifies the coordinates to use for the stability analysis. "cartesian" is recommended as it prevents divisions by 0, which occur when at least one of the oscillator has a 0 ampliutude. Default is "cartesian".

• rewrite_polar (bool, optional) – Rewrite the Jacobian’s determinant and trace in polar coordinates (if computed using cartesian ones). This is time consuming and the current back substitutions from cartesian to polar coordinates are not always sufficient. Default is False.

• eigenvalues (bool, optional) – Compute the eigenvalues of the Jacobian. Default is False.

• bifurcation_curves (bool, optional) – Compute the bifurcation curves. Default is False.

• trace_curves (bool, optional) – Compute the trace curves, i.e. the zeros of the Jacobian’s trace. Default is False.

• analyse_blocks (bool, optional) – Analyse the diagonal blocks of the Jacobian rather than the Jacobian itself. This is relevant if the Jacobian is block-diagonal.

• kwargs_bif (dict, optional) – Passed to bifurcation_curves() Default is dict().

substitution_solve_dof(solve_dof)

Set every oscillator amplitude to 0 except the one to solve for.

Notes

If one wants to solve for \(a_i\), then the system is evaluated for \(a_j=0, \; \forall j \neq i\).

trace_curves(trJ, var_a=False, var_sig=True, solver=<function solve_poly2>)

Compute curves capturing the variations of the trace of the Jacobian. For 2 by 2 Jacobians, a negative trace with a positive determinant indicates a stable response.

Parameters:

• trJ (sympy.Expr) – The trace of the matrix.

• var_a (bool, optional) – Consider the \(i^{\textrm{th}}\) oscillator’s amplitude \(a_i\) as the variable and find the trace curve as an expression for \(a_i\). Default is False.

• var_sig (bool, optional) – Consider the detuning \(\sigma\) as the variable and find the trace curve as an expression for \(\sigma\). Default is True.

• solver (function, optional) – The solver to use to compute the trace curves. Available are solver called as solve(expr, x), which solve expr=0 for x. solve() can be used but is sometimes slow. Default is solve_poly2().

Returns:

• tr_a (list) – The trace curves for \(a_i^2\).

• tr_sig (list) – The trace curves for \(\sigma\).

class oscilate.MMS.steady_state.Substitutions_SS(mms)

Bases: object

Substitutions used in the steady state evaluations.

visualisation

Started on Tue Feb 15 17:25:59 2022

@author: Vincent MAHE

Analyse systems of coupled nonlinear equations using the Method of Multiple Scales (MMS). This sub-module evaluates the symbolic expressions for given numerical parameters and allows to plot the resulting numerical results.

oscilate.MMS.visualisation.numpise_ARC(mms, ss, dyn, param)

Evaluate the amplitude-response curves at given numerical values. This transforms the sympy expressions to numpy arrays.

Parameters:

• mms (Multiple_scales_system) – The MMS object.

• ss (Steady_state) – The MMS results evaluated at steady state.

• dyn (Dynamical_system) – The initial dynamical system.

• param (dict) –

  A dictionnary whose values are tuples with 2 elements:

  1. The sympy symbol of a parameter,

  2. The numerical value(s) taken by that parameter.

  The key of the amplitude vector must be "a". The key of the angular frequency must be "omega".

Returns:

ARC – The amplitude-response curves data.

Return type:

dict

oscilate.MMS.visualisation.numpise_FRC(mms, ss, dyn, param, bbc=True, forced=True, bif=True)

Evaluate the frequency-response and bifurcation curves at given numerical values. This transforms the sympy expressions to numpy arrays.

Parameters:

• mms (Multiple_scales_system) – The MMS object.

• ss (Steady_state) – The MMS results evaluated at steady state.

• dyn (Dynamical_system) – The initial dynamical system.

• param (dict) –

  A dictionary whose values are tuples with 2 elements:

  1. The sympy symbol of a parameter,

  2. The numerical value(s) taken by that parameter.

  The key of the amplitude vector must be "a". The key of the forcing amplitude must be "F".

• bbc (bool, optional) – Evaluate the backbone curve. Default is True.

• forced (bool, optional) – Evaluate the forced response. Default is True.

• bif (bool, optional) – Evaluate the bifurcation curves. Default is True.

Returns:

FRC – The frequency-response curves data.

Return type:

dict

oscilate.MMS.visualisation.numpise_F_ARC(mms, ss, param)

Numpise the forced response’s forcing amplitude \(F\).

Parameters:

• mms (Multiple_scales_system)

• ss (Steady_state)

• param (dict) – See sympy_to_numpy().

Returns:

F – Numpised forced response’s forcing amplitude.

Return type:

numpy.ndarray

oscilate.MMS.visualisation.numpise_omega_FRC(mms, ss, param)

Numpise the forced response’s frequency \(\omega\).

Parameters:

• mms (Multiple_scales_system)

• ss (Steady_state)

• param (dict) – See sympy_to_numpy().

Returns:

omega – Numpised forced response’s frequency.

Return type:

numpy.ndarray

oscilate.MMS.visualisation.numpise_omega_bbc(mms, ss, param)

Numpise the backbone curve’s frequency \(\omega_{\textrm{bbc}}\).

Parameters:

• mms (Multiple_scales_system)

• ss (Steady_state)

• param (dict) – See sympy_to_numpy().

Returns:

omega_bbc – Numpised backbone curve’s frequency.

Return type:

numpy.ndarray

oscilate.MMS.visualisation.numpise_omega_bif(mms, ss, param)

Numpise the bifurcation curves’ frequency \(\omega_{\textrm{bif}}\).

Parameters:

• mms (Multiple_scales_system)

• ss (Steady_state)

• param (dict) – See sympy_to_numpy().

Returns:

omega_bif – Numpised bifurcation curves’ frequency.

Return type:

list of numpy.ndarray

oscilate.MMS.visualisation.numpise_phase(mms, ss, dyn, param, omega, F)

Numpise the phase \(\beta_i\).

Parameters:

• mms (Multiple_scales_system)

• ss (Steady_state)

• dyn (Dynamical_system)

• param (dict) – See sympy_to_numpy().

• omega (numpy.ndarray or list of numpy.ndarray) – The frequency array.

• F (numpy.ndarray) – The forcing amplitude array.

Returns:

phase – Numpised phase.

Return type:

list of numpy.ndarray

oscilate.MMS.visualisation.plot_ARC(ARC, **kwargs)

Plots the amplitude-response curves (ARC), both forcing amplitude-amplitude and forcing amplitude-phase.

Parameters:

ARC (dict) – Dictionary containing the amplitude response curves.

Returns:

• fig1 (Figure) – The amplitude plot \(a(F)\).

• fig2 (Figure) – The phase plot \(\beta(F)\).

oscilate.MMS.visualisation.plot_FRC(FRC, **kwargs)

Plots the frequency response curves (FRC), both frequency-amplitude and frequency-phase. Also includes the stability information if given.

Parameters:

FRC (dict) – Dictionary containing the frequency response curves and the bifurcation curves.

Returns:

• fig1 (Figure) – The amplitude plot \(a(\omega)\).

• fig2 (Figure) – The phase plot \(\beta(\omega)\).

sympy_functions Module

Started on Wed Apr 9 13:39:41 2025

@author: Vincent MAHE

Sympy functions useful to the use of MMS functions.

oscilate.sympy_functions.check_solvability(poly, x)

Check the solvability of a polynomial \(p(x)\).

Parameters:

• poly (sympy.Expr) – The polynomial considered.

• x (sympy.Symbol) – The variable to solve for.

Returns:

bool – True is solvable, False otherwise.

Return type:

bool,

oscilate.sympy_functions.get_block_diagonal_indices(matrix, block_sizes)

Generate a list of \((i, j)\) indices for all elements in the diagonal blocks of a block-diagonal matrix.

Parameters:

• matrix (sympy.Matrix) – The matrix to check for block diagonality.

• block_sizes (int or list of int) – Size(s) of the diagonal blocks.

Returns:

indices – A list of tuples (i, j) representing the indices of elements in the diagonal blocks.

Return type:

list

oscilate.sympy_functions.get_exponent(expr, x)

Get the exponent of \(x\) in an expression of the type \(\lambda x^n\) where \(\lambda\) is a constant while \(n\) is an integer or rational.

Parameters:

• expr (sympy.Expr) – The expression in which one wants to identify the exponent of x.

• x (sympy.Symbol) – The variable whose exponent is to be known.

oscilate.sympy_functions.is_block_diagonal(matrix, block_sizes)

Check if a matrix is block-diagonal given block sizes.

Parameters:

• matrix (sympy.Matrix) – The matrix to check for block diagonality.

• block_sizes (int or list of int) – Size(s) of the diagonal blocks.

Returns:

bool

Return type:

True if the matrix is block-diagonal, False otherwise.

oscilate.sympy_functions.polynomial_terms(poly, x)

Identify the terms of a polynomial.

Parameters:

• poly (sympy.Expr) – The polynomial considered.

• x (sympy.Symbol) – The variable to solve for.

Returns:

dic_x – The polynomial terms.

Return type:

dict

Notes

If the expression given for poly is of the form

\[p(x) = q(x) x^{-n},\]

where the powers of \(x\) in \(q(x)\) are all superior or equal to \(0\), then an auxiliary polynomial

\[P(x) = \dfrac{p(x)}{x^{-n}}\]

is constructed. It is the terms of that positive powers polynomial that are returned.

oscilate.sympy_functions.solve_poly2(poly, x)

Finds the roots of a polynomial of degree 2.

Parameters:

• poly (sympy.Expr) – polynomial whose roots are to be computed

• x (sympy.Symbol) – Variable of the polynomial

Returns:

x_sol – list containing the two roots of the polynomial

Return type:

list of sympy.Expr

Notes

The polynomial of degree 2 takes the form

\[p(x) = a x^2 + bx + c.\]

Note that \(b\) can be null but not \(a\) nor \(c\).

It is a workaround to using solve() or solveset(). These two work but can be very long when coefficients \(a,\; b,\; c\) are expressions involving many parameters. Note that solve() is significantly slower than solveset().

oscilate.sympy_functions.sub_deep(expr, sub)

Performs deep substitutions of an expression.

Parameters:

• expr (sympy.Expr) – Expression on which substitutions are to be performed.

• sub (list of tuples) – The substitutions to perform.

Returns:

expr_sub – The expression with substitutions performed.

Return type:

sympy.Expr

Notes

Deep substitutions are needed when a substitution involves terms that can still be substituted. For instance, one wants to substitute \(a_1\) and \(a_2\) by expressions, but \(a_1\) is actually a function of \(a_2\), so at least 2 substitutions are required.

oscilate.sympy_functions.sympy_to_numpy(expr_sy, param)

Transform a sympy expression into a numpy array.

Parameters:

• expr_sy (sympy.Expr) – A sympy expression.

• param (dict) –

  A dictionnary whose values are tuples with 2 elements:

  1. The sympy symbol of a parameter

  2. The numerical value(s) taken by that parameter

Returns:

expr_np – The numerical values taken by the sympy expression evaluated.

Return type:

numpy.ndarray