Chapter 6
Forced Oscillators of Systems with Finite Degrees of Freedom

In this chapter, we discuss nonlinear systems having finite degrees of freedom. The discussion is limited to weakly nonlinear systems, and approximate solutions are obtained by using the method of multiple scales. In the case of strongly nonlinear systems, perturbation methods can be used in cases for which a basic exact nonlinear solution exists. For the other cases, recourse is often made to numerical methods and / or geometrical methods to obtain a qualitative description of the behavior of the system, including its stability.

In contrast with a single-degree-of-freedom system, which has only a single linear natural frequency and a single mode of oscillation, an n-degree-of-freedom system has n linear natural frequencies and n corresponding modes. Let us denote these frequencies by (Subscript[ω, 1]), (Subscript[ω, 2]), …, 
      (Subscript[ω, n]) and assume that all of them are real and different from zero. An important case occurs whenever two or more of these frequencies are commensurate or nearly commensurate. Examples of near-commensurability are

{(Subscript[ω, 2]) ≈ 2  (Subscript[ω, 1]), 
      (Subscript[ω, 3]) ≈ (Subscript[ω, 2]) ± (Subscript[ω, 1]), 
      (Subscript[ω, 2]) ≈ 3  (Subscript[ω, 1]), 
      (Subscript[ω, 3]) ≈ 2  (Subscript[ω, 2]) ± (Subscript[ω, 1]), 
      (Subscript[ω, 4]) ≈ 
        (Subscript[ω, 3]) ± (Subscript[ω, 2]) ± (Subscript[ω, 1])}

