Chapter 5
Higher-Order Approximations for Systems with Internal Resonances

Off[General::spell1]

5.1 Euler-Lagrange Equations

In this chapter, we use different methods to determine approximate solutions of nonlinear systems possessing internal resonances to orders higher than the order at which the influence of the internal resonance first appears. To describe the methods with minimum algebra, we consider the free oscillations of a two-degree-of-freedom conservative system possessing a two-to-one internal resonance and having simple quadratic nonlinearities. In particular, we consider a system governed by the Lagrangian

Lagrangian = 
      (1/2) ((Subscript[u, 1])^′)[t]^2 + (1/2) ((Subscript[u, 2])^′)[t]^2 - 
        (1/2) (Subscript[ω, 1]^2) (Subscript[u, 1])[t]^2 - 
        (1/2) (Subscript[ω, 2]^2) (Subscript[u, 2])[t]^2 + δ (Subscript[u, 1])[t]^2 (Subscript[u, 2])[t];

where (Subscript[ω, 2]) ≈ 2  (Subscript[ω, 1]). Writing down the Euler-Lagrange equations corresponding to the Lagrangian, we obtain the following second-order equations of motion:

EOM1 = 
      Table[D[D[Lagrangian, ((Subscript[u, i])^′)[t]], t] - 
            D[Lagrangian, (Subscript[u, i])[t]] == 0, {i, 2}]

Approximate solutions of general two-degree-of-freedom system with quadratic and cubic nonlinearities having a two-to-one internal resonance can be obtained by replacing the Lagrangian with the Lagrangian corresponding to the general system.

As discussed in the next section, treating this second-order form, EOM1, of the governing equations by using the method of multiple scales may lead to results that violate the conservative nature of the system under consideration. To determine an approximate solution that preserves the conservative nature of the system, we treat a first-order form of the governing equations. To transform EOM1 into a system of first-order equations, we form a modified Lagrangian by introducing two more states (Subscript[v, 1])[t] and (Subscript[v, 2])[t] such that

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

Then, the Lagrangian becomes

Lag1 = Lagrangian /. velRule

Substituting Lag1 into Euler-Lagrange equations, we obtain the following first-order equations

eq51a = 
      Table[D[D[Lag1, (Subscript[v, i])[t]], t] - D[Lag1, (Subscript[u, i])[t]] == 0, {i, 2}]

Using the velRule and eq51a, we obtain the following first-order equations of motion:

EOM2 = 
      {velRule /. Rule -> Equal, eq51a} // Transpose // Flatten

The corresponding linear system of EOM2 can be defined as \(TraditionalForm\`x\& .  = A\ x\), where x is the state vector; that is,

states = Table[{(Subscript[u, i])[t], (Subscript[v, i])[t]}, {i, 2}] // Flatten

Using these states, we find that the coefficient matrix A of the linear part of EOM2 is given by

matrixA = 
      Outer[D, #[[2]]& /@ Solve[EOM2, D[states, t]][[1]] /. 
          δ -> 0, states]

whose adjoint is defined by

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

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

The eigenvalues and eigenvectors of matrixA can be obtained as

matrixA // Eigensystem

Hence, the right eigenvectors of matrixA corresponding to the eigenvalues I (Subscript[ω, 1]) and I (Subscript[ω, 2]), respectively, are

rightVec = Eigenvectors[matrixA][[{2, 4}]]

Moreover, the left eigenvectors of matrixA corresponding to the eigenvalues I (Subscript[ω, 1]) and I (Subscript[ω, 2]), respectively, are

leftVec = Eigenvectors[hermitian[matrixA]][[{1, 3}]]

whose complex conjugate is

ccleftVec = leftVec /. conjugateRule

5.2 Method of Multiple Scales

We use the method of multiple scales to obtain a second-order uniform expansion of the solution of EOM1 in Section 5.2.1, EOM2 in Section 5.2.2, and the corresponding first-order complex-valued equations of EOM1 in Section 5.2.3. Again, we need the three time scales (Subscript[T, 0]) = t, (Subscript[T, 1]) = ϵ t,  and (Subscript[T, 2]) = (ϵ^2) t, which we symbolize by

Needs["Utilities`Notation`"]

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

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

In terms of the new time scales, we can express the time derivatives as

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

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}; 
    \)\)

5.2.1 Second-Order Real-Valued System

Using method of multiple scales, we assume that the solution of EOM1 can be expressed in the form

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

multiScales = {(Subscript[u, i])_[t] -> (Subscript[u, i])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
        Derivative[n_][(Subscript[u, i])_][t] :> 
          dt[n][(Subscript[u, i])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]]};

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

\(eq521a = 
      \((\(EOM1 /. multiScales\) /. solRule // ExpandAll)\) /. 
        ϵ\^\(n_ /; n > 3\) -> 0; 
eq521a /. displayRule\)

Equating coefficients of like powers of ϵ, we obtain

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

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

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

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

eqOrder[1] /. displayRule // TableFormeqOrder[2] /. displayRule // TableFormeqOrder[3] /. displayRule // TableForm

First-Order Equations: Linear System

To obtain the solution of eqOrder[1] by using DSolve, we transform eqOrder[1] into a set of ordinary-differential equations as

order1Eq = eqOrder[1] /. (Subscript[u, i_, j_]) -> ((Subscript[u, i, j])[#1] &)

The particular solution of order1Eq can be expressed as

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

Hence, we write the first-order solution in function form as

\(sol1 = 
      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 // Evaluate], {i, 2}]\)

Second-Order Equations

Substituting sol1 into eqOrder[2] yields

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

    order2Eq /. displayRule\)

To describe quantitatively the nearness of the two-to-one internal resonance (Subscript[ω, 2]) ≈ 2  (Subscript[ω, 1]), we introduce the detuning parameter σ defined by

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

In eliminating secular terms, we need to express (Subscript[ω, 2]) in terms of (Subscript[ω, 1]) in some of the equations and (Subscript[ω, 1]) in terms of (Subscript[ω, 2]) in the other equations. To accomplish this, we let

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

omgRule = 
      Solve[ResonanceCond, Drop[omgList, {#}]][[1]] & /@
         {1, 2}

To convert small-divisor terms into secular terms, we define the rule

expRule1[j_] := 
      Exp[arg_] :> 
        Exp[Expand[arg /. omgRule[[j]]] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]

Substituting expRule1[i] into the right-hand sides of order2Eq and collecting the coefficients of (E^(I (Subscript[ω, i]) (Subscript[T, 0]))), we obtain the solvability conditions as

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

Because the system is conservative, the modulation equations given by SCond1 must be derivable from a Lagrangian. These equations are the Euler-Lagrange equations corresponding to the Lagrangian

Lag=\(E\^\(I\ T\_1\ σ\)\ δ\ A\_2[T\_1, T\_2]\ 
            A\&_\_1[T\_1, T\_2]\^2\)+\(E\^\(\(-I\)\ T\_1\ σ\)\ δ\ A\_1[T\_1, T\_2]\^2\ 
            A\&_\_2[T\_1, T\_2]\)-I (Subscript[ω, 1]) \(A\&_\_1[T\_1, T\_2]\) (Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]-I (Subscript[ω, 2]) \(A\&_\_2[T\_1, T\_2]\) (Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]+
	I (Subscript[ω, 1]) (Subscript[A, 1])[(Subscript[T, 1]), (Subscript[T, 2])] (Subscript[\(A\&_\), 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]+I (Subscript[ω, 2]) (Subscript[A, 2])[(Subscript[T, 1]), (Subscript[T, 2])] (Subscript[\(A\&_\), 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])];

because

alSCond1=Table[-D[D[Lag,(Subscript[\(A\&_\), k]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]],(Subscript[T, 1])]+\(D[Lag, A\&_\_k[T\_1, T\_2]]\)==0,{k, 2}]

Equations alSCond1 are identical to equations SCond1 as shown below

alSCond1 === SCond1

Alternatively, we can determine the modulation equations by using the time-averaged Lagrangian. To this end, we first define the rule

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

Substituting multiScales, solRule, sol1, and expRule2 into the Lagrangian and selecting the slow-varying terms, we obtain the following time-averaged Lagrangian:

\(TAL1 = 
      Select[\((
              \(\((\(Lagrangian /. multiScales\) /. solRule // Expand)\) /. 
                    ϵ\^\(n_ /; n > 3\) -> 0\) /. sol1 // Expand)
              \) /. expRule2, FreeQ[#, T\_0]\ &] /. ϵ -> 1; 

    TAL1 /. displayRule\)

One can easily verify that TAL1 is identical to Lag. In fact,

TAL1 == Lag

Next, we rewrite SCond1 as

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]]

whose complex conjugate is

ccSCond1Rule1 = SCond1Rule1 /. conjugateRule;

To solve for the particular solution of order2Eq, we express the detuning parameter σ in terms of the (Subscript[ω, i]); that is,

sigRule = Solve[ResonanceCond, σ][[1]]

Moreover, we express the scale (Subscript[T, 1]) in terms of the scale (Subscript[T, 0]) by using the rule

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

With the solvability conditions and the expRule3, order2Eq becomes

\(order2Eqm = 
      \(\(order2Eq /. SCond1Rule1\) /. ccSCond1Rule1\) /. expRule3; 

    order2Eqm /. displayRule\)

Again, to use DSolve, we first transform order2Eqm into a set of ordinary-differential equations and then obtain the particular solution of order2Eqm as

sol2p = 
      DSolve[order2Eqm /. (Subscript[u, i_, j_]) -> ((Subscript[u, i, j])[#1] &), {
                (Subscript[u, 1, 2])[(Subscript[T, 0])], (Subscript[u, 2, 2])[(Subscript[T, 0])]}, (Subscript[T, 0])][[1]] /. 
          C[_] -> 0 // Simplify

The second-order solution can be expressed in pure function form as

\(sol2 = 
      Table[u\_\(i, 2\) -> 
          Function[{T\_0, T\_1, T\_2}, 
            \(B\_i[T\_1, T\_2]\ Exp[I\ ω\_i\ T\_0] + 
                  B\&_\_i[T\_1, T\_2]\ Exp[\(-I\)\ ω\_i\ T\_0] + 
                  sol2p[\([i, 2]\)] // Expand\) // Evaluate], {i, 2}]; 

    sol2 /. displayRule\)

As shown below, the obtained modulation equations will violate the conservative nature of system unless we include an appropriate part of the homogeneous solutions of the second-order problem. This is the reason why we included the homogeneous solutions with the undetermined functions (Subscript[B, i])[t]. They will be determined so as to preserve the conservative nature of the system.

Third-Order Equations

Substituting sol1 and sol2 into eqOrder[3] yields

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

Substituting expRule1[i] into the right-hand sides of order3Eq and collecting the coefficients of (E^(I (Subscript[ω, i]) (Subscript[T, 0]))),  we obtain the solvability conditions as

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

    SCond2 /. displayRule\)

Using the solvability conditions at second order, SCond1Rule1, we obtain

\(SCond1Rule2 = \(D[SCond1Rule1, T\_1] /. SCond1Rule1\) /. ccSCond1Rule1; 
    
SCond1Rule2 /. displayRule\)

Hence, we rewrite the solvability conditions at third order as

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]]/.SCond1Rule2//ExpandAll;
SCond2Rule1 /. displayRule

Reconstitution

Using SCond1Rule1 and SCond2Rule1, we reconstitute the modulation equations as

\(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 // Expand\) // 
              Collect[#, ϵ]&)\) // Thread; 

    moduEq /. displayRule\)

Without including the (Subscript[B, i]) in the second-order solution, the following two coefficients:

\(Coefficient[moduEq[\([1, 2]\)], 
          A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_2[T\_1, T\_2]] /. 
        ω\_2 -> 2\ ω\_1 // Simplify\)

\(Coefficient[moduEq[\([2, 2]\)], 
          A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_1[T\_1, T\_2]] /. 
        ω\_2 -> 2\ ω\_1 // Simplify\)

are not the same, which violates the conservative nature of the system in that the modulation equations must be derivable from a Lagrangian. Using the flexibility given by the homogeneous parts of the solutions of the second-order problem, we choose the functions (Subscript[B, i])[t] so that these coefficients are the same. To this end, we assume that

\(\(BForm = {
        B\_1 -> Function[{T\_1, T\_2}, 
            c\_1\ E\^\(I\ T\_1\ σ\)\ A\&_\_1[T\_1, T\_2]\ 
              A\_2[T\_1, T\_2]\ ], 
        B\&_\_1 -> 
          Function[{T\_1, T\_2}, 
            c\_1\ E\^\(\(-I\)\ T\_1\ σ\)\ A\_1[T\_1, T\_2]\ 
              A\&_\_2[T\_1, T\_2]\ ], 
        B\_2 -> Function[{T\_1, T\_2}, 
            c\_2\ E\^\(\(-I\)\ T\_1\ σ\)\ A\_1[T\_1, T\_2]\^2\ ]}; 
    \)\)

where (Subscript[c, 1]) and (Subscript[c, 2]) are undetermined constants. Substituting BForm and the solvability conditions into the right-hand sides of moduEq yields

\(moduEqMod = 
      \((\(Expand[\(\(# /. BForm\) /. SCond1Rule1\) /. ccSCond1Rule1]\ &\)\ /@
            \ moduEq)\) // ExpandAll; 
moduEqMod /. displayRule\)

In order that the moduEqMod be derivable from a Lagrangian, the following condition must be satisfied:

\(\(Coefficient[moduEqMod[\([1, 2]\)], 
              A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_2[T\_1, T\_2]] - 
            Coefficient[moduEqMod[\([2, 2]\)], 
              A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_1[T\_1, T\_2]] == 
          0 /. ϵ -> 1\) /. ω\_2 -> 2  ω\_1\)

This provides a compatibility condition. However, one extra constant is floating !!!

To obtain a consistent second-order uniform expansion, we apply the method of multiple scales to the system of four first-order real-valued equations EOM2 in Section 5.2.2 and to a system of two first-order complex-valued equations in Section 5.2.3.

5.2.2 First-Order Real-Valued System

As shown in the preceding section, treating the second-order form of the equations governing a system may lead to inconsistent results unless appropriate parts of the homogeneous solutions of the perturbation equations are included. Moreover, the consistent expansions might contain arbitrary constants that need to be chosen judiciously. In this section, we show that treatment of an equivalent set of first-order real-valued equations, namely EOM2, leads to consistent results without floating constants.

We seek a second-order uniform expansion of the solution of EOM2 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}] &)};

Expressing the time derivative in EOM2 in terms of the time scales (Subscript[T, 0]), (Subscript[T, 1]),  and (Subscript[T, 2]), we have

\(\(multiScales = {u_\_i_[t] -> u\_i[T\_0, T\_1, T\_2], 
        \(\(Derivative[n_]\)[u_\_i_]\)[t] :> 
          \(dt[n]\)[u\_i[T\_0, T\_1, T\_2]]}; \)\)

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

eq522a = 
      (EOM2 /. multiScales /. solRule // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we obtain

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

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

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

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

eqOrder[1] /. displayRule // TableFormeqOrder[2] /. displayRule // TableFormeqOrder[3] /. displayRule // TableForm

First-Order Equations

The homogeneous solution of eqOrder[1] can be obtained by using

eigenForm = 
      #[[2]] Exp[#[[1]] (Subscript[T, 0])] & /@ 
        Transpose[Eigensystem[matrixA]]

\(sol1hForm = 
      Flatten[Table[
            I\ \(ω\_i\) {\(-A\&_\_i[T\_1, T\_2]\), A\_i[T\_1, T\_2]}, {
              i, 2}]] . eigenForm\)

Since the first-order equations are homogeneous, the solution of eqOrder[1] can be expressed as

states1 = {(Subscript[u, 1, 1]), (Subscript[v, 1, 1]), (Subscript[u, 2, 1]), (Subscript[v, 2, 1])};

sol1 = 
      Table[states1[[i]] -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, sol1hForm[[i]] // Evaluate], {i, 
          4}]

Second-Order Equations

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

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

    order2Eq /. displayRule\)

Next, we substitute the expRule1[i] into the right-hand sides of order2Eq to transform the small-divisor terms into secular terms, collect the terms responsible for secular terms, and obtain

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

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

The solvability conditions of the second-order equations demands that ST11 and ST12 be orthogonal to every solution of the adjoint homogeneous problems, namely the components of ccleftVec. Imposing these conditions, we obtain

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

or

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]];
SCond1Rule1 /. displayRule

whose complex conjugate is

ccSCond1Rule1 = SCond1Rule1 /. conjugateRule;

Using these solvability conditions and expRule3, we rewrite order2Eq as

\(order2Eqm = 
      \(\(order2Eq /. SCond1Rule1\) /. ccSCond1Rule1\) /. expRule3; 

    order2Eqm /. displayRule\)

Transforming order2Eqm into a set of ordinary-differential equations and using DSolve, we obtain the particular solutions as

\(\((sol2p = 
          \(DSolve[
                  order2Eqm /. u_\_\(i_, 2\) -> \((u\_\(i, 2\)[#1]\ &)\), {
                    u\_\(1, 2\)[T\_0], v\_\(1, 2\)[T\_0], u\_\(2, 2\)[T\_0], 
                    v\_\(2, 2\)[T\_0]}, T\_0]\)[\([1]\)] /. C[_] -> 0 // 
            Simplify; )\) // Timing\)

Hence, the second-order solutions of order2Eqm can be expressed as

\(states2 = states1 /. u_\_\(i_, 1\) -> u\_\(i, 2\)\)

\(sol2 = 
      Table[states2[\([i]\)] -> 
          Function[{T\_0, T\_1, T\_2}, 
            Expand[sol2p[\([i, 2]\)]] /. Exp[a_] :> Exp[a // Expand] // 
              Evaluate], {i, 4}]; 
sol2 /. displayRule\)

Third-Order Equations

Substituting the first- and second-order solutions into the right-hand sides of the third-order equations, eqOrder[3], yields

order3Eqrhs = 
      Expand[#[[2]] /. sol1 /. sol2] & /@ eqOrder[3];

To eliminate the terms that lead to secular terms from order3Eqrhs, we use expRule1[i] to transform the terms that convert small-divisor terms into secular terms, calculate the vectors proportional to the (E^(I (Subscript[ω, i]) (Subscript[T, 0]))),  and obtain

\(ST21 = 
      Coefficient[order3Eqrhs /. expRule1[1], Exp[I\ ω\_1\ T\_0]]; 

    ST21 /. displayRule\)

\(ST22 = 
      Coefficient[order3Eqrhs /. expRule1[2], Exp[I\ ω\_2\ T\_0]]; 

    ST22 /. displayRule\)

Hence, demanding that ST21 and ST22 be orthogonal to the components of ccleftVec yields the solvability conditions

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

or

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]];
SCond2Rule1 /. displayRule

Reconstitution

We reconstitute SCond1Rule1 and SCond2Rule1 and obtain

\(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 // Expand)\) // 
        Thread; 
moduEq /. displayRule\)

These equations are derivable from a Lagrangian because

\(\(Coefficient[moduEq[\([1, 2]\)], 
              A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_2[T\_1, T\_2]] - 
            Coefficient[moduEq[\([2, 2]\)], 
              A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_1[T\_1, T\_2]] == 
          0 /. ϵ -> 1\) /. ω\_2 -> 2  ω\_1\)

Therefore, treatment of the first-order version of the governing equations yields consistent results without floating constants.

5.2.3 First-Order Complex-Valued System

In this section, we show that treatment of a first-order complex -valued version of the governing equations yields consistent results without floating constants. We start by transforming the second-order form EOM1 of the governing equations into a system of first-order complex-valued equations by introducing the transformation

\(\(transfRule = {u\_k_[t] -> ζ\_k[t] + ζ\&_\_k[t], 
        \(u\_k_\^′\)[t] -> 
          I\ ω\_k\ \((ζ\_k[t] - ζ\&_\_k[t])\)}; \)\)

Then, we define the dependent variables

depVar = {(Subscript[u, 1])[t], ((Subscript[u, 1])^′)[t], (Subscript[u, 2])[t], 
        ((Subscript[u, 2])^′)[t]};

\(\(nmodes = {ζ\_1[t], ζ\&_\_1[t], ζ\_2[t], 
        ζ\&_\_2[t]}; \)\)

Substituting transfRule into depVar and solving for the nmodes, we obtain

zetaRule = 
      Solve[depVar == (depVar /. transfRule) // Thread, nmodes][
        [1]]

It follows from EOM1 that the acceleration is given by

acceleration = 
      Solve[EOM1, {((Subscript[u, 1])^″)[t], 
            ((Subscript[u, 2])^″)[t]}][[1]]

Then, differentiating the dependent variables in zetaRule and using transfRule and acceleration, we transform the system of two second-order equations EOM1 into the following two first-order complex-valued equations:

\(EOM3 = 
      \(\(ζ\_#\^′\)[t] == 
            \((\(D[ζ\_#[t] /. zetaRule, t] /. acceleration\) /. 
                  transfRule // Expand)\)\ &\)\ /@\ {1, 2}\)

To determine a second-order uniform expansion of the solution of EOM3 using the method of multiple scales, we expand ζ and \(TraditionalForm\`ζ\&_\) in the form

\(\(solRule = {
        ζ\_i_ -> 
          \((Sum[ϵ\^j\ ζ\_\(i, j\)[#1, #2, #3], {j, 3}]\ &)\), 
        ζ\&_\_i_ -> 
          \((Sum[ϵ\^j\ ζ\&_\_\(i, j\)[#1, #2, #3], {j, 3}]\ &)
            \)}; \)\)

eq523a = 
      (EOM3 /. multiScales /. solRule // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we obtain

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

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

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

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

eqOrder[1] /. displayRule // TableFormeqOrder[2] /. displayRule // TableFormeqOrder[3] /. displayRule // TableForm

First-Order Equations

Since eqOrder[1] are homogeneous, we can write their solutions as

\(sol1 = 
      Table[{ζ\_\(i, 1\) -> 
              Function[{T\_0, T\_1, T\_2}, 
                A\_i[T\_1, T\_2]\ Exp[I\ ω\_i\ T\_0] // Evaluate], 
            ζ\&_\_\(i, 1\) -> 
              Function[{T\_0, T\_1, T\_2}, 
                A\&_\_i[T\_1, T\_2]\ Exp[\(-I\)\ ω\_i\ T\_0] // 
                  Evaluate]}, {i, 2}] // Flatten\)

Second-Order Equations

Substituting the first-order solution into the second-order equations, eqOrder[2], we have

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

Converting the terms that produce small-divisor terms by using the expRule1 and then eliminating the terms that lead to secular terms from order2Eq yields

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

    SCond1 /. displayRule\)

or

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]];
SCond1Rule1 /. displayRule

Substituting SCond1Rule1 and expRule3 into order2Eq yields

\(order2Eqm = \(order2Eq /. SCond1Rule1\) /. expRule3; 

    order2Eqm /. displayRule\)

We transform order2Eqm into a set of ordinary-differential equations and obtain their particular solutions as

sol2p = 
      DSolve[order2Eqm /. 
              (Subscript[ζ, i_, 2]) -> ((Subscript[ζ, i, 2])[#1] &), {
              (Subscript[ζ, 1, 2])[(Subscript[T, 0])], (Subscript[ζ, 2, 2])[(Subscript[T, 0])]}, (Subscript[T, 0])][
          [1]] /. C[_] -> 0;

expRule4 = Exp[a_] :> Exp[Expand[a]];

We write the second-order solution in function form as

\(sol2 = 
      Table[{ζ\_\(i, 2\) -> 
              Function[{T\_0, T\_1, T\_2}, 
                \((sol2p[\([i, 2]\)] // Expand)\) /. expRule4 // Evaluate], 
            ζ\&_\_\(i, 2\) -> 
              Function[{T\_0, T\_1, T\_2}, 
                \(\((sol2p[\([i, 2]\)] // Expand)\) /. expRule4\) /. 
                    conjugateRule // Evaluate]}, {i, 2}] // Flatten; 

    sol2 /. displayRule\)

Third-Order Equations

Substituting the first- and second-order solutions into the right-hand sides of the third-order equations, eqOrder[3], we obtain

order3Eqrhs = 
      Expand[#[[2]] /. sol1 /. sol2] & /@ eqOrder[3];

Converting the terms that produce small-divisor terms and eliminating the terms that produce secular terms in (Subscript[ζ,  i, 3]) from order3Eqrhs demands the following conditions:

\(SCond2 = 
      \(Coefficient[order3Eqrhs[\([#]\)] /. expRule1[#], 
              Exp[I\ ω\_#\ T\_0]] == 0\ &\)\ /@\ {1, 2}; 

    SCond2 /. displayRule\)

or

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]];
SCond2Rule1 /. displayRule

Reconstitution

We reconstitute SCond1Rule1 and SCond2Rule1 and obtain

\(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 // Expand)\) // 
        Thread; 
moduEq /. displayRule\)

5.3 Method of Normal Forms

To apply the method of normal forms, we find it convenient to start with the first-order complex-valued form EOM3

EOM3={(Subscript[ζ, 1]^′)[t]==\(I\ ω\_1\ ζ\_1[t] - 
                \(I\ δ\ ζ\_1[t]\ ζ\_2[t]\)\/ω\_1 - 
                \(I\ δ\ ζ\_2[t]\ 
                    ζ\&_\_1[t]\)\/ω\_1 - 
                \(I\ δ\ ζ\_1[t]\ 
                    ζ\&_\_2[t]\)\/ω\_1 - 
                \(I\ δ\ ζ\&_\_1[t]\ 
                    ζ\&_\_2[t]\)\/ω\_1\),(Subscript[ζ, 2]^′)[t]==\(\(-\(\(I\ δ\ ζ\_1[t]\^2\)\/\(2\ 
                        ω\_2\)\)\) + I\ ω\_2\ ζ\_2[t] - 
                \(I\ δ\ ζ\_1[t]\ 
                    ζ\&_\_1[t]\)\/ω\_2 - 
                \(I\ δ\ ζ\&_\_1[t]\^2\)\/\(2\ 
                    ω\_2\)\)};

of the governing equations. Then, according to the method of normal forms, we introduce the near-identity transformation

\(basicTerms = 
      Table[{η\_k[t], η\&_\_k[t]}, {k, 2}] // Flatten\)

\(zetaRule = {
        ζ\_i_ -> 
          Function[t, 
            ϵ\ η\_i[t] + 
                Sum[ϵ\^\(j + 1\)\ 
                    h\_\(i, j\)[Sequence\ @@\ basicTerms], {j, 2}] // 
              Evaluate], 
        ζ\&_\_i_ -> 
          Function[t, 
            ϵ\ η\&_\_i[t] + 
                Sum[ϵ\^\(j + 1\)\ 
                    h\&_\_\(i, j\)[Sequence\ @@\ basicTerms], {j, 2}] // 
              Evaluate]}\)

that results in the simplest possible equations

\(etaRule = 
      \(η\_i_\^′\)[t] -> 
        I\ ω\_i\ η\_i[t] + 
          Sum[ϵ\^j\ g\_\(i, j\)[t], {j, 2}]\)

where the overbar denotes the complex conjugate

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

Substituting the expansion for the (Subscript[ζ,  i]), the zetaRule, into EOM3, using the etaRule, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3),  we have

eq53a = 
      (EOM3 /. zetaRule /. etaRule /. 
              (etaRule /. conjugateRule) // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0;

Second-Order Equations

Collecting the coefficients of (ϵ^2) in eq53a yields

order2expr = 
      Coefficient[Subtract @@ #, (ϵ^2)] & /@ eq53a;

To determine the resonance terms in order2expr, we note that all of the possible forms of its nonhomogeneous terms are

possibleQTerms = 
      Outer[Times, basicTerms, basicTerms] // Flatten // Union

It follows from the linear parts of the etaRule that to the first approximation

\(\(form = {η\_i_[t] -> E\^\(I\ ω\_i\ t\), 
        η\&_\_i_[t] -> E\^\(\(-I\)\ ω\_i\ t\)}; \)\)

Hence, the possible resonance terms are given by

ResonantQTerm[1] = 
      ((E^(-I (Subscript[ω, 1]) t)) possibleQTerms /. form /. 
                  (Subscript[ω, 2]) -> 2 (Subscript[ω, 1]) /. (E^(_ t)) -> 0) 
            possibleQTerms // Union // Rest

ResonantQTerm[2] = 
      ((E^(-I (Subscript[ω, 2]) t)) possibleQTerms /. form /. 
                  (Subscript[ω, 1]) -> (1/2) (Subscript[ω, 2]) /. (E^(_ t)) -> 0) 
            possibleQTerms // Union // Rest

Then, the nonresonance terms and their associated coefficients can be defined as

Do[NRQT[i] = Complement[possibleQTerms, ResonantQTerm[i]]; 
      coeffsQ[i] = Table[(Subscript[Γ, i, j]), {j, Length[NRQT[i]]}], {
        i, 2}]

Thus, the (Subscript[h, i, 1]) and their complex conjugates have the form

\(\(hFormQ = 
      Table[{h\_\(i, 1\) -> 
              \((Evaluate[
                    coeffsQ[i] . NRQT[i] /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\), 
            h\&_\_\(i, 1\) -> 
              \((Evaluate[
                    \(coeffsQ[i] . NRQT[i] /. conjugateRule\) /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\)}, {i, 
            2}] // Flatten; \)\)

Substituting for the (Subscript[h, i, 1]) in order2expr, equating the coefficients of the possible nonresonance terms to zero, and solving the resulting algebraic equations for the (Subscript[Γ, i,  j]), we obtain

coeffsQRule = 
      Table[Solve[
              Coefficient[order2expr[[i]] /. hFormQ, NRQT[i]] == 0 // 
                Thread, coeffsQ[i]][[1]], {i, 2}] // Flatten

Using these coefficients, we rewrite hFormQ as

hSolQ = 
      Table[(Subscript[h, i, 1])[Sequence @@ basicTerms] -> 
          ((Subscript[h, i, 1])[Sequence @@ basicTerms] /. hFormQ /. 
              coeffsQRule), {i, 2}]

or in pure function form as

\(\(hRuleQ = 
      Table[{h\_\(i, 1\) -> 
              \((Evaluate[
                    hSolQ[\([i, 2]\)] /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\), 
            h\&_\_\(i, 1\) -> 
              \((Evaluate[
                    \(hSolQ[\([i, 2]\)] /. conjugateRule\) /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\)}, {i, 
            2}] // Flatten; \)\)

We choose the (Subscript[g, i, 1]) to eliminate the resonance terms in order2expr according to

gRuleQ = 
      Table[(Subscript[g, i, 1])[t] -> 
          - Coefficient[order2expr[[i]], ResonantQTerm[i]] . 
              ResonantQTerm[i], {i, 2}]

Third-Order Equations

Collecting the coefficients of (ϵ^3) in eq53a and using hRuleQ and gRuleQ, we have

order3expr = 
      Map[Coefficient[Subtract @@ #, (ϵ^3)] &, eq53a] /. hRuleQ
               /. gRuleQ /. (gRuleQ /. conjugateRule) // Expand;

The nonhomogeneous terms in order3expr are proportional to

possibleCTerms = 
      Outer[Times, possibleQTerms, basicTerms] // Flatten // Union

Next, we determine the resonance terms according to

ResonantCTerm[1] = 
      ((E^(-I (Subscript[ω, 1]) t)) possibleCTerms /. form /. 
                  (Subscript[ω, 2]) -> 2 (Subscript[ω, 1]) /. (E^(_ t)) -> 0) 
            possibleCTerms // Union // Rest

ResonantCTerm[2] = 
      ((E^(-I (Subscript[ω, 2]) t)) possibleCTerms /. form /. 
                  (Subscript[ω, 1]) -> (1/2) (Subscript[ω, 2]) /. (E^(_ t)) -> 0) 
            possibleCTerms // Union // Rest

Then, the nonresonance terms and their associated coefficients are given by

Do[NRCT[i] = Complement[possibleCTerms, ResonantCTerm[i]]; 
      coeffsC[i] = Table[(Subscript[Λ, i, j]), {j, Length[NRCT[i]]}], {
        i, 2}]

Consequently, the (Subscript[h, i, 2]) and their complex conjugates have the forms

\(\(hFormC = 
      Table[{h\_\(i, 2\) -> 
              \((Evaluate[
                    coeffsC[i] . NRCT[i] /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\), 
            h\&_\_\(i, 2\) -> 
              \((Evaluate[
                    \(coeffsC[i] . NRCT[i] /. conjugateRule\) /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\)}, {i, 
            2}] // Flatten; \)\)

Substituting for the (Subscript[h, i, 2]) in order3expr, equating the coefficient of each possible nonresonance term to zero, and solving the resulting algebraic equations for the (Subscript[Λ, i,  j]), we obtain

coeffsCRule = 
      Table[Solve[
              Coefficient[order3expr[[i]] /. hFormC, NRCT[i]] == 0 // 
                Thread, coeffsC[i]][[1]], {i, 2}] // Flatten

Using these coefficients, we rewrite hFormC as

hSolC = 
      Table[(Subscript[h, i, 2])[Sequence @@ basicTerms] -> 
          ((Subscript[h, i, 2])[Sequence @@ basicTerms] /. hFormC /. 
              coeffsCRule), {i, 2}]

or in pure function form as

\(\(hRuleC = 
      Table[{h\_\(i, 2\) -> 
              \((Evaluate[
                    hSolC[\([i, 2]\)] /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\), 
            h\&_\_\(i, 2\) -> 
              \((Evaluate[
                    \(hSolC[\([i, 2]\)] /. conjugateRule\) /. 
                      Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\)}, {i, 
            2}] // Flatten; \)\)

Choosing the (Subscript[g, i, 2]) to eliminate the resonance terms from order3expr yields

gRuleC = 
      Table[(Subscript[g, i, 2])[t] -> 
          - Coefficient[order3expr[[i]], ResonantCTerm[i]] . 
              ResonantCTerm[i], {i, 2}]

Combining etaRule, gRuleQ, and gRuleC and letting (Subscript[η, i]) = (Subscript[A, i])[t] (E^(I (Subscript[ω, i]) t)), we obtain the modulation equation

\(moduEq = 
      Table[\(2\ I\ ω\_k\ \(E\^\(\(-I\)\ ω\_k\ t\)\) 
                  \((\(η\_k\^′\)[t] - 
                      \((\(\(\(η\_k\^′\)[t] /. etaRule\) /. gRuleQ
                            \) /. gRuleC)\))\) == 0 /. {
                η\_i_ -> \((A\_i[#]\ E\^\(I\ ω\_i\ #\)&)\), 
                η\&_\_i_ -> 
                  \((A\&_\_i[#]\ E\^\(\(-I\)\ ω\_i\ #\)&)\)}\) /. 
            Exp[a_] :> 
              Exp[a /. 
                    ω\_2 -> 
                      2  ω\_1 + \(ϵ\^2\) σ // 
                  Expand] // ExpandAll, {k, 2}]\)

Using zetaRule, hRuleQ, and hRuleC, we obtain the following second-order uniform expansion of the solution of EOM3:

\(solution = 
      Table[u\_k[t] == 
          Collect[\(\(ζ\_k[t] + ζ\&_\_k[t] /. zetaRule\) /. hRuleQ
                \) /. hRuleC, ϵ], {k, 2}]\)

5.4 Generalized Method of Averaging

We define the nonlinear part of EOM1, a transformation rule, and the states vector as

\(NLTerms = {f\_1 -> 2\ ϵ\ δ\ u\_1[t]\ u\_2[t], 
        f\_2 -> ϵ\ δ\ u\_1[t]\^2}; 

    transformRule = u\_i_ -> \((a\_i[#]\ Sin[ϕ\_i[#]]\ &)\); 

    states = {a\_1[t], a\_2[t], ϕ\_1[t], ϕ\_2[t]}; \)

Substituting the above rules into the four first-order equations obtained by using the method of variation of parameters, we have

eq54a = 
      D[states, t] == 
          ({((Subscript[f, 1])/(Subscript[ω, 1])) Cos[(Subscript[ϕ, 1])[t]], 
                      ((Subscript[f, 2])/(Subscript[ω, 2])) Cos[(Subscript[ϕ, 2])[t]], 
                      (Subscript[ω, 1]) - 
                        ((Subscript[f, 1])/((Subscript[a, 1])[t] (Subscript[ω, 1]))) Sin[(Subscript[ϕ, 1])[t]], 
                      (Subscript[ω, 2]) - 
                        ((Subscript[f, 2])/((Subscript[a, 2])[t] (Subscript[ω, 2]))) Sin[(Subscript[ϕ, 2])[t]]} /. 
                    NLTerms /. transformRule // TrigReduce // Expand)
             // Thread

We seek a second-order approximate solution of eq54a in the form

basicTerms = {(Subscript[η, 1])[t], (Subscript[η, 2])[t], (Subscript[φ, 1])[t], 
        (Subscript[φ, 2])[t]};

solRule = { 
        (Subscript[a, i])_ -> 
          (Evaluate[
                (Subscript[η, i])[t] + 
                    Sum[(ϵ^j) 
                        (Subscript[a, i, j])[Sequence @@ basicTerms], {j, 2}] /. 
                  t -> #] &), 
        (Subscript[ϕ, i])_ -> 
          (Evaluate[
                (Subscript[φ, i])[t] + 
                    Sum[(ϵ^j) 
                        (Subscript[ϕ, i, j])[Sequence @@ basicTerms], {j, 
                        2}] /. t -> #] &)};

where the (Subscript[η, i])[t] and (Subscript[φ, i])[t] are expanded in power series in ϵ as

basicDRule = 
      D[basicTerms, t] -> {
            Sum[(ϵ^i) (Subscript[A, 1, i])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], {i, 
                2}], Sum[
              (ϵ^i) (Subscript[A, 2, i])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], {i, 
                2}], (Subscript[ω, 1]) + 
              Sum[(ϵ^i) 
                  (Subscript[Φ, 1, i])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], {i, 
                  2}], (Subscript[ω, 2]) + 
              Sum[(ϵ^i) 
                  (Subscript[Φ, 2, i])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], {i, 
                  2}]} // Thread

The functions (Subscript[a, i,  j]) and (Subscript[ϕ, i,  j]) are fast varying functions of the (Subscript[φ, i]), while it follows from basicDRule that the (Subscript[η, i]), and hence the (Subscript[A, i,  j]) and (Subscript[Φ, i,  j]), are slowly varying functions of t.

To the second approximation, we differentiate the (Subscript[a, i])[t] and (Subscript[ϕ, i])[t] once with respect to t, use solRule and basicDRule, expand the result for small ϵ, discard terms of order higher than (ϵ^2), and obtain

eq54bLHS = 
      CoefficientList[
            Expand[#[[1]] /. solRule /. basicDRule] /. 
              (ϵ^(n_ /; n > 2)) -> 0, ϵ] & /@ eq54a;

Next, we substitute solRule into the right-hand sides of eq54a, expand the result for small ϵ, keep terms up to O((ϵ^2)), and rewrite their right-hand sides as

eq54bRHS = 
      CoefficientList[
            Series[#[[2]] /. solRule, {ϵ, 0, 2}] // Normal, 
            ϵ] & /@ eq54a;

Equating coefficients of like powers of ϵ in eq54bLHS and eq54bRHS, we obtain

eqEps = 
      MapThread[Equal, {eq54bLHS, eq54bRHS}, 2] // Transpose // Rest;

To express the nearness of the two-to-one internal resonance, we introduce the detuning parameter σ defined according to

omgRule = {(Subscript[ω, 2]) -> 2  (Subscript[ω, 1]) + ϵ σ};

Hence, σ can be related to the (Subscript[ω, i]) by

sigRule = Solve[omgRule /. Rule -> Equal, σ][[1]]

To use the method of averaging, we define the following rules:

psiRule = {(Subscript[φ, 1])[t] -> (Subscript[ω, 1]) t + (Subscript[β, 1]), 
        (Subscript[φ, 2])[t] -> (Subscript[ω, 2]) t + (Subscript[β, 2])};

betaRule = 
      Solve[psiRule /. Rule -> Equal, {(Subscript[β, 1]), (Subscript[β, 2])}][[1]]

expRule1 = 
      f_[a_] :> 
        f[Expand[a /. psiRule /. omgRule] /. ϵ t -> (Subscript[t, 1])];

expRule2 = 
      f_[a_] :> 
        f[Expand[a /. (Subscript[t, 1]) -> ϵ t /. sigRule /. betaRule]];

Next, we separate the fast and slowly varying terms in the first-order equations, eqEps[[1]]. To this end, we define a slow state vector and a fast state vector according to

SVT[j_] := {(Subscript[A, 1, j])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], 
        (Subscript[A, 2, j])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], 
        (Subscript[Φ, 1, j])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]], 
        (Subscript[Φ, 2, j])[(Subscript[η, 1])[t], (Subscript[η, 2])[t]]}solVar[j_] := {(Subscript[a, 1, j]), (Subscript[a, 2, j]), (Subscript[ϕ, 1, j]), 
        (Subscript[ϕ, 2, j])}

Then, the slowly varying parts of eqEps[[1]] are given by

SVT1Rule = 
      Table[Solve[
              eqEps[[1, i]] /. 
                        Thread[solVar[1] -> (0 &)] /. expRule1 /. 
                    f_[_ t + _. ] -> 0 /. expRule2 /. 
                f_[a_] :> f[Expand[a]], SVT[1][[i]]] [[1]], {
            i, 4}] // Flatten

whereas the fast varying parts are given by

FVT1 = 
      Table[Subtract @@ eqEps[[1, i]] /. 
            Thread[solVar[1] -> (0 &)], {i, 4}] /. SVT1Rule

To determine the solution corresponding to these fast varying terms, we use the method of undetermined coefficients. To accomplish this, we first identify the possible forms of the nonhomogeneous terms. The result is

FVT1Forms = 
      Flatten[Cases[#, (Cos | Sin)[a_] -> {Cos[a], Sin[a]}, Infinity]]
             & /@ FVT1

The fast-varying component of the first-order solution can be expressed as a linear combination of these forms; that is,

sol1Form = 
      MapIndexed[(coeffs1[#2[[1]]] = Array[c, Length[#1]]) . #1 &, 
        FVT1Forms]

Substituting sol1Form into eqEps[[1]], using SVT1Rule, collecting the coefficients of FVT1Forms, solving the resulting algebraic equation for the undetermined coefficients, and then substituting the result back into sol1Form, we obtain

sol1rhs = 
      Table[sol1Form[[i]] /. 
          Solve[Coefficient[
                    Subtract @@ eqEps[[1, i]] /. SVT1Rule /. 
                      solVar[1][[i]] -> 
                        (Evaluate[
                              sol1Form[[i]] /. 
                                Thread[basicTerms -> {#1, #2, #3, #4}]] &)
                          , FVT1Forms[[i]]] == 0 // Thread, coeffs1[i]]
              [[1]], {i, 4}]

or in a pure function form as

sol1Rule = 
      Table[solVar[1][[i]] -> 
          (Evaluate[
                sol1rhs[[i]] /. Thread[basicTerms -> {#1, #2, #3, #4}]]
                 &), {i, 4}]

Substituting the slow- and fast-varying components of the first-order solution, SVT1Rule and sol1Rule, into the second-order equations, we have

(order2Expr = 
          Table[Subtract @@ eqEps[[2, i]] /. SVT1Rule /. 
                  sol1Rule // Expand // TrigReduce, {i, 4}]; ) // Timing

Since we are seeking an expansion valid up to O((ϵ^2)), we do not need to solve for the (Subscript[a, i,  2]) and (Subscript[ϕ, i,  2]). All we need to do is to investigate the above expressions to determine the slowly varying parts and hence determine the (Subscript[A, i,  2]) and (Subscript[Φ, i,  2]). The result is

SVT2Rule = 
      Table[Solve[
            order2Expr[[i]] == 0 /. 
                      Thread[solVar[2] -> (0 &)] /. expRule1 /. 
                  f_[_ t + _. ] -> 0 /. expRule2 /. 
              f_[a_] :> f[Expand[a]], SVT[2][[i]]], {i, 4}] // Flatten

Hence, we find that, to the second approximation, the solution is given by

eq54c = {(Subscript[u, 1])[t], (Subscript[u, 2])[t]} == 
          (Normal[
                    Series[{(Subscript[u, 1])[t], (Subscript[u, 2])[t]} /. transformRule /. 
                        solRule, {ϵ, 0, 1}]] /. sol1Rule // 
                TrigReduce // Expand) // Thread

where the amplitudes and phases are given by

ampEqs = 
      basicDRule[[{1, 2}]] /. SVT1Rule /. SVT2Rule /. psiRule
             /. f_[a_] :> f[Collect[a, t]] /. Rule -> Equal

phaseEqs = 
      basicDRule[[{3, 4}]] /. SVT1Rule /. SVT2Rule /. psiRule
             /. f_[a_] :> f[Collect[a, t]] /. Rule -> Equal