Depending on the order of the nonlinearity in the system, these commensurate relationships of frequencies can cause the corresponding modes to be strongly coupled, and an internal or autoparametric resonance is said to exist. For example, if the system has quadratic nonlinearities, then to first order an internal resonance may exist if (Subscript[ω, m]) ≈ 2  (Subscript[ω, k]) or \(TraditionalForm
      \`ω\_q ≈ ω\_p ± ω\_m\). For a system with cubic nonlinearities, to first order an internal resonance may exist if (Subscript[ω, m]) ≈ 3  (Subscript[ω, k]) or \(TraditionalForm
      \`ω\_q ≈ 
        \(ω\_p ± ω\_m\) ± 
          ω\_k\) or (Subscript[ω, m]) ≈ (Subscript[ω, k]) or \(TraditionalForm
      \`ω\_q ≈ 
        2  ω\_p ± ω\_m\). When an internal resonance exists in a free undamped system, energy imparted initially to one of the modes involved in the internal resonance will be continuously exchanged among the modes involved in that internal resonance. If damping is present in the system, then the energy will decay with time as it is being continuously exchanged.

In a conservative nongyroscopic system, if the linear motion is oscillatory, then the nonlinear motion is bounded. For a conservative gyroscopic multidegree-of-freedom system, the nonlinear motion may be unbounded and hence unstable if an internal resonance exists.

If an external harmonic excitation of frequency Ω acts on a multidegree-of-freedom system, then in addition to all of the primary and secondary resonances (p Ω ≈ q (Subscript[ω, m]), with p and q being integers) of a single-degree-of-freedom system, there might exist other resonance combinations of the frequencies in the form

p Ω == 
      (Subscript[a, 1]) (Subscript[ω, 1]) + (Subscript[a, 2]) (Subscript[ω, 2]) + ⋯ + 
        (Subscript[a, N]) (Subscript[ω, N])

where p and the (Subscript[a, n]) are positive or negative integers such that

| p | +(∑_(n = 1)^N) | (Subscript[a, n]) |  == M

where M is the order of the nonlinearity plus one and N is the number of degrees of freedom. The type of combination resonance which might exist in a system depends on the order of its nonlinearity. For a system having quadratic nonlinearities, to first order the combination resonances that might exist involve two frequencies in addition to Ω;  that is, \(TraditionalForm
      \`Ω ≈ ω\_m + ω\_k\) or \(TraditionalForm
      \`Ω ≈ ω\_m - ω\_k\). The first of these is called a summed combination resonance or a combination resonance of the additive type, whereas the second is called a difference combination resonance or a combination resonance of the difference type. These types of combination resonances were predicted theoretically by Malkin (1956) and found experimentally by Yamamoto (1957, 1960). For a system having cubic nonlinearities, to first order the combination resonances that might exist involve either two or three of the natural frequencies in addition to Ω; that is,

{Ω ≈ 
        (Subscript[ω, p]) ± (Subscript[ω, m]) ± (Subscript[ω, k]), 
      Ω ≈ 2  (Subscript[ω, p]) ± (Subscript[ω, m]), 
      Ω ≈ (Subscript[ω, p]) ± 2  (Subscript[ω, m]), 
      2  Ω ≈ (Subscript[ω, m]) ± (Subscript[ω, k])}

Preliminaries

Off[General::spell1]

Needs["Utilities`Notation`"]

To use the method of multiple scales, we introduce different time scales, symbolize them as

Symbolize[(Subscript[T, 0])];Symbolize[(Subscript[T, 1])];Symbolize[(Subscript[T, 2])];

and form a list of them as follows:

timeScales = {(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])};

In terms of the time scales (Subscript[T, n]), the time derivatives become

\(\(\ \(dt[1]\)[expr_] := 
      Sum[ϵ\^i\ D[expr, timeScales[\([i + 1]\)]], {i, 0, maxOrder}]; 
    
\ \(dt[2]\)[expr_] := 
      \((\(dt[1]\)[\(dt[1]\)[expr]] // Expand)\) /. 
        ϵ\^\(i_ /; i > maxOrder\) -> 0; \)\)

In the course of the analysis, we need the complex conjugates of A and Γ. We define them using the following rule:

\(\(conjugateRule = {A -> A\&_, A\&_ -> A, 
        Γ -> Γ\&_, 
        Γ\&_ -> Γ, 
        Complex[0, n_] -> Complex[0, \(-n\)]}; \)\)

To represent some of the expressions in a more concise way, we introduce the following display rule:

\(\(displayRule = {
        \(\(Derivative[a__]\)[u_\_\(i_, j_\)]\)[__] :> 
          SequenceForm[
            Times\ @@\ MapIndexed[D\_\(#2[\([1]\)] - 1\)\%#1\ &, {a}], 
            u\_\(i, j\)], 
        \(\(Derivative[a__]\)[A_\_i_]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\_i], u_\_\(i_, j_\)[__] :> u\_\(i, j\), A_\_i_[__] :> A\_i}; 
    \)\)

6.1 Externally Excited Linearly Uncoupled Systems

System of Equations

We consider the response of the following two-degree-of-freedom system with quadratic and cubic geometric nonlinearities to an external (additive) excitation:

\(eq61a = {

        \(u\_1\^″\)[t] + ω\_1\%2\ u\_1[t] + 
            2  μ\&^\_1\ \(u\_1\^′\)[t] + 
            D[V[u\_1[t], u\_2[t]], u\_1[t]] == 
          F\_1\ Cos[Ω\ t + τ\_1], 

        \(u\_2\^″\)[t] + ω\_2\%2\ u\_2[t] + 
            2  μ\&^\_2\ \(u\_2\^′\)[t] + 
            D[V[u\_1[t], u\_2[t]], u\_2[t]] == 
          F\_2\ Cos[Ω\ t + τ\_2]}\)

where the potential function V possesses cubic and quartic nonlinearities. The most general such potential function can be expressed as follows:

basicModes = {(Subscript[u, 1])[t], (Subscript[u, 2])[t]};

cubicTerms = 
      Nest[Outer[Times, basicModes, #]&, basicModes, 2] // Flatten // 
        Union

quarticTerms = 
      Outer[Times, basicModes, cubicTerms] // Flatten // Union

potential = 
      V -> (Evaluate[
              Sum[(Subscript[δ, i]) cubicTerms[[i]], {i, 
                      Length[cubicTerms]}] + 
                  Sum[(Subscript[α, i]) quarticTerms[[i]], {i, 
                      Length[quarticTerms]}] /. 
                Thread[basicModes -> {#1, #2}]] &)

Approximate solutions of such a system can be obtained using a number of techniques, including the method of normal forms, the method of multiple scales, the Krylov-Bogoliubov-Mitropolsky technique, and the method of averaging. With these methods, one determines a set of ordinary-differential equations (modulation equations) that govern the time evolution of the amplitudes and phases of the modes participating in the response. If these equations are time-dependent, one usually transforms them into a set of autonomous equations. In this chapter, we describe how the method of multiple scales can be used to determine uniform expansions of this system, including the modulation equations.

Secondary resonances of nonlinear systems include subharmonic, superharmonic, and combination resonances. One or more of these resonances might be activated in the presence or absence of internal resonances. When a combination resonance is activated in the presence of an internal resonance, fractional harmonics might be excited.

For a consistent expansion, we first transform eq61a into a system of four first-order equations. To this end, we introduce the two states (Subscript[v, 1])[t] and (Subscript[v, 2])[t] defined by

vel = {((Subscript[u, 1])^′)[t] -> (Subscript[v, 1])[t], 
        ((Subscript[u, 2])^′)[t] -> (Subscript[v, 2])[t]};

Substituting the velocity and acceleration terms, using vel and potential, into eq61a and combining the result with vel, we transform eq61a into the following set of four first-order equations:

eq61b = 
      {vel /. Rule -> Equal, 
            eq61a /. D[vel, t] /. vel /. potential} // Transpose // 
        Flatten

We seek a uniform second-order expansion of the solution of eq61b in the form

solRule = {
        (Subscript[u, i])_ -> 
          (Sum[(ϵ^j) (Subscript[u, i, j])[#1, #2, #3], {j, 3}] &), 
        (Subscript[v, i])_ -> 
          (Sum[(ϵ^j) (Subscript[v, i, j])[#1, #2, #3], {j, 3}] &)};

and set

maxOrder = 2;

The damping and forcing terms in eq61b need to be scaled so that they balance the influence of the nonlinearity. The scaling of the forcing depends on the type of resonance. Next, we treat different resonance cases in the following sections.

6.1.1 Ω ≈ (Subscript[ω, 2]) and (Subscript[ω, 2]) ≈ 2  (Subscript[ω, 1])

In this case, we have a combination of a primary resonance of the second mode, a subharmonic resonance of order one-half of the first mode, and a two-to-one internal resonance. In order to bring the effects of damping, forcing, and nonlinearity at the same order, we let

\(\(scaling = {μ\&^\_n_ -> ϵ\ μ\_n, 
        F\_1 -> ϵ\ f\_1, F\_2 -> ϵ\^2\ f\_2}; \)\)

To describe quantitatively the nearness of the resonances, we introduce the two detunings (Subscript[σ, 1]) and (Subscript[σ, 2]) defined by

ResonanceCond = {
        (Subscript[ω, 2]) == 2  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
        Ω == (Subscript[ω, 2]) + ϵ (Subscript[σ, 2])};

Then, we define the following lists:

omgList = {(Subscript[ω, 1]), (Subscript[ω, 2])};

omgRule = 
      Solve[ResonanceCond, {Ω, #} // Flatten][[1]] &
             /@ omgList // Reverse

Using the time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), we express the time derivative as

multiScales = {(Subscript[u, i])_[t] -> (Subscript[u, i]) @@ timeScales, 
        (Subscript[v, i])_[t] -> (Subscript[v, i]) @@ timeScales, 
        Derivative[1][u_][t] -> dt[1][u @@ timeScales], 
        t -> (Subscript[T, 0])};

Substituting the scaling, multiScales, and solRule into eq61b, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq611a = 
      (eq61b /. scaling /. multiScales /. solRule // TrigToExp
               // ExpandAll) /. (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we obtain

eqEps = 
      Rest[Thread[CoefficientList[Subtract @@ #, ϵ] == 0]] &
             /@ eq611a // Transpose;

To place the linear operator on one side and the nonhomogeneous terms on the other side, we define

\(eqOrder[i_] := 
      \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f\_s_ -> 0\) /. 
              u_\_\(k_, 1\) -> u\_\(k, i\))\) == 
          \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f\_s_ -> 0\) /. 
                u_\_\(k_, 1\) -> u\_\(k, i\))\) - 
            \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([i]\)])\) // Thread\)

Using the eqOrder[i] and displayRule, we rewrite eqEps in a concise way as

eqOrder[1] /. displayRuleeqOrder[2] /. displayRuleeqOrder[3] /. displayRule

First-Order Problem: Linear System

The first-order problem, eqOrder[1], consists of two sets of uncoupled linear nonhomogeneous differential equations. Hence, the general solution of each set can be obtained by using the principle of superposition as the sum of any particular solution and the solution of the homogeneous equations. To determine the general solution of the homogeneous sets, we rewrite the homogeneous part of eqOrder[1] as

\(linearSys = \(#[\([1]\)]\ &\)\ /@\ eqOrder[1]; 

    linearSys /. displayRule\)

Next, we seek a solution of the linearSys in the form (Subscript[u, i, 1]) = P (E^(I (Subscript[ω, i]) (Subscript[T, 0]))) and (Subscript[v, i, 1]) = Q (E^(I (Subscript[ω, i]) (Subscript[T, 0]))) and obtain

expr1 = 
      linearSys /. {(Subscript[u, i_, 1]) -> (P (E^(I (Subscript[ω, i]) #1)) &), 
              (Subscript[v, i_, 1]) -> (Q (E^(I (Subscript[ω, i]) #1)) &)} /. 
          Exp[a_] -> 1 // Partition[#, 2] &

The coefficient matrices of expr1 are

coefMat = Outer[D, #, {P, Q}] & /@ expr1

and their adjoints are defined by

hermitian[mat_?MatrixQ] := mat /. conjugateRule // Transpose

Hence, the right and left eigenvectors of coefMat are

rightVec = 
      #/#[[1]] & /@ 
        (NullSpace[#][[1]] & /@ coefMat)

leftVec = NullSpace[hermitian[#]][[1]] & /@ coefMat

whose complex conjugates are

ccleftVec = leftVec /. conjugateRule

To determine particular solutions of the partial-differential equations eqOrder[1] using the Mathematica function DSolve, we first transform them into a set of ordinary-differential equations as

\(order1Eq = eqOrder[1] /. u_\_\(i_, 1\) -> \((u\_\(i, 1\)[#1]\ &)\)\)

Then, the particular solutions of order1Eq are

sol1p = 
      DSolve[order1Eq, {(Subscript[u, 1, 1])[(Subscript[T, 0])], (Subscript[v, 1, 1])[(Subscript[T, 0])], 
                (Subscript[u, 2, 1])[(Subscript[T, 0])], (Subscript[v, 2, 1])[(Subscript[T, 0])]}, (Subscript[T, 0])][[1]] /. 
          C[_] -> 0 // Simplify

To simplify some of the expressions, we let

fRule = {
        (Subscript[f, i])_ -> 
          2 (Subscript[Λ, i]) ((Subscript[ω, i]^2) - (Ω^2))};

and express (Subscript[u, i, 1]) and (Subscript[v, i, 1]) in pure function form as

\(sol1u = 
      Table[u\_\(i, 1\) -> 
          Function[{T\_0, T\_1, T\_2}, 
            A\_i[T\_1, T\_2]\ Exp[I\ ω\_i\ T\_0] + 
                A\&_\_i[T\_1, T\_2]\ Exp[\(-I\)\ ω\_i\ T\_0] + 
                \((\(\(u\_\(i, 1\)[T\_0] /. sol1p\) /. fRule // Simplify\) // 
                    Expand)\) // Evaluate], {i, 2}]\)

sol1v = 
      Table[(Subscript[v, i, 1]) -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            D[(Subscript[u, i, 1]) @@ timeScales, (Subscript[T, 0])] /. sol1u // Evaluate], {i, 
          2}]

Then, we combine them using the Mathematica function Join as

sol1 = Join[sol1u, sol1v];

Second-Order Problem

Substituting the first-order solution sol1 into the second-order problem, eqOrder[2], yields

\(order2Eq = eqOrder[2] /. sol1 // ExpandAll; 

    order2Eq /. displayRule\)

Any particular solution of order2Eq contains secular terms and small-divisor terms, making the expansion nonuniform. For a uniform expansion, we choose (Subscript[D, 1]) (Subscript[A, 1]) and (Subscript[D, 1]) (Subscript[A, 2]) to eliminate these terms. To accomplish this, we first convert the small-divisor terms into secular terms using the rule

expRule1[i_] := 
      Exp[a_] :> 
        Exp[Expand[a /. omgRule[[i]]] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]

To eliminate the terms that produce secular terms (i.e., determine the solvability conditions) from order2Eq, we collect the terms proportional to (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))) and (E^(I (Subscript[ω, 2]) (Subscript[T, 0]))) and obtain

\(ST11 = 
      \(Coefficient[order2Eq[\([#, 2]\)] /. expRule1[1], 
            Exp[I\ ω\_1\ T\_0]]\ &\)\ /@\ {1, 2}; 

    ST11 /. displayRule\)

\(ST12 = 
      \(Coefficient[order2Eq[\([#, 2]\)] /. expRule1[2], 
            Exp[I\ ω\_2\ T\_0]]\ &\)\ /@\ {3, 4}; 

    ST12 /. displayRule\)

Then, the solvability conditions, conditions for the elimination of the terms that produce secular terms, demand that ST11 and ST12 be orthogonal to every solution of the corresponding adjoint problem, namely, the components of ccleftVec. The result is

\(SCond1 = {ccleftVec[\([1]\)] . ST11 == 0, 
        ccleftVec[\([2]\)] . ST12 == 0}; 
SCond1 /. displayRule\)

Solving SCond1 for (Subscript[D, 1]) (Subscript[A, 1]) and (Subscript[D, 1]) (Subscript[A, 2]) yields

SCond1Rule1=Solve[SCond1,{(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]//ExpandAll;
SCond1Rule1 /. displayRule

With these conditions, order2Eq are solvable. However, their solutions are not unique. To render them unique, we demand that they be orthogonal to solutions of the corresponding adjoint problems; that is, we damand that {(Subscript[u, 1, 2 ]), (Subscript[v, 1, 2])} be orthogonal to ccleftVec[[1]] and {(Subscript[u, 2, 2 ]), (Subscript[v, 2, 2])} be orthogonal to ccleftVec[[2]].

Next, we use the solvability conditions to eliminate (Subscript[D, 1]) (Subscript[A, 1]) and (Subscript[D, 1]) (Subscript[A, 2]) from order2Eq. To simplify the resulting equations using Mathematica, we replace the scale (Subscript[T, 1]) in SCond1Rule1 with ϵ (Subscript[T, 0]). To this end, we express the detuning parameters (Subscript[σ, i]) in terms of the (Subscript[ω, i]) and Ω as

sigRule = 
      Solve[ResonanceCond, {(Subscript[σ, 1]), (Subscript[σ, 2])}][[1]]

and replace (Subscript[T, 1]) with ϵ (Subscript[T, 0]) using the rule

expRule2 = 
      Exp[a_] :> Exp[a /. (Subscript[T, 1]) -> ϵ (Subscript[T, 0]) /. sigRule // Expand];

Using the expRule2, we rewrite SCond1Rule1 as

\(SCond1Rule2 = SCond1Rule1 /. expRule2; 
SCond1Rule2 /. displayRule\)

whose complex conjugate is

ccSCond1Rule2 = SCond1Rule2 /. conjugateRule;

Substituting this form for the solvability conditions into order2Eq yields

\(order2Eqm = \(order2Eq /. SCond1Rule2\) /. ccSCond1Rule2 // ExpandAll; 
    
order2Eqm /. displayRule\)

Next, we use the method of undetermined coefficients to determine the particular solutions of order2Eqm. To accomplish this, we first identify the form of the nonhomogeneous terms. To this end, we let

\(basicH = 
      Table[{A\_i[T\_1, T\_2]\ E\^\(I\ ω\_i\ T\_0\), 
            A\&_\_i[T\_1, T\_2]\ E\^\(\(-I\)\ ω\_i\ T\_0\)}, {i, 
            2}] // Flatten\)

collectForm = 
      Join[basicH, {
          (Subscript[Λ, 1]) 
            (E^(I Ω (Subscript[T, 0]) + I (Subscript[τ, 1]))), 
          (Subscript[Λ, 1]) 
            (E^(-I Ω (Subscript[T, 0]) - I (Subscript[τ, 1])))}]

Then, the possible form of the nonhomogeneous terms in order2Eqm is

possibleTerms = 
      Join[collectForm, {(E^(I Ω (Subscript[T, 0]) + I (Subscript[τ, 2]))), 
          (E^(-I Ω (Subscript[T, 0]) - I (Subscript[τ, 2])))}, 
        Outer[Times, collectForm, collectForm] // Flatten // Union]

Using the method of undetermined coefficients, we assume that the (Subscript[u, i, 2]) and (Subscript[v, i, 2]) are linear combinations of these possible terms, substitute the result into the governing equations, order2Eqm, equate the coefficients of each of the possible terms on both sides of each equation, and obtain a system of nonhomogeneous algebraic equations governing the unknown coefficients. Associated with each possible term is a pair of two algebraic equations. The solutions of all pairs, except those associated with resonance terms, are unique. The solutions of the pairs corresponding to the possible resonance terms are rendered unique by requiring them to be orthogonal to the components of ccleftVec. The possible resonance terms are proportional to (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 2]) (Subscript[T, 0]))),  and their complex conjugates; that is,

ResonantTerms[i_] := 
      (# /. {a_ /; a =!= 0 -> 1} & /@ 
                ((E^(-I (Subscript[ω, i]) (Subscript[T, 0]))) possibleTerms /. 
                      expRule1[i] /. Exp[_ (Subscript[T, 0]) + _. ] -> 0)) 
            possibleTerms // Union // Rest

RT = Array[ResonantTerms, 2]

Hence, the coefficients of RT in the nonhomogeneous terms in order2Eqm[[1]] are

r1Rule = 
      MapIndexed[
          (Subscript[r, 1, #2[[1]]]) -> Coefficient[order2Eqm[[1, 2]], #1] &, 
          RT[[1]]] /. Exp[_ (Subscript[T, 0]) + _. ] -> 0

and the coefficients of RT in the nonhomogeneous terms in order2Eqm[[2]] are

r2Rule = 
      MapIndexed[
          (Subscript[r, 2, #2[[1]]]) -> Coefficient[order2Eqm[[2, 2]], #1] &, 
          RT[[1]]] /. Exp[_ (Subscript[T, 0]) + _. ] -> 0

It follows from

Table[ccleftVec[[1]] . {(Subscript[r, 1, i]), (Subscript[r, 2, i])}, {i, 3}] /. 
        r1Rule /. r2Rule

that the parts r1Rule and r2Rule are orthogonal to the solution of the adjoint, the first component of ccleftVec, as they should.

Similarly, the coefficients of RT in the nonhomogeneous terms in order2Eqm[[3]] are

r3Rule = 
      MapIndexed[
        (Subscript[r, 3, #2[[1]]]) -> Coefficient[order2Eqm[[3, 2]], #1] &, 
        RT[[2]]]

and the coefficients of RT in the nonhomogeneous terms in order2Eqm[[4]] are

r4Rule = 
      MapIndexed[
          (Subscript[r, 4, #2[[1]]]) -> Coefficient[order2Eqm[[4, 2]], #1] &, 
          RT[[2]]] /. Exp[_ (Subscript[T, 0]) + _. ] -> 0

Again, it follows from

Table[ccleftVec[[2]] . {(Subscript[r, 3, i]), (Subscript[r, 4, i])}, {i, 4}] /. 
        r3Rule /. r4Rule

that the vector {(Subscript[r, 3,  i ]), (Subscript[r, 4,  i])} is orthogonal to the solution of the adjoint, the second component of ccleftVec, as it should.

The complement of the resonance terms RT in possibleTerms yields the nonresonance terms; that is,

NRT = 
      Complement[possibleTerms, Join[#, # /. conjugateRule]] & /@ RT

To implement the method of undetermined coefficients, we associate with each possible resonance term an unknown coefficient using the rules

RTsymbolList1 = 
      Table[(Subscript[Γ, 1, j]), {j, Length[RT[[1]]]}]

RTsymbolList2 = 
      Table[(Subscript[Γ, 3, j]), {j, Length[RT[[2]]]}]

And, we associate with each possible nonresonance term two unknown coefficients using the rules

NRTsymbolList1[i_] = 
      Table[(Subscript[Λ, i, j]), {j, Length[NRT[[1]]]}]

NRTsymbolList2[i_] = 
      Table[(Subscript[Λ, i + 2, j]), {j, Length[NRT[[2]]]}]

Accounting for the orthogonality conditions, {(Subscript[u, 1, 2 ]), (Subscript[v, 1, 2])} must be orthogonal to ccleftVec[[1]] and {(Subscript[u, 2, 2 ]), (Subscript[v, 2, 2])} must be orthogonal to ccleftVec[[2]], for the resonance terms, RT, we express (Subscript[u, i, 2]) and (Subscript[v, i, 2]) in the following function form:

{ratio1, ratio2} = Table[-ccleftVec[[i, 1]], {i, 2}]

sol2Form = {
        (Subscript[u, 1, 2]) -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            RTsymbolList1 . RT[[1]] + 
                (RTsymbolList1 . RT[[1]] /. conjugateRule) + 
                NRTsymbolList1[1] . NRT[[1]] // Evaluate], 
        (Subscript[v, 1, 2]) -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            ratio1*RTsymbolList1 . RT[[1]] + 
                (ratio1*RTsymbolList1 . RT[[1]] /. conjugateRule) + 
                NRTsymbolList1[2] . NRT[[1]] // Evaluate], 
        (Subscript[u, 2, 2]) -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            RTsymbolList2 . RT[[2]] + 
                (RTsymbolList2 . RT[[2]] /. conjugateRule) + 
                NRTsymbolList2[1] . NRT[[2]] // Evaluate], 
        (Subscript[v, 2, 2]) -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            ratio2*RTsymbolList2 . RT[[2]] + 
                (ratio2*RTsymbolList2 . RT[[2]] /. conjugateRule) + 
                NRTsymbolList2[2] . NRT[[2]] // Evaluate]};

Substituting sol2Form into order2Eqm and collecting the coefficients of like NRT, we obtain two sets of algebraic equations as

eq611b = 
      (Coefficient[Subtract @@ # /. sol2Form, NRT[[1]]] /. 
                    Exp[_ (Subscript[T, 0]) + _. ] -> 0 & /@ 
                order2Eqm[[{1, 2}]] // Flatten) == 0 // Thread

eq611c = 
      (Coefficient[Subtract @@ # /. sol2Form, NRT[[2]]] /. 
                    Exp[_ (Subscript[T, 0]) + _. ] -> 0 & /@ 
                order2Eqm[[{3, 4}]] // Flatten) == 0 // Thread

Solving eq611b for the NRTsymbolList1, we obtain

coef11 = 
      Solve[eq611b, Array[NRTsymbolList1, 2] // Flatten][[1]]

Solving eq611c for the NRTsymbolList2, we obtain

coef12 = 
      Solve[eq611c, Array[NRTsymbolList2, 2] // Flatten][[1]]

Substituting sol2Form into order2Eqm and collecting the coefficients of like RT, we obtain two sets of algebraic equations as

eq1 = 
      MapIndexed[
          Coefficient[order2Eqm[[1, 1]] /. sol2Form, #1] == 
              (Subscript[r, 1, #2[[1]]]) &, RT[[1]]] /. 
        Exp[_ (Subscript[T, 0]) + _. ] -> 0

eq2 = 
      MapIndexed[
          Coefficient[order2Eqm[[3, 1]] /. sol2Form, #1] == 
              (Subscript[r, 3, #2[[1]]]) &, RT[[2]]] /. 
        Exp[_ (Subscript[T, 0]) + _. ] -> 0

Solving eq1 and eq2 for the RTsymbolList1 and RTsymbolList2, we obtain

coef21 = Solve[eq1, RTsymbolList1][[1]] /. r1Rule

coef22 = Solve[eq2, RTsymbolList2][[1]] /. r3Rule

Combining these coefficients using the function Join yields

coeffs = 
      Join[coef11, coef12, coef21, coef22, coef21 /. conjugateRule, 
          coef22 /. conjugateRule] /. 
        (omgRule[[1]] /. ϵ -> 0);

Then, we express the solution of the second-order equations in pure function form as

\(sol2 = 
      sol2Form /. 
        Function[{T\_0, T\_1, T\_2}, a_] :> 
          Function[{T\_0, T\_1, T\_2} // Evaluate, 
            \(a /. coeffs // Expand\) // Evaluate]; 
sol2 /. displayRule\)

Third-Order Problem

Substituting the first- and second-order solutions into the third-order problem, eqOrder[3], yields

order3Eq = eqOrder[3] /. sol1 /. sol2 // ExpandAll;

Substituting the expRule1 into the right-hand sides of order3Eq and collecting the terms that could produce secular terms, the terms proportional to (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))) and (E^(I (Subscript[ω, 2]) (Subscript[T, 0]))), we have

\(ST21 = 
      \(Coefficient[order3Eq[\([#, 2]\)] /. expRule1[1], 
            Exp[I\ ω\_1\ T\_0]]\ &\)\ /@\ {1, 2}; 

    ST21 /. displayRule\)

\(ST22 = 
      \(Coefficient[order3Eq[\([#, 2]\)] /. expRule1[2], 
            Exp[I\ ω\_2\ T\_0]]\ &\)\ /@\ {3, 4}; 

    ST22 /. displayRule\)

Then, the solvability conditions demand that ST21 and ST22 be orthogonal to solutions of their corresponding adjoints, the components of ccleftVec. The result is

\(SCond2 = {ccleftVec[\([1]\)] . ST21 == 0, 
          ccleftVec[\([2]\)] . ST22 == 0} // ExpandAll; 

    SCond2 /. displayRule\)

Solving SCond2 for (Subscript[D, 2]) (Subscript[A, 1]) and (Subscript[D, 2]) (Subscript[A, 2]), we have

SCond2Rule1=Solve[SCond2,{(Subscript[A, 1]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]//ExpandAll;
SCond2Rule1 /. displayRule

Reconstitution

Using the method of reconstitution, \(TraditionalForm
      \`A\_k\^′ = 
        ϵ\ \(D\_1\) A\_k + ϵ\^2\ \(D\_2\) A\_k + \ 
          ⋯\), we combine the partial-differential equations SCond1Rule1 and SCond2Rule1 into the following two ordinary-differential equations governing the modulation of the complex-valued functions (Subscript[A, k]) :

\(moduEq = 
      Table[2\ I\ ω\_k\ A\_k\^′, {k, 2}] == 
          \((\(Table[
                    2\ I\ ω\_k\ \(dt[1]\)[\ A\_k[T\_1, T\_2]], {k, 
                      2}] /. SCond1Rule1\) /. SCond2Rule1 // 
              Collect[#, ϵ]\ &)\) // Thread; 

    moduEq /. displayRule\)

6.1.2 The function MMS

According to the procedures described in the previous section, we can build a function named MMS (Method of Multiple Scales) specifically for eq61b. A more general function (a Package) can be created by considering as arguments the governing equations, symbols for dependent variables, independent variable, excitation amplitudes and frequencies, and all other related quantities which allow the program to identify their respective meanings. We will use MMS to solve for different resonance cases in the following sections.

MMS[scaling_List, ResonanceCond : {__Equal}]:=
   Module[{},
		omgList = {(Subscript[ω, 1]), (Subscript[ω, 2])};
     omgRule = 
              Solve[ResonanceCond, {Ω, #} // Flatten][
                      [1]] & /@ omgList // Reverse;
     multiScales = {(Subscript[u, i])_[t] -> (Subscript[u, i]) @@ timeScales, 
                (Subscript[v, i])_[t] -> (Subscript[v, i]) @@ timeScales, 
                Derivative[1][u_][t] -> dt[1][u @@ timeScales], 
                t -> (Subscript[T, 0])};
		eqa = (
                  eq61b /. scaling /. multiScales /. solRule // 
                      TrigToExp // ExpandAll) /. 
                (ϵ^(n_ /; n > 3)) -> 0;
		eqEps = 
              Rest[Thread[
                        CoefficientList[Subtract @@ #, ϵ] == 0]] &
                     /@ eqa // Transpose;
		\(eqOrder[i_] := 
              \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f\_s_ -> 0\) /. 
                      u_\_\(k_, 1\) -> u\_\(k, i\))\) == 
                  \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. 
                          f\_s_ -> 0\) /. u_\_\(k_, 1\) -> u\_\(k, i\))\) - 
                    \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([i]\)])\) // Thread\);
	
		  		
		linearSys = #[[1]] & /@ eqOrder[1];
		expr1 = 
              linearSys /. {
                      (Subscript[u, i_, 1]) -> (P (E^(I (Subscript[ω, i]) #1)) &), 
                      (Subscript[v, i_, 1]) -> (Q (E^(I (Subscript[ω, i]) #1)) &)}
                       /. Exp[a_] -> 1 // Partition[#, 2] &;
		coefMat = Outer[D, #, {P, Q}] & /@ expr1;
		hermitian[mat_?MatrixQ] := mat /. conjugateRule // Transpose;
		rightVec = 
              #/#[[1]] & /@ 
                (NullSpace[#][[1]] & /@ coefMat);
		leftVec = 
              NullSpace[hermitian[#]][[1]] & /@ coefMat;
		ccleftVec = leftVec /. conjugateRule;
		\(order1Eq = 
              eqOrder[1] /. u_\_\(i_, 1\) -> \((u\_\(i, 1\)[#1]\ &)\)\);
		sol1p = 
              DSolve[
                      order1Eq, {(Subscript[u, 1, 1])[(Subscript[T, 0])], (Subscript[v, 1, 1])[(Subscript[T, 0])], 
                        (Subscript[u, 2, 1])[(Subscript[T, 0])], (Subscript[v, 2, 1])[(Subscript[T, 0])]}, (Subscript[T, 0])][
                    [1]] /. C[_] -> 0 // Simplify;
		fRule = {
                (Subscript[f, i])_ -> 
                  2 (Subscript[Λ, i]) 
                    ((Subscript[ω, i]^2) - (Ω^2))};
		\(sol1u = 
              Table[u\_\(i, 1\) -> 
                  Function[{T\_0, T\_1, T\_2}, 
                    A\_i[T\_1, T\_2]\ Exp[I\ ω\_i\ T\_0] + 
                        A\&_\_i[T\_1, T\_2]\ Exp[\(-I\)\ ω\_i\ T\_0] + 
                        \((\(\(u\_\(i, 1\)[T\_0] /. sol1p\) /. fRule // 
                              Simplify\) // Expand)\) // Evaluate], {i, 2}]
                \);
		sol1v = 
              Table[(Subscript[v, i, 1]) -> 
                  Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
                    D[(Subscript[u, i, 1]) @@ timeScales, (Subscript[T, 0])] /. sol1u // 
                      Evaluate], {i, 2}];
		sol1 = Join[sol1u, sol1v];
				
		  
		order2Eq = eqOrder[2] /. sol1 // ExpandAll;
		expRule1[i_] := 
              Exp[a_] :> 
                Exp[Expand[a /. omgRule[[i]]] /. 
                    (ϵ^n)_.  (Subscript[T, 0]) :> timeScales[[n + 1]]];
		ST11 = 
              Coefficient[order2Eq[[#, 2]] /. expRule1[1], 
                    Exp[I (Subscript[ω, 1]) (Subscript[T, 0])]] & /@ {1, 2};
		ST12 = 
              Coefficient[order2Eq[[#, 2]] /. expRule1[2], 
                    Exp[I (Subscript[ω, 2]) (Subscript[T, 0])]] & /@ {3, 4};
		SCond1 = {ccleftVec[[1]] . ST11 == 0, 
                ccleftVec[[2]] . ST12 == 0};
		SCond1Rule1=Solve[SCond1,{(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]//ExpandAll;
		sigRule = 
              Solve[ResonanceCond, {(Subscript[σ, 1]), (Subscript[σ, 2])}][[1]]
                ;
		expRule2 = 
              Exp[a_] :> 
                Exp[a /. {(Subscript[T, 1]) -> ϵ (Subscript[T, 0]), 
                          (Subscript[T, 2]) -> (ϵ^2) (Subscript[T, 0])} /. sigRule // Expand]
                  ;
		SCond1Rule2 = SCond1Rule1 /. expRule2;
		ccSCond1Rule2 = SCond1Rule2 /. conjugateRule;
		order2Eqm = 
              order2Eq /. SCond1Rule2 /. ccSCond1Rule2 // ExpandAll;
		\(If[Union[\(#[\([2]\)]\ &\)\ /@\ order2Eqm] === {0}, 
		\ \ \ 
              sol2 = {u\_\(1, 2\) -> \((0\ &)\), v\_\(1, 2\) -> \((0\ &)\), 
                  u\_\(2, 2\) -> \((0\ &)\), v\_\(2, 2\) -> \((0\ &)\)}, 	

              \ \ \ 		
              basicH = 
                Table[{A\_i[T\_1, T\_2]\ E\^\(I\ ω\_i\ T\_0\), 
                      A\&_\_i[T\_1, T\_2]\ 
                        E\^\(\(-I\)\ ω\_i\ T\_0\)}, {i, 2}] // 
                  Flatten; 
	\ \ \ 	
              collectForm = 
                Join[basicH, 
                  If[list1 = 
                      List\ @@\ 
                          \(Plus\ @@\ 
                            \((\(#[\([2]\)]\ &\)\ /@\ eqOrder[1])\)\) /. 
                        c_\ f\_i_\ E\^a_ -> Λ\_i\ E\^a; 
                    Head[list1] === List, list1, {}]]; 
\ \ \ 		
              possibleTerms = 
                Join[collectForm, 
                  If[list1 = 
                      List\ @@\ 
                          \(Plus\ @@\ 
                            \((\(#[\([2]\)]\ &\)\ /@\ eqOrder[2] /. 
                                u_\_\(i_, j_\) -> \((0\ &)\))\)\) /. 
                        c_\ f\_i_\ E\^a_ -> E\^a; Head[list1] === List, 
                    list1, {}], 
                  \(Outer[Times, collectForm, collectForm] // Flatten\) // 
                    Union]; 
\ \ \ 		
              ResonantTerms[i_] := 
                \(\((\(# /. {a_ /; a =!= 0 -> 1}\ &\)\ /@\ 
                          \((\(E\^\(\(-I\)\ ω\_i\ T\_0\)\ 
                                  possibleTerms /. expRule1[i]\) /. 
                              Exp[_\ T\_0 + _. ] -> 0)\))\)\ possibleTerms // 
                    Union\) // Rest; 
	\ \ \ 	
              RT = Array[ResonantTerms, 2]; 
\ \ \ 		
              r1Rule = 
                MapIndexed[
                    r\_\(1, #2[\([1]\)]\) -> 
                        Coefficient[order2Eqm[\([1, 2]\)], #1]\ &, 
                    RT[\([1]\)]] /. Exp[_\ T\_0 + _. ] -> 0; 

              		 (*\ r2Rule = 
                  MapIndexed[
                      r\_\(2, #2[\([1]\)]\) -> 
                          Coefficient[order2Eqm[\([2, 2]\)], #1]\ &, 
                      RT[\([1]\)]] /. Exp[_\ T\_0 + _. ] -> 0; \ *) 

              \ \ \ 		
              r3Rule = 
                MapIndexed[
                    r\_\(3, #2[\([1]\)]\) -> 
                        Coefficient[order2Eqm[\([3, 2]\)], #1]\ &, 
                    RT[\([2]\)]] /. Exp[_\ T\_0 + _. ] -> 0; 

              		 (*\ r4Rule = 
                  MapIndexed[
                      r\_\(4, #2[\([1]\)]\) -> 
                          Coefficient[order2Eqm[\([4, 2]\)], #1]\ &, 
                      RT[\([2]\)]] /. Exp[_\ T\_0 + _. ] -> 0; \ *) 

              \ \ \ 		
              NRT = \(Complement[possibleTerms, Join[#, # /. conjugateRule]]
                      \ &\)\ /@\ RT; 
\ \ \ 		
              RTsymbolList1 = 
                Table[Γ\_\(1, j\), {j, Length[RT[\([1]\)]]}]; 

              \ \ \ 		
              RTsymbolList2 = 
                Table[Γ\_\(3, j\), {j, Length[RT[\([2]\)]]}]; 

              \ \ \ 		
              NRTsymbolList1[i_] = 
                Table[Λ\_\(i, j\), {j, Length[NRT[\([1]\)]]}]; 
              
\ \ \ 		
              NRTsymbolList2[i_] = 
                Table[Λ\_\(i + 2, j\), {j, 
                    Length[NRT[\([2]\)]]}]; 

              \ \ \ 		{ratio1, ratio2} = 
                Table[\(-ccleftVec[\([i, 1]\)]\), {i, 2}]; 
\ \ \ 		
              sol2Form = {
                  u\_\(1, 2\) -> 
                    Function[{T\_0, T\_1, T\_2}, 
                      RTsymbolList1 . RT[\([1]\)] + 
                          \((RTsymbolList1 . RT[\([1]\)] /. conjugateRule)
                            \) + NRTsymbolList1[1] . NRT[\([1]\)] // 
                        Evaluate], 
                  v\_\(1, 2\) -> 
                    Function[{T\_0, T\_1, T\_2}, 
                      ratio1*RTsymbolList1 . RT[\([1]\)] + 
                          \((ratio1*RTsymbolList1 . RT[\([1]\)] /. 
                              conjugateRule)\) + 
                          NRTsymbolList1[2] . NRT[\([1]\)] // Evaluate], 
                  u\_\(2, 2\) -> 
                    Function[{T\_0, T\_1, T\_2}, 
                      RTsymbolList2 . RT[\([2]\)] + 
                          \((RTsymbolList2 . RT[\([2]\)] /. conjugateRule)
                            \) + NRTsymbolList2[1] . NRT[\([2]\)] // 
                        Evaluate], 
                  v\_\(2, 2\) -> 
                    Function[{T\_0, T\_1, T\_2}, 
                      ratio2*RTsymbolList2 . RT[\([2]\)] + 
                          \((ratio2*RTsymbolList2 . RT[\([2]\)] /. 
                              conjugateRule)\) + 
                          NRTsymbolList2[2] . NRT[\([2]\)] // Evaluate]}; 

              \ \ \ 		
              eqb = \((
                      \(Coefficient[Subtract\ @@\ # /. sol2Form, 
                                NRT[\([1]\)]] /. Exp[_\ T\_0 + _. ] -> 0\ &
                            \)\ /@\ order2Eqm[\([{1, 2}]\)] // Flatten)\) == 
                    0 // Thread; 
\ \ \ 		
              eqc = \((
                      \(Coefficient[Subtract\ @@\ # /. sol2Form, 
                                NRT[\([2]\)]] /. Exp[_\ T\_0 + _. ] -> 0\ &
                            \)\ /@\ order2Eqm[\([{3, 4}]\)] // Flatten)\) == 
                    0 // Thread; 
\ \ \ 		
              coef11 = 
                \(Solve[eqb, Array[NRTsymbolList1, 2] // Flatten]\)[\([1]\)]; 
              
\ \ \ 		
              coef12 = 
                \(Solve[eqc, Array[NRTsymbolList2, 2] // Flatten]\)[\([1]\)]; 
              
\ \ \ 		
              eqd = MapIndexed[
                    Coefficient[order2Eqm[\([1, 1]\)] /. sol2Form, #1] == 
                        r\_\(1, #2[\([1]\)]\)\ &, RT[\([1]\)]] /. 
                  Exp[_\ T\_0 + _. ] -> 0; 
\ \ \ 		
              eqe = MapIndexed[
                    Coefficient[order2Eqm[\([3, 1]\)] /. sol2Form, #1] == 
                        r\_\(3, #2[\([1]\)]\)\ &, RT[\([2]\)]] /. 
                  Exp[_\ T\_0 + _. ] -> 0; 
\ \ \ 		
              coef21 = \(Solve[eqd, RTsymbolList1]\)[\([1]\)] /. r1Rule; 

              \ \ \ 		
              coef22 = \(Solve[eqe, RTsymbolList2]\)[\([1]\)] /. r3Rule; 

              \ \ \ 		
              coeffs = 
                Join[coef11, coef12, coef21, coef22, coef21 /. conjugateRule, 
                    coef22 /. conjugateRule] /. 
                  \((omgRule[\([1]\)] /. ϵ -> 0)\); 
\ \ \ 		
              sol2 = sol2Form /. 
                  Function[{T\_0, T\_1, T\_2}, a_] :> 
                    Function[{T\_0, T\_1, T\_2} // Evaluate, 
                      \(a /. coeffs // Expand\) // Evaluate]
	\ \ ]\);
		
		  
		order3Eq = eqOrder[3] /. sol1 /. sol2 // ExpandAll;
		ST21 = 
              Coefficient[order3Eq[[#, 2]] /. expRule1[1], 
                    Exp[I (Subscript[ω, 1]) (Subscript[T, 0])]] & /@ {1, 2};
		ST22 = 
              Coefficient[order3Eq[[#, 2]] /. expRule1[2], 
                    Exp[I (Subscript[ω, 2]) (Subscript[T, 0])]] & /@ {3, 4};
		SCond2 = {ccleftVec[[1]] . ST21 == 0, 
                  ccleftVec[[2]] . ST22 == 0} // ExpandAll;
		SCond2Rule1=Solve[SCond2,{(Subscript[A, 1]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]//ExpandAll;
		
		  
		moduEq = 
              Table[2 I (Subscript[ω, k]) ((Subscript[A, k])^′), {k, 2}] == 
                  (Table[
                            2 I (Subscript[ω, k]) 
                              dt[1][ (Subscript[A, k])[(Subscript[T, 1]), (Subscript[T, 2])]], {k, 2}] /. 
                          SCond1Rule1 /. SCond2Rule1 // 
                      Collect[#, ϵ] &) // Thread;
		Print["\!\(\*
StyleBox[\(The\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(second\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(-\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(order\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(approximate\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(solution\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(:\),
FontColor->RGBColor[0, 0, 1]]\)"];
		Print[
              Table[(Subscript[u, i])[t] == 
                  ((Subscript[u, i]) @@ timeScales /. solRule /. 
                            (ϵ^3) -> 0 /. sol1 /. sol2 /. 
                      displayRule), {i, 2}]];
		If[Or @@ 
                  Table[((Subscript[F, i])/(Subscript[f, i]) /. scaling) === ϵ, {i, 2}],
			Print["\!\(\*
StyleBox[\(where\),
FontColor->RGBColor[0, 0, 1]]\)"];
			\(Print[
                    \*""(Subscript[Λ, i])==((Subscript[f, i])/(2 ((Subscript[ω, i]^2) - (Ω^2))))""]\)
		]; 
		Print["\n\!\(\*
StyleBox[\(The\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(modulation\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(equations\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(:\),
FontColor->RGBColor[0, 0, 1]]\)"];
		Print[moduEq /. displayRule]
]

As an example, we check the case in Section 6.1.1:

\(\(scaling1 = {μ\&^\_n_ -> ϵ\ μ\_n, 
        F\_1 -> ϵ\ f\_1, F\_2 -> ϵ\^2\ f\_2}; \)\)

ResonanceCond1 = {
        (Subscript[ω, 2]) == 2  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
        Ω == (Subscript[ω, 2]) + ϵ (Subscript[σ, 2])};

MMS[scaling1, ResonanceCond1] // Timing

6.1.3 Ω ≈ (Subscript[ω, 1]) and (Subscript[ω, 2]) ≈ 2  (Subscript[ω, 1])

In this case, we have a combination of a primary resonance of the first mode, a superharmonic resonance of the second mode, and a two-to-one internal resonance. In order to bring the effects of damping, forcing, and nonlinearity at the same order, we let

\(\(scaling3 = {μ\&^\_n_ -> ϵ\ μ\_n, 
        F\_1 -> ϵ\^2\ f\_1, F\_2 -> ϵ\ f\_2}; \)\)

To describe quantitatively the nearness of the resonances, we introduce the detuning parameters (Subscript[σ, 1]) and (Subscript[σ, 2]) defined according to

ResonanceCond3 = {
        (Subscript[ω, 2]) == 2  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
        Ω == (Subscript[ω, 1]) + ϵ (Subscript[σ, 2])};

Using MMS, we obtain the second-order approximate solution and the two equations governing the modulation of the complex-valued functions (Subscript[A, k]):

MMS[scaling3, ResonanceCond3] // Timing

6.1.4 \(TraditionalForm
      \`Ω ≈ ω\_1 + ω\_2\) and (Subscript[ω, 2]) ≈ 2  (Subscript[ω, 1])

In this case, we have a combination resonance and a two-to-one internal resonance. We let

\(\(scaling4 = {μ\&^\_n_ -> ϵ\ μ\_n, 
        F\_1 -> ϵ\ f\_1, F\_2 -> ϵ\ f\_2}; \)\)

ResonanceCond4 = {
        (Subscript[ω, 2]) == 2  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
        Ω == 
          (Subscript[ω, 1]) + (Subscript[ω, 2]) + ϵ (Subscript[σ, 2])};

Using MMS, we obtain the second-order approximate solution and the two equations governing the modulation of the complex-valued functions (Subscript[A, k]):

MMS[scaling4, ResonanceCond4] // Timing

6.1.5 Ω ≈ (Subscript[ω, 2]) and (Subscript[ω, 2]) ≈ 3  (Subscript[ω, 1])

In this case, we have a primary resonance of the second mode, a subharmonic resonance of the first mode, and a three-to-one internal resonance. We let

\(\(scaling5 = {μ\&^\_n_ -> ϵ\^2\ μ\_n, 
        F\_1 -> ϵ\ f\_1, F\_2 -> ϵ\^3\ f\_2}; \)\)

ResonanceCond5 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMS, we obtain the second-order approximate solution and the two equations governing the modulation of the complex-valued functions (Subscript[A, k]):

MMS[scaling5, ResonanceCond5] // Timing

6.1.6 \(TraditionalForm
      \`Ω ≈ ω\_2 - 2  ω\_1\) and (Subscript[ω, 2]) ≈ 3  (Subscript[ω, 1])

In this case, we have a primary resonance of the first mode, a superharmonic resonance of the second mode, and a three-to-one internal resonance. We let

\(\(scaling6 = {μ\&^\_n_ -> ϵ\^2\ μ\_n, 
        F\_1 -> ϵ\^3\ f\_1, F\_2 -> ϵ\ f\_2}; \)\)

ResonanceCond6 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 
          (Subscript[ω, 2]) - 2  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 2])};

Using MMS, we obtain the second-order approximate solution and the two equations governing the modulation of the complex-valued functions (Subscript[A, k]):

MMS[scaling6, ResonanceCond6] // Timing

6.1.7 \(TraditionalForm
      \`Ω ≈ ω\_1 + ω\_2\) and (Subscript[ω, 2]) ≈ (Subscript[ω, 1])

In this case, we have a subharmonic resonance of either mode and a one-to-one internal resonance. We let

\(\(scaling7 = {μ\&^\_n_ -> ϵ\ μ\_n, 
        F\_1 -> ϵ\ f\_1, F\_2 -> ϵ\ f\_2}; \)\)

ResonanceCond7 = {
        (Subscript[ω, 2]) == (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
        Ω == 
          (Subscript[ω, 1]) + (Subscript[ω, 2]) + ϵ (Subscript[σ, 2])};

Using MMS, we obtain the second-order approximate solution and the two equations governing the modulation of the complex-valued functions (Subscript[A, k]):

MMS[scaling7, ResonanceCond7] // Timing

6.1.8 \(TraditionalForm
      \`Ω ≈ 2  ω\_2 - ω\_1\) and (Subscript[ω, 2]) ≈ (Subscript[ω, 1])

In this case, we have a primary resonance of either mode and a one-to-one internal resonance. We let

\(\(scaling8 = {μ\&^\_n_ -> ϵ\^2\ μ\_n, 
        F\_1 -> ϵ\^3\ f\_1, F\_2 -> ϵ\^3\ f\_2}; \)\)

ResonanceCond8 = {
        (Subscript[ω, 2]) == (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 
          2  (Subscript[ω, 2]) - (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 2])};

Using MMS, we obtain the second-order approximate solution and the two equations governing the modulation of the complex-valued functions (Subscript[A, k]):

MMS[scaling8, ResonanceCond8] // Timing

6.2 Parametrically Excited Linearly Coupled Systems

6.2.1 Two-Pendulum Oscillator

As an example, we consider a parametrically excited two-pendulum oscillator with a three-to-one internal resonance:

eq621a = {
        ((Subscript[θ, 1])^″)[t] + 
            α ((Subscript[θ, 2])^″)[t] + (Subscript[θ, 1])[t] == 
          (1/2) α ((((Subscript[θ, 2])[t] - (Subscript[θ, 1])[t]))^2) 
              ((Subscript[θ, 2])^″)[t] + 
            (1/6) (Subscript[θ, 1])[t]^3 + 
            α ((Subscript[θ, 2])[t] - (Subscript[θ, 1])[t]) 
              ((Subscript[θ, 2])^′)[t]^2 - 
            2  (Subscript[F, 1]) (Ω^2) (Subscript[θ, 1])[t] 
              Cos[Ω t], 
        ((Subscript[θ, 2])^″)[t] + 
            ((Subscript[θ, 1])^″)[t] + (Subscript[θ, 2])[t] == 
          (1/2) ((((Subscript[θ, 2])[t] - (Subscript[θ, 1])[t]))^2) 
              ((Subscript[θ, 1])^″)[t] + 
            (1/6) (Subscript[θ, 2])[t]^3 - 
            ((Subscript[θ, 2])[t] - (Subscript[θ, 1])[t]) 
              ((Subscript[θ, 1])^′)[t]^2 - 
            2  (Subscript[F, 1]) (Ω^2) (Subscript[θ, 2])[t] 
              Cos[Ω t]};

where

params = {α -> (16/25) + (ϵ^2) (Subscript[σ, 1])};

for small ϵ. Here, the excitation frequency, Ω, is assumed to be close to twice the natural frequency of the second mode, (Subscript[ω, 2]), which is three times the natural frequency of the first mode, (Subscript[ω, 1]). Hence, we define the following lists:

omgList = {(Subscript[ω, 1]), (Subscript[ω, 2])};

ResonanceCond = {
        Ω == 2 (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2]), 
        (Subscript[ω, 2]) == 3 (Subscript[ω, 1])};

omgRule = 
      Solve[ResonanceCond, {Ω, #} // Flatten][[1]] &
             /@ omgList // Reverse

In order to bring the effects of forcing and nonlinearity at the same order, we let

scaling = {(Subscript[F, 1]) -> (ϵ^2) (Subscript[f, 1])};

For a consistent expansion, we first transform eq621a into a system of four first-order equations. To this end, we introduce the two states (Subscript[v, 1])[t] and (Subscript[v, 2])[t] defined by

vel = {((Subscript[θ, 1])^′)[t] -> (Subscript[v, 1])[t], 
        ((Subscript[θ, 2])^′)[t] -> (Subscript[v, 2])[t]};

Substituting the velocity and acceleration terms, using vel, into eq621a and combining the result with vel, we transform eq621a into the following set of four first-order equations:

eq621b = 
      {vel /. Rule -> Equal, eq621a /. D[vel, t] /. vel} // Transpose
           // Flatten

We seek a uniform second-order expansion of the solution of eq621b in the form

solRule = {
        (Subscript[θ, i])_ -> 
          (Sum[(ϵ^j) (Subscript[θ, i, j])[#1, #2, #3], {j, 3}] &)
            , (Subscript[v, i])_ -> 
          (Sum[(ϵ^j) (Subscript[v, i, j])[#1, #2, #3], {j, 3}] &)};

maxOrder = 2;

Using the time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), we express the dependent variables and their time derivatives as

multiScales = {(Subscript[θ, i])_[t] -> (Subscript[θ, i]) @@ timeScales, 
        (Subscript[v, i])_[t] -> (Subscript[v, i]) @@ timeScales, 
        Derivative[1][u_][t] -> dt[1][u @@ timeScales], 
        t -> (Subscript[T, 0])};

Substituting the params, scaling, multiScales, and solRule into eq621b, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq621c = 
      (eq621b /. params /. scaling /. multiScales /. 
                solRule // TrigToExp // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we obtain

eqEps = 
      Rest[Thread[CoefficientList[Subtract @@ #, ϵ] == 0]] &
             /@ eq621c // Transpose;

To place the linear operator on one side and the nonhomogeneous terms on the other side, we define

\(eqOrder[i_] := 
      \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f\_s_ -> 0\) /. 
              u_\_\(k_, 1\) -> u\_\(k, i\))\) == 
          \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f\_s_ -> 0\) /. 
                u_\_\(k_, 1\) -> u\_\(k, i\))\) - 
            \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([i]\)])\) // Thread\)

Using the eqOrder[i] and displayRule, we rewrite eqEps in a concise way as

eqOrder[1] /. displayRuleeqOrder[2] /. displayRuleeqOrder[3] /. displayRule

First-Order Problem: Linear System

The first-order problem, eqOrder[1], consists of a set of coupled linear homogeneous differential equations. Hence, the general solution is the solution of the homogeneous equations. To determine the general solution of the homogeneous set, we rewrite eqOrder[1] as

\(linearSys = \(#[\([1]\)]\ &\)\ /@\ eqOrder[1]; 

    linearSys /. displayRule\)

Next, we seek a solution of the linearSys in the form

assumedForm = {
        (Subscript[θ, i_, 1]) -> ((Subscript[c, i]) (E^(I ω #)) &), 
        (Subscript[v, i_, 1]) -> ((Subscript[d, i]) (E^(I ω #)) &)};

Substituting assumedForm into linearSys and collecting coefficients of (Subscript[c, i]) and (Subscript[d, i]), we obtain the coefficient matrix as

coefMat = 
      Outer[D, (E^(-I ω (Subscript[T, 0]))) linearSys /. assumedForm // 
          Expand, {(Subscript[c, 1]), (Subscript[d, 1]), (Subscript[c, 2]), (Subscript[d, 2])}]

whose adjoint is defined by

hermitian[mat_?MatrixQ] := mat /. conjugateRule // Transpose

The natural frequencies of the two modes involved are

evals = 
      Cases[ω /. Solve[Det[coefMat] == 0, ω], n_ /; n > 0]

Hence, the right and left eigenvectors of coefMat corresponding to evals are

rightVec = 
      #/#[[1]] & /@ 
        (NullSpace[coefMat /. ω -> #][[1]] & /@ evals)

leftVec = 
      NullSpace[hermitian[coefMat] /. ω -> #][[1]] & /@ 
        evals

whose complex conjugates are

ccleftVec = leftVec /. conjugateRule

Based on the rightVec, we can assume the general solution form as

sol1Form = 
      Transpose[rightVec] . {(Subscript[A, 1])[(Subscript[T, 1]), (Subscript[T, 2])] (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))), 
          (Subscript[A, 2])[(Subscript[T, 1]), (Subscript[T, 2])] (E^(I (Subscript[ω, 2]) (Subscript[T, 0])))}

where (Subscript[ω, 1]) = (Sqrt[5])/3 and (Subscript[ω, 2]) = (Sqrt[5]). Therefore, the general solution of eqOrder[1] can be expressed in pure function form as

\(sol1 = {θ\_\(1, 1\), v\_\(1, 1\), θ\_\(2, 1\), 
            v\_\(2, 1\)} -> 
          \((\(Function[{T\_0, T\_1, T\_2}, # + \((# /. conjugateRule)\) // 
                    Evaluate]\ &\)\ /@\ sol1Form)\) // Thread; 

    sol1 /. displayRule\)

where (Subscript[A, 1]) and (Subscript[A, 2]) are to be determined from the solvability conditions at the next levels of approximation.

Second-Order Problem

Substituting the first-order solution sol1 into the second-order problem, eqOrder[2], yields

\(order2Eq = eqOrder[2] /. sol1 // ExpandAll; 

    order2Eq /. displayRule\)

For a uniform expansion, we choose (Subscript[D, 1]) (Subscript[A, 1]) and (Subscript[D, 1]) (Subscript[A, 2]) to eliminate the secular or small-divisor terms. To accomplish this, we first convert the small-divisor terms into secular terms using the rule

expRule1[i_] := 
      Exp[a_] :> 
        Exp[Expand[a /. omgRule[[i]]] /. (ϵ^2) (Subscript[T, 0]) -> (Subscript[T, 2])]

To eliminate the terms that produce secular terms (i.e., determine the solvability conditions) from order2Eq, we collect the terms proportional to (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))) and (E^(I (Subscript[ω, 2]) (Subscript[T, 0]))) and obtain

ST11 = 
      Coefficient[order2Eq[[#, 2]] /. expRule1[1], 
            (E^(I (Subscript[ω, 1]) (Subscript[T, 0])))] & /@ Range[4]

ST12 = 
      Coefficient[order2Eq[[#, 2]] /. expRule1[2], 
            (E^(I (Subscript[ω, 2]) (Subscript[T, 0])))] & /@ Range[4]

Then, the solvability conditions demand that ST11 and ST12 be orthogonal to every solution of the corresponding adjoint problem, namely, the components of ccleftVec. The result is

SCond11 = ccleftVec[[1]] . ST11 == 0

SCond12 = ccleftVec[[2]] . ST12 == 0

Solving SCond11 and SCond12 for (Subscript[D, 1]) (Subscript[A, 1]) and (Subscript[D, 1]) (Subscript[A, 2]), respectively, yields

SCond11Rule=Solve[SCond11,(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]

SCond12Rule=Solve[SCond12,(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]

SCond1Rule = Join[SCond11Rule, SCond12Rule];

whose complex conjugate is

ccSCond1Rule = SCond1Rule /. conjugateRule;

Substituting SCond1Rule and ccSCond1Rule into order2Eq, we obtain

\(order2Eqm = \(order2Eq /. SCond1Rule\) /. ccSCond1Rule // ExpandAll; 

    order2Eqm /. displayRule\)

which is a set of homogeneous equations. Then, we express the solution of the second-order equations in pure function form as

sol2 = {(Subscript[θ, 1, 2]) -> (0 &), (Subscript[v, 1, 2]) -> (0 &), 
        (Subscript[θ, 2, 2]) -> (0 &), (Subscript[v, 2, 2]) -> (0 &)};

Third-Order Problem

Substituting the first- and second-order solutions into the third-order problem, eqOrder[3], yields

order3Eq = eqOrder[3] /. sol1 /. sol2 // ExpandAll;

Substituting the expRule1 into the right-hand sides of order3Eq and collecting the terms that could produce secular terms, the terms proportional to (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))) and (E^(I (Subscript[ω, 2]) (Subscript[T, 0]))), we have

\(ST21 = 
      \(Coefficient[#[\([2]\)] /. expRule1[1], E\^\(I\ ω\_1\ T\_0\)]\ &
          \)\ /@\ order3Eq; 
ST21 /. displayRule\)

\(ST22 = 
      \(Coefficient[#[\([2]\)] /. expRule1[2], E\^\(I\ ω\_2\ T\_0\)]\ &
          \)\ /@\ order3Eq; 
ST22 /. displayRule\)

Then, solvability conditions demand that ST21 and ST22 be orthogonal to their corresponding adjoints. The result is

SCond21=Solve[ccleftVec[[1]] . ST21 == 0,(Subscript[A, 1]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]/.{(Subscript[ω, 1]) -> (Sqrt[5])/3, (Subscript[ω, 2]) -> (Sqrt[5])}//ExpandAll;
SCond21 /. displayRule

SCond22=Solve[ccleftVec[[2]] . ST22 == 0,(Subscript[A, 2]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]/.{Ω -> 2 (Sqrt[ 5]), (Subscript[ω, 1]) -> (Sqrt[5])/3, 
              (Subscript[ω, 2]) -> (Sqrt[5])}//ExpandAll;
SCond22 /. displayRule

SCond2Rule = Join[SCond21, SCond22];

Reconstitution

Using the method of reconstitution, \(TraditionalForm
      \`A\_k\^′ = 
        ϵ\ \(D\_1\) A\_k + ϵ\^2\ \(D\_2\) A\_k + \ 
          ⋯\), we combine the partial-differential equations SCond1Rule and SCond2Rule into the following two ordinary-differential equations governing the modulation of the complex-valued functions (Subscript[A, k]) :

\(moduEq = 
      Table[2\ I\ ω\_k\ A\_k\^′, {k, 2}] == 
          \((\(Table[
                  2\ I\ ω\_k\ \(dt[1]\)[\ A\_k[T\_1, T\_2]], {k, 
                    2}] /. SCond1Rule\) /. SCond2Rule)\) // Thread; 

    moduEq /. displayRule\)

6.2.2 The function MMSC

According to the procedures described in the previous section, we build a function named MMSC (Method of Multiple Scales for Linearly Coupled Systems) to automate the process.

MMSC[eqs_List, depVar_List, scaling_List, ResonanceCond : {__Equal}]:=
   Module[{},
		omgList = {(Subscript[ω, 1]), (Subscript[ω, 2])};
     omgRule = 
              Solve[ResonanceCond, {Ω, #} // Flatten][
                      [1]] & /@ omgList // Reverse;
		mydepVar = {(Subscript[u, 1]), (Subscript[v, 1]), (Subscript[u, 2]), (Subscript[v, 2])};
		solRule = {
                (Subscript[u, i])_ -> 
                  (Sum[ϵ^j (Subscript[u, i, j])[#1, #2, #3], {j, 3}] &)
                    , (Subscript[v, i])_ -> 
                  (Sum[ϵ^j (Subscript[v, i, j])[#1, #2, #3], {j, 3}] &)}
                ;
		maxOrder = 2;
     multiScales = {(Subscript[u, i])_[t] -> (Subscript[u, i]) @@ timeScales, 
                (Subscript[v, i])_[t] -> (Subscript[v, i]) @@ timeScales, 
                Derivative[1][u_][t] -> dt[1][u @@ timeScales], 
                t -> (Subscript[T, 0])};
		eqa = (
                  eqs /. Thread[depVar -> mydepVar] /. scaling /. 
                          multiScales /. solRule // TrigToExp // 
                    ExpandAll) /. ϵ^(n_ /; n > 3) -> 0;
		eqEps = 
              Rest[Thread[
                        CoefficientList[Subtract @@ #, ϵ] == 0]] &
                     /@ eqa // Transpose;
	  \(eqOrder[i_] := 
              \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f\_s_ -> 0\) /. 
                      u_\_\(k_, 1\) -> u\_\(k, i\))\) == 
                  \((\(\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. 
                          f\_s_ -> 0\) /. u_\_\(k_, 1\) -> u\_\(k, i\))\) - 
                    \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([i]\)])\) // Thread\);
	
		  		
		linearSys = #[[1]] & /@ eqOrder[1];
		assumedForm = {
                (Subscript[u, i_, 1]) -> ((Subscript[c, i]) E^(I ω #) &), 
                (Subscript[v, i_, 1]) -> ((Subscript[d, i]) E^(I ω #) &)};
		coefMat = 
              Outer[D, 
                E^(-I ω (Subscript[T, 0])) linearSys /. assumedForm // 
                  Expand, {(Subscript[c, 1]), (Subscript[d, 1]), (Subscript[c, 2]), (Subscript[d, 2])}];
		hermitian[mat_?MatrixQ] := mat /. conjugateRule // Transpose;
		evals = 
              Cases[ω /. Solve[Det[coefMat] == 0, ω], 
                n_ /; n > 0];
		values = 
              Append[omgRule[[1]] /. ϵ -> 0 /. 
                  (Subscript[ω, 1]) -> evals[[1]], 
                (Subscript[ω, 1]) -> evals[[1]]];
		rightVec = 
              #/#[[1]] & /@ 
                (NullSpace[coefMat /. ω -> #][[1]] & /@
                     evals);
     leftVec = 
              NullSpace[hermitian[coefMat] /. ω -> #][[1]] &
                   /@ evals;
		ccleftVec = leftVec /. conjugateRule;
		\(order1Eq = 
              eqOrder[1] /. u_\_\(i_, 1\) -> \((u\_\(i, 1\)[#1]\ &)\)\);
		sol1p = 
              DSolve[
                      order1Eq, {(Subscript[u, 1, 1])[(Subscript[T, 0])], (Subscript[v, 1, 1])[(Subscript[T, 0])], 
                        (Subscript[u, 2, 1])[(Subscript[T, 0])], (Subscript[v, 2, 1])[(Subscript[T, 0])]}, (Subscript[T, 0])][
                    [1]] /. C[_] -> 0 // Simplify;
		fRule = {
                (Subscript[f, i])_ -> 
                  2 (Subscript[Λ, i]) 
                    ((Subscript[ω, i]^2) - Ω^2)};
		sol1Form = 
              (#[[2]] & /@ sol1p) + 
                ((# + (# /. conjugateRule)) & /@ 
                    (Transpose[rightVec] . {
                          (Subscript[A, 1])[(Subscript[T, 1]), (Subscript[T, 2])] E^(I (Subscript[ω, 1]) (Subscript[T, 0])), 
                          (Subscript[A, 2])[(Subscript[T, 1]), (Subscript[T, 2])] E^(I (Subscript[ω, 2]) (Subscript[T, 0]))}))
                  ;
		sol1 = {(Subscript[u, 1, 1]), (Subscript[v, 1, 1]), (Subscript[u, 2, 1]), (Subscript[v, 2, 1])} -> 
                  (Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, # // Evaluate] & /@ 
                      sol1Form) // Thread;
				
		  
		order2Eq = eqOrder[2] /. sol1 // ExpandAll;
		expRule1[i_] := 
              Exp[a_] :> 
                Exp[Expand[a /. omgRule[[i]]] /. 
                    ϵ^n_.  (Subscript[T, 0]) :> timeScales[[n + 1]]];
		ST11 = 
              Coefficient[#[[2]] /. expRule1[1], 
                    E^(I (Subscript[ω, 1]) (Subscript[T, 0]))] & /@ order2Eq;
		ST12 = 
              Coefficient[#[[2]] /. expRule1[2], 
                    E^(I (Subscript[ω, 2]) (Subscript[T, 0]))] & /@ order2Eq;
		SCond1 = {ccleftVec[[1]] . ST11 == 0, 
                ccleftVec[[2]] . ST12 == 0};
		SCond1Rule1=Solve[SCond1,{(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]//ExpandAll;
		  
		sigRule = Solve[ResonanceCond, (Subscript[σ, 2])][[1]];
		expRule2 = 
              Exp[a_] :> 
                Exp[a /. {(Subscript[T, 1]) -> ϵ (Subscript[T, 0]), 
                          (Subscript[T, 2]) -> ϵ^2 (Subscript[T, 0])} /. sigRule // Expand]
                  ;
		SCond1Rule2 = SCond1Rule1 /. expRule2;
		ccSCond1Rule2 = SCond1Rule2 /. conjugateRule;
		order2Eqm = 
              order2Eq /. SCond1Rule2 /. ccSCond1Rule2 // ExpandAll;
		\(If[Union[\(#[\([2]\)]\ &\)\ /@\ order2Eqm] === {0}, 
	\ \ \ 	
              sol2 = {u\_\(1, 2\) -> \((0\ &)\), v\_\(1, 2\) -> \((0\ &)\), 
                  u\_\(2, 2\) -> \((0\ &)\), v\_\(2, 2\) -> \((0\ &)\)}, 

              		\ \ \ 
              Return["The particular solution of the second order equations needs to be solved, which is not included in this function."]

              		 (*
\ \ \ 		
                basicH = 
                  Table[{A\_i[T\_1, T\_2]\ E^\((I\ ω\_i\ T\_0)\), 
                        A\&_\_i[T\_1, T\_2]\ 
                          E^\((\(-I\)\ ω\_i\ T\_0)\)}, {i, 2}] // 
                    Flatten; 
\ \ \ 		
                collectForm = 
                  Join[basicH, 
                    If[list1 = 
                        List\ @@\ 
                            \(Plus\ @@\ 
                              \((\(#[\([2]\)]\ &\)\ /@\ eqOrder[1])\)\) /. 
                          c_\ f\_i_\ E^a_ -> Λ\_i\ E^a; 
                      Head[list1] === List, list1, {}]]; 
\ \ \ 		
                possibleTerms = 
                  Join[collectForm, 
                    If[list1 = 
                        List\ @@\ 
                            \(Plus\ @@\ 
                              \((\(#[\([2]\)]\ &\)\ /@\ eqOrder[2] /. 
                                  u_\_\(i_, j_\) -> \((0\ &)\))\)\) /. 
                          c_\ f\_i_\ E^a_ -> E^a; Head[list1] === List, 
                      list1, {}], 
                    \(Outer[Times, collectForm, collectForm] // Flatten\) // 
                      Union]; 
	\ \ \ 	
                ResonantTerms[i_] := 
                  \(\((\(# /. {a_ /; a =!= 0 -> 1}\ &\)\ /@\ 
                            \((\(E^\((\(-I\)\ ω\_i\ T\_0)\)\ 
                                    possibleTerms /. expRule1[i]\) /. 
                                Exp[_\ T\_0 + _. ] -> 0)\))\)\ 
                        possibleTerms // Union\) // Rest; 
\ \ \ 	\ \ \ 
                RT = {Array[ResonantTerms, 2], 
                      Array[ResonantTerms, 2] /. conjugateRule} // Flatten; 

                \ \ \ 		NRT = Complement[possibleTerms, RT]; 
\ \ \ 		
                rRule = 
                  Table[MapIndexed[
                        r\_\(i, #2[\([1]\)]\) -> 
                            Coefficient[order2Eqm[\([i, 2]\)], #1]\ &, RT], {
                        i, 4}] // Flatten; 
\ \ \ 		
                list1 = Join[ccleftVec, Conjugate[ccleftVec]]; 
\ \ \ \ \ 	
                r4Rule = 
                  Table[\(Solve[
                          list1[\([j]\)] . Table[r\_\(i, j\), {i, 4}] == 0, 
                          r\_\(4, j\)]\)[\([1]\)], {j, 4}] // Flatten; 

                \ \ \ 		
                RTsymbolList[i_] = 
                  Table[Γ\_\(i, j\), {j, Length[RT]}]; 

                \ \ \ 		
                NRTsymbolList[i_] = 
                  Table[Λ\_\(i, j\), {j, Length[NRT]}]; 

                		\ \ \ 
                sol2Form = {u\_\(1, 2\), v\_\(1, 2\), u\_\(2, 2\), 
                        v\_\(2, 2\)} -> 
                      \((\(Function[{T\_0, T\_1, T\_2}, 
                              RTsymbolList[#] . RT + 
                                  NRTsymbolList[#] . NRT // Evaluate]\ &\)\ /@
                          \ Range[4])\) // Thread; 
\ \ \ 		
                coef1 = 
                  \(Solve[
                          \(Coefficient[Subtract\ @@\ # /. sol2Form, NRT] == 
                                  0 /. Exp[_\ T\_0 + _. ] -> 0\ &\)\ /@\ 
                            order2Eqm, Array[NRTsymbolList, 4] // Flatten]\)[
                        \([1]\)] /. values // ExpandAll; 
\ \ \ 		
                eq1 = Table[
                        MapIndexed[
                            Coefficient[
                                  order2Eqm[\([k, 1]\)] /. sol2Form, #1] == 
                                r\_\(k, #2[\([1]\)]\)\ &, RT] /. 
                          Exp[_\ T\_0 + _. ] -> 0, {k, 4}] /. r4Rule // 
                    Flatten; 
\ \ \ \ \ \ \ \ 
                coef2 = 
                  \(\((\(Solve[eq1, Join\ @@\ Array[RTsymbolList, 4]]\)[
                              \([1]\)] // Simplify)\) /. values\) /. rRule // 
                    ExpandAll; 
	\ \ \ 	
                sol2 = sol2Form /. 
                    Function[{T\_0, T\_1, T\_2}, a_] :> 
                      Function[{T\_0, T\_1, T\_2} // Evaluate, 
                        \(\(a /. coef1\) /. coef2 // Expand\) // Evaluate]
                        \ \ \ *) 
	\ \ ]\);
		
		  
		order3Eq = eqOrder[3] /. sol1 /. sol2 // ExpandAll;
		ST21 = 
              Coefficient[#[[2]] /. expRule1[1], 
                    Exp[I (Subscript[ω, 1]) (Subscript[T, 0])]] & /@ order3Eq;
		ST22 = 
              Coefficient[#[[2]] /. expRule1[2], 
                    Exp[I (Subscript[ω, 2]) (Subscript[T, 0])]] & /@ order3Eq;
		SCond2 = {ccleftVec[[1]] . ST21 == 0, 
                  ccleftVec[[2]] . ST22 == 0} // ExpandAll;
		SCond2Rule1=Solve[SCond2,{(Subscript[A, 1]^((0, 1)))[(Subscript[T, 1]), 
                              (Subscript[T, 2])],(Subscript[A, 2]^((0, 1)))[(Subscript[T, 1]), 
                              (Subscript[T, 2])]}][[1]]/.values//ExpandAll;
		
		  
		moduEq = 
              Table[2 I (Subscript[ω, k]) ((Subscript[A, k])^′), {k, 2}] == 
                  (Table[
                            2 I (Subscript[ω, k]) 
                              dt[1][ (Subscript[A, k])[(Subscript[T, 1]), (Subscript[T, 2])]], {k, 2}] /. 
                          (SCond1Rule1 /. values) /. SCond2Rule1 // 
                      Collect[#, ϵ] &) // Thread;
		Print["\!\(\*
StyleBox[\(The\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(second\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(-\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(order\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(approximate\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(solution\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(:\),
FontColor->RGBColor[0, 0, 1]]\)"];
		Print[
              Table[(Subscript[u, i])[t] == 
                    ((Subscript[u, i]) @@ timeScales /. solRule /. 
                              ϵ^3 -> 0 /. sol1 /. sol2 /. 
                        displayRule), {i, 2}] /. Thread[mydepVar -> depVar]]
              ;
		If[Or @@ 
                  Table[((Subscript[F, i])/(Subscript[f, i]) /. scaling) === ϵ, {i, 2}],
			Print["\!\(\*
StyleBox[\(where\),
FontColor->RGBColor[0, 0, 1]]\)"];
			\(Print[
                    \*""(Subscript[Λ, i])==(Subscript[f, i])/(2((Subscript[ω, i]^2)-Ω^2))""]\)
		]; 
		Print["\n\!\(\*
StyleBox[\(The\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(modulation\),
FontColor->RGBColor[0, 0, 1]]\) \!\(\*
StyleBox[\(equations\),
FontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\(:\),
FontColor->RGBColor[0, 0, 1]]\)"];
		Print[moduEq /. displayRule]
]

As an example, we check the case in Section 6.2.1:

eqtest={
		(Subscript[θ, 1]^′)[t]==(Subscript[v, 1])[t],
		(Subscript[v, 1]^′)[t]+α (Subscript[v, 2]^′)[t]+(Subscript[θ, 1])[t]==-2 (Ω^2) Cos[t Ω] (Subscript[F, 1]) 
                    (Subscript[θ, 1])[t]+(1/6) (Subscript[θ, 1])[t]^3+α (Subscript[v, 2])[t]^2 
                    (-(Subscript[θ, 1])[t] + (Subscript[θ, 2])[t])+(1/2) α (((-(Subscript[θ, 1])[t] + (Subscript[θ, 2])[t]))^2) (Subscript[v, 2]^′)[t],
		(Subscript[θ, 2]^′)[t]==(Subscript[v, 2])[t],
		(Subscript[v, 2]^′)[t]+(Subscript[v, 1]^′)[t]+(Subscript[θ, 2])[t]==-2 (Ω^2) Cos[t Ω] (Subscript[F, 1]) 
                    (Subscript[θ, 2])[t]+(1/6) (Subscript[θ, 2])[t]^3-(Subscript[v, 1])[t]^2 (-(Subscript[θ, 1])[t] + (Subscript[θ, 2])[t])+(1/2) (((-(Subscript[θ, 1])[t] + (Subscript[θ, 2])[t]))^2) (Subscript[v, 1]^′)[t]
	}/.{α -> (16/25) + (ϵ^2) (Subscript[σ, 1])};

scaling1 = {(Subscript[F, 1]) -> (ϵ^2) (Subscript[f, 1])};

ResonanceCond1 = {(Subscript[ω, 2]) == 3  (Subscript[ω, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

MMSC[eqtest, {(Subscript[θ, 1]), (Subscript[v, 1]), (Subscript[θ, 2]), (Subscript[v, 2])}, scaling1, 
        ResonanceCond1] // Timing