Chapter 3
Systems with Quadratic and Cubic Nonlinearities

Off[General::spell1]

3.1 Nondimensional Equation of Motion

We consider the free oscillations of a particle of mass m under the action of gravity and restrained by a nonlinear spring. The equation of motion is

\(eq31a = 
      m\ D[\(x\^*\)[\(t\^*\)], {\(t\^*\), 2}] + f[\(x\^*\)[\(t\^*\)]] == 
        m\ g\)

where g is the gravitational acceleration and \(TraditionalForm\`f[\(x\^*\)]\) is the restraining force of the spring. We assume that \(TraditionalForm\`f[\(x\^*\)]\) is an odd cubic function of \(TraditionalForm\`\(\(x\^*\); \)\) that is,

eq31b = eq31a /. f[a_] -> (Subscript[k, 1]) a + (Subscript[k, 3]) (a^3)

The equilibrium positions \(TraditionalForm\`\(x\^*\) = \(\(x\_s\%*\) = constant\)\) can be obtained by dropping the acceleration term. The result is

\(eq31c = eq31b /. \(x\^*\) -> \((\(x\_s\%*\)\ &)\)\)

In this chapter, we investigate small oscillations about one of the equilibrium positions. To this end, we shift this equilibrium position to the origin by using the transformation

\(\(xRule = {\(x\^*\) -> \((\(x\_s\%*\) + \(u\^*\)[#]\ &)\)}; \)\)

Substituting this transformation into eq31b and using the equilibrium condition eq31c, we obtain

eq31d = 
      Expand[Subtract @@ eq31b /. xRule] == 0 /. 
        (eq31c /. Equal -> Rule)

As before, we introduce the following dimensionless quantities:

\(\(dimenRule = {\(u\^*\) -> \((\(x\_s\%*\)\ u[#/T]\ &)\), 
        \(t\^*\) -> T\ t}; \)\)

into eq31d and obtain

\(eq31e = 
      \((eq31d[\([1]\)]\ T\^2/\((m\ \(x\_s\%*\))\) /. dimenRule // Expand)
          \) == 0\)

We rewrite eq31e as

eq31f = 
      (Collect[eq31e[[1]], 
              u[t]] /. {_ u[t]^n_.  -> (Subscript[α, n]) u[t]^n}) == 0

where the (Subscript[α, j]) are given by

\(alphas = 
      \(\((α\_# -> Coefficient[eq31e[\([1]\)], u[t]\^#])\)\ &\)\ /@
        \ {1, 2, 3}\)

We choose T so that (Subscript[α, 1]) = 1 and hence the natural frequency is unity; that is,

TRule = Solve[(Subscript[α, 1]) == 1 /. alphas, T][[2]]

Consequently, (Subscript[α, 2]) and (Subscript[α, 3]) can be rewritten as

alphas = alphas /. TRule // Rest

With (Subscript[α, 1]) being unity, eq31f becomes

eq31g = eq31f /. (Subscript[α, 1]) -> 1

In contrast with the DuffingEq, eq31g contains a quadratic as well as a cubic term.

In the next section, we determine a second-order straightforward expansion to the solutions of eq31g for small but finite amplitudes. In Section 3.3, we determine a uniform second-order expansion by using the Lindstedt-Poincaré technique. In Section 3.4, we determine a uniform second-order expansion by using the method of multiple scales. In Section 3.5, we show that the first approximation obtained with the method of averaging yields an incomplete solution. In Section 3.6, we introduce the generalized method of averaging and obtain a uniform second-order expansion. Finally, in Sections 3.7 and 3.8, we introduce the Krylov-Bogoliubov-Mitropolsky technique and the method of normal forms, respectively.

3.2 Straightforward Expansion

To carry out a straightforward expansion for small but finite amplitudes for eq31g, we need to introduce a small parameter ϵ as a bookkeeping parameter. In terms of this parameter, we seek a third-order expansion in the form

uExpRule[maxOrder_] := 
      u -> (Sum[(ϵ^i) (Subscript[u, i])[#], {i, 1, maxOrder}] &)

where maxOrder is the order of the expansion sought. In what follows, we let maxOrder = 3 for an expansion of order three. Clearly, ϵ is a measure of the amplitude of oscillation and can be set equal to unity in the final solution if the amplitude is taken to be small, as described below. Substituting this expansion into eq31g, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we have

eq32a = 
      (eq31g /. uExpRule[3] // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0

Equating coefficients of like powers of ϵ in eq32a yields

eqEps = 
      Table[Coefficient[eq32a[[1]], (ϵ^i)] == 0, {i, 3}]

which can be solved sequentially for (Subscript[u, 1]), (Subscript[u, 2]),  and (Subscript[u, 3]).

The general solution of the first-order equation, eqEps[[1]], can be expressed as

sol1 = 
      DSolve[{eqEps[[1]], (Subscript[u, 1])[0] == a Cos[β], 
              ((Subscript[u, 1])^′)[0] == -a Sin[β]}, (Subscript[u, 1])[t], t][
          [1]] // Simplify

where a and β are constants. Substituting sol1, (Subscript[u, 1])[t], into the second-order equation, eqEps[[2]], and solving the resulting equation for (Subscript[u, 2])[t], we obtain

sol2 = 
      DSolve[eqEps[[2]] /. sol1, (Subscript[u, 2])[t], t][[1]] /. C[_] -> 0 // 
        Simplify

As before, we do not include the solution of the homogeneous problem for (Subscript[u, 2])[t].

Substituting sol1 and sol2 into the third-order equation, eqEps[[3]], and solving the resulting equation (Subscript[u, 3])[t], we have

sol3 = 
      DSolve[eqEps[[3]] /. sol1 /. sol2, (Subscript[u, 3])[t], t][[1]] /. 
          C[_] -> 0 // Simplify

Simply letting C[_] → 0 did not remove all of the solutions of the homogeneous equation in this case. Therefore, we clear it one more time and obtain

sol3 = sol3 /. Cos[t + _] -> 0 // ExpandAll

Combining the first-, second-, and third-order solutions, we obtain, to the third approximation, that

sol = 
      u[t] == (u[t] /. uExpRule[3] /. Flatten[{sol1, sol2, sol3}])

We note that the dependence of u on ϵ and a is in the combination ϵ a. Thus, one can set ϵ = 1 and consider a as the perturbation parameter.

The straightforward expansion sol breaks down for \(TraditionalForm
      \`t ≥ O(\(ϵ\^\(-1\)\) a\^\(-1\))\) because the third-order term is the same order or larger than the second-order term owing to the presence of the mixed-secular term. Next, we use the Lindstedt-Poincaré technique to determine a third-order uniform solution.

3.3 The Lindstedt-Poincaré Technique

To apply the Lindstedt-Poincaré technique to eq31g, we introduce the transformation τ = ω t and rewrite it as

eq33a = 
      eq31g /. u -> (u[ω #] &) /. t -> τ/ω

To determine a uniform expansion of order maxOrder of the solution of eq33a, we expand both ω and u[τ] in powers of ϵ as

omgRule[maxOrder_] := 
      ω -> 1 + Sum[(ϵ^i) (Subscript[ω, i]), {i, maxOrder}]uExpRule[maxOrder_] := 
      u -> (Sum[(ϵ^i) (Subscript[u, i])[#], {i, maxOrder + 1}] &)

As discussed in the preceding chapter, the first term in the expansion of ω is taken to be the linear natural frequency of the system, which is unity in this case. We note that the order of the expansion for u is larger than that of ω and the uniform expansion required. As discussed below, we eliminate the terms that lead to secular terms from the equation at order maxOrder+1 but do not include in the solution of the resulting equation in the expansion.

For a uniform second order expansion, we let

maxOrder = 2;

Substituting the expansions of u and ω into eq33a, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq33b = 
      (eq33a /. {uExpRule[maxOrder], omgRule[maxOrder]} // ExpandAll) /. 
        (ϵ^(n_ /; n > maxOrder + 1)) -> 0

eqEps = 
      Table[Coefficient[eq33b[[1]], (ϵ^i)] == 0, {i, 3}]

The general solution of the first-order equation, eqEps[[1]], can be expressed as

sol1 = 
      DSolve[{eqEps[[1]], (Subscript[u, 1])[0] == a Cos[β], 
            ((Subscript[u, 1])^′)[0] == -a Sin[β]}, (Subscript[u, 1]), τ][
        [1]]

where a and β are constants. Instead of using (Subscript[u, 1])[τ], we use (Subscript[u, 1]) in the second argument to DSolve so that sol1 can be used in the terms involved with derivatives.

Substituting sol1, (Subscript[u, 1])[τ], into the second-order equation, eqEps[[2]], moving the nonhomogeneous terms to the right-hand side of the resulting equation, we obtain

eq33c=(Subscript[u, 2])[τ]+(Subscript[u, 2]^″)[τ]-# & /@ Reverse[eqEps[[2]] /. sol1 // Simplify]

Expanding the right-hand side of eq33c in a Fourier series using trigonometric identities yields

eq33c[[2]] = eq33c[[2]] // TrigReduce

Eliminating the terms, Cos[β + τ] and Sin[β + τ], demands that (Subscript[ω, 1]) = 0. Then, the particular solution of eq33c can be expressed as

sol2 = 
      DSolve[eq33c /. (Subscript[ω, 1]) -> 0, (Subscript[u, 2])[τ], τ][[1]] /. 
          C[_] -> 0 // Simplify

or in a pure function format as

sol21 = (Subscript[u, 2]) -> Function[τ, sol2[[1, 2]] // Evaluate]

Substituting sol1 and sol21 into the third-order equation, eqEps[[3]], and using the fact that (Subscript[ω, 1]) = 0, we obtain

eq33d=(Subscript[u, 3])[τ]+(Subscript[u, 3]^″)[τ]-# & /@ Reverse[
          eqEps[[3]] /. (Subscript[ω, 1]) -> 0 /. sol1 /. sol21 // 
            Simplify]

Expanding the right-hand side of eq33d in a Fourier series using trigonometric identities, we have

eq33d[[2]] = eq33d[[2]] // TrigReduce

Eliminating the terms that lead to secular terms from eq33d[[2]] demands that

omg2Rule = 
      Solve[Coefficient[eq33d[[2]], Cos[τ + β]] == 0, 
            (Subscript[ω, 2])][[1]] // ExpandAll

As discussed above, for a second-order uniform expansion, we do not need to solve for (Subscript[u, 3])[τ]. Combining the first- and second-order solutions, we obtain, to the second approximation, that

u[τ] == 
      (u[τ] /. uExpRule[2] /. sol1 /. sol21 // Simplify // 
          Expand)

where

τ == ω t /. omgRule[2] /. omg2Rule /. 
      (Subscript[ω, 1]) -> 0

The above expansion is uniform to second order because secular terms do not appear in it and the correction term (the term proportional to (ϵ^2)) is small compared with the first term.

Returning to sol in the previous section, we note that the first secular term appears at O((ϵ^3)). Consequently, we could have concluded that (Subscript[ω, 1]) = 0 before carrying out the expansion because the term ϵ (Subscript[ω, 1]) in omgRule[2] creates secular terms at O((ϵ^2)) and not at O((ϵ^3)), as needed to eliminate the secular term from sol.

3.4 The Method of Multiple Scales

We use the method of multiple scales to attack directly eq31g in Section 3.4.1, the corresponding first-order real-valued equations of eq31g in Section 3.4.2, and the corresponding first-order complex-valued equation of eq31g in Section 3.4.3. To obtain a second-order uniform expansion by using the method of multiple scales, 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])];

and list as

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

For a uniform second-order expansion, we let

maxOrder = 2;

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

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

3.4.1 Second-Order Real-Valued System

Using the derivative rule, we transform eq31g into the partial-differential equation

eq341a = 
      eq31g /. {u[t] -> u[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
          Derivative[n_][u][t] -> dt[n][u[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]]}

Again, to determine a uniform expansion of order maxOrder, we expand u to order maxOrder+1, eliminate the secular terms from the maxOrder+1 equation, but do not include its solution in the final approximate solution. Hence, we seek a uniform expansion of the solution of eq341a in the form

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

Substituting this expansion into eq341a, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we have

eq341b = 
      (eq341a /. solRule // ExpandAll) /. 
        (ϵ^(n_ /; n > maxOrder + 1)) -> 0

Equating coefficients of like powers of ϵ in eq341b, we obtain

eqEps = 
      CoefficientList[eq341b[[1]], ϵ] == 0 // Thread // Rest

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

eqOrder[i_] := 
      (eqEps[[1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, i])) - # & /@  
          eqEps[[i]] // Reverse

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

Array[eqOrder, 3] /. displayRule // TableForm

where \(TraditionalForm
      \`D\_n\ u\_i = ∂u\_i[T\_0, T\_1, T\_2]/∂T\_n\).

The general solution of the first-order equation, eqOrder[1], can be expressed as

\(\(sol1 = {
        u\_1 -> \((A[#2, #3]\ Exp[I\ #1] + A\&_[#2, #3]\ Exp[\(-I\)\ #1]\ &)
            \)}; \)\)

where \(TraditionalForm\`A\&_\) is the complex conjugate of A defined by the conjugate rule

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

Then, the second-order equation, eqOrder[2], becomes

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

Eliminating the term (E^(I (Subscript[T, 0]))) that produces secular terms in (Subscript[u, 2]) from the right-hand side of eq341c demands that

SCond1 = Coefficient[eq341c[[2]], Exp[I (Subscript[T, 0])]] == 0

or

SCond1Rule=Solve[SCond1,Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]
(* either copy and paste Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])] from SCond1 or use Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])] as the second argument to Solve *)

Eliminating the term (E^(-I (Subscript[T, 0]))) that produces secular terms in (Subscript[u, 2]) from the right-hand side of eq341c yields the complex conjugate of SCond1Rule; that is,

ccSCond1Rule = SCond1Rule /. conjugateRule

It follows from SCond1 that A = A[(Subscript[T, 2])]. Substituting the solvability conditions into eq341c, we have

eq341d = eq341c /. SCond1Rule /. ccSCond1Rule

The particular solution of eq341d can be obtained by using DSolve; the result is

u2Sol = 
      DSolve[eq341d, (Subscript[u, 2])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], timeScales][[1]] /. 
            C[_] -> (0 &) // TrigToExp // ExpandAll

whose right-hand side can be used directly to express (Subscript[u, 2]) in a pure function form as

sol2 = {
        (Subscript[u, 2]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, u2Sol[[1, 2]] // Evaluate]}

Substituting the first- and second-order solutions, sol1 and sol2, into the third-order equation, eqOrder[3], we obtain

\(eq341e = \(eqOrder[3] /. sol1\) /. sol2 // ExpandAll; 

    eq341e /. displayRule\)

Eliminating the terms that produce secular terms in (Subscript[u, 3]) from the right-hand side of eq341e demands that

SCond2 = Coefficient[eq341e[[2]], Exp[I (Subscript[T, 0])]] == 0

Using SCond1Rule and the fact that A = A[(Subscript[T, 2])], we can rewrite SCond2 as

\(SCond = 
      \(SCond2 /. D[SCond1Rule, T\_1]\) /. {A -> \((A[#2]\ &)\), 
          A\&_ -> \((A\&_[#2]\ &)\)}\)

Expressing A in the polar form

\(\(ruleA = {A -> \((\(1\/2\) a[#]\ Exp[I\ β[#]]\ &)\), 
        A\&_ -> \((\(1\/2\) a[#]\ Exp[\(-I\)\ β[#]]\ &)\)}; \)\)

where a and β are real and using the rule

realRule = {Re[s_] -> s, Im[s_] -> 0};

we write SCond as

eq341f = 
      SCond[[1]] Exp[-I β[(Subscript[T, 2])]] /. ruleA // Expand

Separating the real and imaginary parts in eq341f, we obtain the modulation equations

ampEq = Im[eq341f] == 0 /. realRule

phaseEq = Re[eq341f] == 0 /. realRule

The solution of ampEq is a = (Subscript[a, 0]) = const. Then, if (Subscript[a, 0]) ≠ 0, the solution of phaseEq is

betaRule = 
      DSolve[phaseEq /. a[(Subscript[T, 2])] -> (Subscript[a, 0]), β[(Subscript[T, 2])], (Subscript[T, 2])][[1]] /. 
          C[1] -> (Subscript[β, 0]) // ExpandAll

where (Subscript[β, 0]) is a constant. Substituting ruleA into sol1 and sol2 and recalling that (Subscript[T, 0]) = t and (Subscript[T, 2]) = (ϵ^2) t, we obtain, to the second approximation, that

\(u[t] == 
      \((\(\(\(\((
                  \(\(\(\(\(u[T\_0, T\_1, T\_2] /. solRule\) /. sol1\) /. sol2
                            \) /. ϵ\^3 -> 0\) /. {A[__] -> A[T\_2], 
                          A\&_[__] -> A\&_[T\_2]}\) /. ruleA // ExpToTrig)
                  \) /. {a[T\_2] -> a\_0}\) /. betaRule\) /. {T\_0 -> t, 
              T\_2 -> ϵ\^2\ t}\) /. 
          Cos[arg_] :> Cos[Collect[arg, t]])\)\)

which is in full agreement with that obtained by using the Lindstedt-Poincaré technique.

3.4.2 First-Order Real-Valued System

In this section, we first transform eq31g into a system of two real-valued first-order equations using the transformation v[t] = (u^′)[t] and obtain

eq342a = {(u^′)[t] == v[t], 
        eq31g /. (u^″)[t] -> (v^′)[t]}

To determine a second-order uniform expansion of the solution of eq342a using the method of multiple scales, we first introduce

multiScales={u[t] -> u[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])],v[t] -> v[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])],u'[t]->dt[1][u[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]],v'[t]->dt[1][v[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]]};

to transform the derivative with respect to t in terms of the derivatives with respect to the three time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]) and then expand u and v in the form

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

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

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

Equating coefficients of like powers of ϵ in eq342b yields

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

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

\(eqOrder[i_] := 
      \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. u_\_1 -> u\_i)\) == 
          \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. u_\_1 -> u\_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 determine the solution of the first-order equations, eqOrder[1], we list their left-hand sides and obtain

linearSys = #[[1]] & /@ eqOrder[1]

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

coefList = 
      (E^(-I (Subscript[T, 0]))) linearSys /. {(Subscript[u, 1]) -> (P (E^(I #))&), 
            (Subscript[v, 1]) -> (Q (E^(I #))&)} // Expand

The coefficient matrix of coefList is

coefMat = Outer[Coefficient, coefList, {P, Q}]

and its adjoint is defined by

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

Hence, the right and left eigenvectors of coefMat are

rightVec = {1, c1} /. 
        Solve[(coefMat . {1, c1})[[1]] == 0, c1][[1]]

leftVec = {1, c1} /. 
        Solve[(hermitian[coefMat] . {1, c1})[[1]] == 0, c1][
          [1]]

Then the solution of the first-order equations can be expressed as

sol1Form = rightVec A[(Subscript[T, 1]), (Subscript[T, 2])] (E^(I (Subscript[T, 0])))

and hence (Subscript[u, 1]) and (Subscript[v, 1]) can be expressed in a pure function form as

order1Sol = {
        (Subscript[u, 1]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol1Form[[1]] + (sol1Form[[1]] /. conjugateRule) // 
              Evaluate], 
        (Subscript[v, 1]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol1Form[[2]] + (sol1Form[[2]] /. conjugateRule) // 
              Evaluate]}

Second-Order Equations

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

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

    order2Eq /. displayRule\)

To eliminate the terms that produce secular terms in (Subscript[u, 2]) and (Subscript[v, 2]), we determine the vector proportional to (E^(I (Subscript[T, 0]))) in the right-hand sides of order2Eq and obtain

STerms1 = 
      Coefficient[#[[2]], Exp[I (Subscript[T, 0])]] & /@ order2Eq

Then, the condition for the elimination of the terms that produce secular terms, solvability condition, demands that STerms1 be orthogonal to the adjoint, left eigenvector. The result is

SCond1=Solve[Conjugate[leftVec] . STerms1 == 0,Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]

whose complex conjugate is

ccSCond1 = SCond1 /. conjugateRule

With this solvability condition, order2Eq becomes

\(order2Eqm = \(order2Eq /. SCond1\) /. ccSCond1; 

    order2Eqm /. displayRule\)

whose particular solution can be expressed as

DSolve[order2Eqm, {(Subscript[u, 2])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], (Subscript[v, 2])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]}, 
      timeScales]

Here, directly applying DSolve to solve for order2Eqm does not work. Since the differential operators on the left-hand sides of order2Eqm are essentially ordinary-differential operators in terms of (Subscript[T, 0]), we can first transform order2Eqm into an ordinary-differential-equation form and then use DSolve to obtain the particular solution as

sol2pForm = 
      DSolve[
                order2Eqm /. {(Subscript[u, 2]) -> ((Subscript[u, 2])[#1] &), 
                    (Subscript[v, 2]) -> ((Subscript[v, 2])[#1] &)}, {(Subscript[u, 2])[(Subscript[T, 0])], (Subscript[v, 2])[(Subscript[T, 0])]}, 
                (Subscript[T, 0])][[1]] /. C[_] -> 0 // TrigToExp // ExpandAll

We can then express the solution of the second-order equations in a pure function form as

order2Sol = {
        (Subscript[u, 2]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol2pForm[[1, 2]] // Evaluate], 
        (Subscript[v, 2]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol2pForm[[2, 2]] // Evaluate]}

Third-Order Equations

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

\(order3Eq = \(eqOrder[3] /. order1Sol\) /. order2Sol // ExpandAll; 

    order3Eq /. displayRule\)

To eliminate the terms that lead to secular terms from order3Eq, we calculate the vector proportional to (E^(I (Subscript[T, 0]))) in their right-hand sides and obtain

STerms2 = 
      Coefficient[#[[2]], Exp[I (Subscript[T, 0])]] & /@ order3Eq

Hence, the solvability condition at this order is

SCond2 = Conjugate[leftVec] . STerms2 == 0 // ExpandAll

This solvability condition is the same as that obtained by attacking the second-order form of the governing equation.

3.4.3 First-Order Complex-Valued System

In this section, we first transform eq31g into a single first-order complex-valued equation using the transformation

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

whose inverse is given by

\(zetaRule = 
      \(Solve[transfRule /. Rule -> Equal, {ζ[t], ζ\&_[t]}]\)[
        \([1]\)]\)

where (Overscript[ζ,  _]) is the complex conjugate of ζ as defined by the complex conjugate rule

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

It follows from eq31g that the acceleration is given by

acceleration = Solve[eq31g, (u^″)[t]][[1]]

Differentiating zetaRule[[1]] once with respect to t and using acceleration and transfRule, we obtain the following first-order complex-valued equation:

eq343a = 
      D[zetaRule[[1]], t] /. acceleration /. transfRule /. 
          Rule -> Equal // ExpandAll

To determine a second-order uniform expansion of the solution of eq343a using the method of multiple scales, we first introduce

multiScales={ζ[t] -> ζ[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])],\(ζ\&_[t] -> ζ\&_[T\_0, T\_1, T\_2]\),ζ'[t]->dt[1][ζ[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]]};

and then expand ζ and \(TraditionalForm\`ζ\&_\) in the form

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

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

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

Equating coefficients of like powers of ϵ yields

eqEps = 
      CoefficientList[Subtract @@ eq343b, ϵ] == 0 // Thread // 
        Rest;

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

eqOrder[i_] := 
      (eqEps[[1, 1]] /. (Subscript[ζ, 1]) -> (Subscript[ζ, i])) - # & /@  
          eqEps[[i]] // Reverse

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

Array[eqOrder, 3] /. displayRule // TableForm

The general solution of the first-order equation, eqOrder[1], can be expressed in a pure function form as

sol1 = 
      DSolve[eqOrder[1], (Subscript[ζ, 1]), timeScales][[1]] /. 
        C[1] -> A

whose complex conjugate is

ccsol1 = sol1 /. conjugateRule

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

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

    order2Eq /. displayRule\)

Eliminating the terms that lead to seculr terms in (Subscript[ζ, 2]) from order2Eq yields

SCond1 = Coefficient[order2Eq[[2]], (E^(I (Subscript[T, 0])))] == 0

or

SCond1Rule=Solve[SCond1,Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]

With this solvability condition, the particular solution of order2Eq can be expressed in a pure function form as

sol2 = 
      DSolve[order2Eq /. SCond1Rule, (Subscript[ζ, 2]), timeScales][[1]] /. 
        C[_] -> (0 &)

whose complex conjugate is

ccsol2 = sol2 /. conjugateRule;

Substituting the first- and second-order solutions in the third-order equation, eqOrder[3], we obtain

\(order3Eq = 
      \(\(\(eqOrder[3] /. sol1\) /. ccsol1\) /. sol2\) /. ccsol2 // 
        ExpandAll; 
order3Eq /. displayRule\)

Eliminating the terms that produce secular terms in (Subscript[ζ, 3]) from order3Eq demands the following condition

\(SCond2 = Coefficient[order3Eq[\([2]\)], E\^\(I\ T\_0\)] == 0; 

    SCond2 /. displayRule\)

which is in full agreement with that obtained by attacking directly eq31g.

3.5 The Method of Averaging

As before, we use the method of variation of parameters to transform the dependent variable u into a and β according to

eq35a = u[t] == ϵ a Cos[t + β];

subject to the constraint

eq35b = D[eq35a, t]

where a and β are functions of t

tdependentRule = {a -> a[t], β -> β[t]};

and ϵ is a small dimensionless parameter that is a measure of the amplitude of oscillation. Differentiating eq35a with respect to t yields

eq35c = D[eq35a /. tdependentRule, t]

Comparing eq35b and eq35c, we conclude that

cond35a = 
      Expand[eq35c[[2]] - (eq35b[[2]] /. tdependentRule)] == 0

Differentiating eq35b once with respect to t yields

eq35d = D[eq35b /. tdependentRule, t]

Substituting eq35a and eq35d into eq31g, we have

cond35b = 
      eq31g /. ({eq35a /. tdependentRule, eq35d} /. Equal -> Rule) // 
        ExpandAll

Solving cond35a and cond35b for (a^′)[t] and (β^′)[t], we have

transformedEq = 
      Solve[{cond35a, cond35b}, {(a^′)[t], 
                (β^′)[t]}][[1]] /. Rule -> Equal // 
        Simplify

where a is assumed to be different from zero in arriving at transformedEq[[2]]. Using trigonometric identities, we expand the right-hand sides of transformedEq in Fourier series and obtain

ampEq1 = TrigReduce /@ transformedEq[[1]]

phaseEq1 = TrigReduce /@ transformedEq[[2]]

Since ϵ is small, (a^′)[t] and (β^′)[t] are slowly varying functions of t. Then, according to the method of averaging, we can neglect the fast-varying terms in ampEq1 and phaseEq1 and obtain the following averaged equations:

ampAvgEq = ampEq1 /. Sin[_.  t + _] -> 0

phaseAvgEq = phaseEq1 /. Cos[_.  t + _] -> 0

Whereas ampAvgEq agrees with ampEq, phaseAvgEq does not agree with phaseEq obtained in the preceding section by using the method of multiple scales. There is a term (5/12) (ϵ^2) (a^2) (Subscript[α, 2]^2) missing from phaseAvgEq. Following the details of the solution in the preceding section, one finds that this term is the result of the interaction of the first- and second-order approximations. This interaction was not taken into account in arriving at ampAvgEq and phaseAvgEq. To include the effect of this interaction, we need to carry out the solutions of ampEq1 and phaseEq1 to higher order. This is accomplished by using the generalized method of averaging, which is discussed next, or its variant the Krylov-Bogoliubov-Mitropolsky technique, which is discussed in Section 3.7.

3.6 The Generalized Method of Averaging

To apply this method, we introduce the variable ϕ = t + β and rewrite ampEq1 and phaseEq1 as

ampEq2 = 
      ampEq1 /. 
        Sin[n_.  t + a_] :> 
          Sin[Expand[n (ϕ[t] - β[t]) + a]]

phaseEq2 = 
      (ϕ^′)[t] == 1 + phaseEq1[[2]] /. 
        Cos[n_.  t + a_] :> 
          Cos[Expand[n (ϕ[t] - β[t]) + a]]

We seek approximate solutions to ampEq2 and phaseEq2 in the form

asol = 
      a[t] == (Subscript[a, 0])[t] + 
          Sum[(ϵ^i) (Subscript[a, i])[(Subscript[a, 0])[t], (Subscript[ϕ, 0])[t]], {i, 2}]

phisol = 
      ϕ[t] == 
        (Subscript[ϕ, 0])[t] + 
          Sum[(ϵ^i) (Subscript[ϕ, i])[(Subscript[a, 0])[t], (Subscript[ϕ, 0])[t]], {i, 2}]

where (Subscript[a, 0])[t] and (Subscript[ϕ, 0])[t] are expanded in power series in ϵ as

a0Eq = 
      ((Subscript[a, 0])^′)[t] == Sum[(ϵ^i) (Subscript[A, i])[(Subscript[a, 0])[t]], {i, 2}]

phi0Eq = 
      ((Subscript[ϕ, 0])^′)[t] == 
        1 + Sum[(ϵ^i) (Subscript[Φ, i])[(Subscript[a, 0])[t]], {i, 2}]

The functions (Subscript[a, 1]), (Subscript[a, 2]), …, (Subscript[a, n]) and (Subscript[ϕ, 1]), (Subscript[ϕ, 2]), …, (Subscript[ϕ, n]) are fast varying functions of (Subscript[ϕ, 0]), while it follows from a0Eq and phi0Eq that (Subscript[a, 0]), and hence the (Subscript[A, n]) and (Subscript[Φ, n]) are slowly varying functions of t.

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

aEq1 = 
      (D[asol, t] /. ({a0Eq, phi0Eq} /. Equal -> Rule) // ExpandAll)
           /. (ϵ^(n_ /; n > 2)) -> 0

phiEq1 = 
      (D[phisol, t] /. ({a0Eq, phi0Eq} /. Equal -> Rule) // ExpandAll)
           /. (ϵ^(n_ /; n > 2)) -> 0

Next, we substitute asol and phisol into ampEq2 and phaseEq2, expand the right-hand sides for small ϵ keeping terms up to O((ϵ^2)), and rewrite their right-hand sides as

ampEqrhs = 
      Series[ampEq2[[2]] /. ({asol, phisol} /. Equal -> Rule), {
            ϵ, 0, 2}] // Normal

phaseEqrhs = 
      Series[phaseEq2[[2]] /. ({asol, phisol} /. Equal -> Rule), {
            ϵ, 0, 2}] // Normal

Equating coefficients of like powers of ϵ from aEq1, ampEqrhs, phiEq1, and phaseEqrhs, we obtain

ampEq3 = 
      CoefficientList[#, ϵ] & /@ 
            (aEq1[[2]] == ampEqrhs) // Thread // Rest

phaseEq3 = 
      CoefficientList[#, ϵ]& /@ 
            (phiEq1[[2]] == phaseEqrhs) // Thread // Rest

Next, we use the method of separation of variables to separate fast and slowly varying terms in the first-order equations, ampEq3[[1]] and phaseEq3[[1]]. The slowly varying parts yield

cond36a = {(Subscript[A, 1])[(Subscript[a, 0])[t]] -> 0, (Subscript[Φ, 1])[(Subscript[a, 0])[t]] -> 0};

whereas the fast varying parts yield

order1Eq = {ampEq3[[1]], phaseEq3[[1]]} /. cond36a

To determine the solutions of order1Eq using DSolve, we transform them from partial-differential to ordinary-differential equations, replace (Subscript[ϕ, 0])[t] with s, and obtain

order1EqTransf = 
      order1Eq /. {(Subscript[a, 1]) -> ((Subscript[a, 1])[#2] &), 
              (Subscript[ϕ, 1]) -> ((Subscript[ϕ, 1])[#2] &)} /. f_[t] -> f /. 
        (Subscript[ϕ, 0]) -> s

Instead of replacing (Subscript[ϕ, 0])[t], we could have symbolized it. Using DSolve, we find that (Subscript[a, 1]) and (Subscript[ϕ, 1]) are given by

rule36[1] = 
      DSolve[order1EqTransf[[1]], (Subscript[a, 1])[s], s][[1]] /. 
            C[_] -> 0 /. {s -> (Subscript[ϕ, 0])[t], (Subscript[a, 0]) -> (Subscript[a, 0])[t]} /. 
        (Subscript[a, 1])[arg_] -> (Subscript[a, 1])[(Subscript[a, 0])[t], arg]

rule36[2] = 
      DSolve[order1EqTransf[[2]], (Subscript[ϕ, 1])[s], s][[1]] /. 
            C[_] -> 0 /. {s -> (Subscript[ϕ, 0])[t], (Subscript[a, 0]) -> (Subscript[a, 0])[t]} /. 
        (Subscript[ϕ, 1])[arg_] -> (Subscript[ϕ, 1])[(Subscript[a, 0])[t], arg]

which can be combined into

cond36b = Join[rule36[1], rule36[2]];

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

order2Eq = 
      {ampEq3[[2]], phaseEq3[[2]]} /. cond36a /. cond36b

Next, we use trigonometric identities to expand the right-hand sides of the second-order equations in Fourier series and obtain

order2Eqrhs = 
      (#[[2]] // TrigReduce // Expand) & /@ order2Eq

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

cond36c = {(Subscript[A, 2])[(Subscript[a, 0])[t]], (Subscript[Φ, 2])[(Subscript[a, 0])[t]]} -> 
          (order2Eqrhs /. {Sin[_] -> 0, Cos[_] -> 0}) // Thread

Substituting the fast-varying components of the first-order solution, cond36b, into the expansions for a and ϕ and discarding terms of order higher than ϵ, we obtain, to the second approximation, that

asolF = asol /. cond36b /. (ϵ^2) -> 0

phisolF = phisol /. cond36b /. (ϵ^2) -> 0

Substituting the solw-varying components of the first- and second-order solutions, cond36a and cond36c, into a0Eq and phi0Eq yields

a0Eq1 = a0Eq /. cond36a /. cond36c

phi0Eq1 = phi0Eq /. cond36a /. cond36c

It follows from the above equations that (Subscript[a, 0]) =  constant and

phi0Rule = 
      DSolve[phi0Eq1 /. (Subscript[a, 0])[t] -> (Subscript[a, 0]), (Subscript[ϕ, 0])[t], t][[1]] /. 
        C[1] -> (Subscript[β, 0])

where (Subscript[β, 0]) is a constant. Substituting these expansions into eq35a, we obtain the following second-order expansion:

eq36a = 
      u[t] == ϵ a[t] Cos[ϕ[t]] /. 
        ({asolF, phisolF} /. Equal -> Rule)

To compare the present solution with those obtained by using the method of multiple scales and the Lindstedt-Poincaré technique, we expand the circular functions in eq36a for small ϵ about (Subscript[ϕ, 0]) and obtain

u[t] == 
      (Series[eq36a[[2]], {ϵ, 0, 2}] // Normal // 
            TrigReduce // Expand)

This expansion and phi0Rule are in full agreement with those obtained by using the method of multiple scales.

3.7 The Krylov-Bogoliubov-Mitropolsky Technique

In this section, we describe a variant of the generalized method of averaging, namely the Krylov-Bogoliubov-Mitropolsky technique. It is often referred to as the asymptotic method. When the nonlinear terms are neglected, the solution of eq31g is

linearSol = u == ϵ a Cos[t + β];

where a and β are constants and ϵ is a small dimensionless parameter that is a measure of the amplitude. When the nonlinear terms are included, we consider linearSol to be the first term in an approximate solution of eq31g but with slowly varying rather than constant a and β. Moreover, we introduce the fast scale ϕ = t + β and use a to represent the slow variations. Thus, we seek a second-order uniform expansion of the solution of eq31g in the form

uSol = 
      u[a[t], ϕ[t]] == 
        ϵ a[t] Cos[ϕ[t]] + 
          Sum[(ϵ^i) (Subscript[u, i])[a[t], ϕ[t]], {i, 2, 3}]

Again, terms of order (ϵ^3) have been included in the expansion, but these terms will not be included in the final result. Since a and β are slowly varying functions of t, we express them in power series of ϵ in terms of a as

apEq = a'[t] == Sum[(ϵ^i) (Subscript[A, i])[a[t]], {i, 2}]

phipEq = 
      ϕ'[t] == 
        1 + Sum[(ϵ^i) (Subscript[Φ, i])[a[t]], {i, 2}]

In what follows, we need (a^″)[t] and (ϕ^″)[t]. To this end, we differentiate apEq and phipEq with respect to t, use these equations to express (a^′)[t] and (ϕ^′)[t] in terms of the (Subscript[A, i])[a[t]] and (Subscript[Φ, i])[a[t]], discard terms of order higher than (ϵ^2), and obtain

appEq = 
      (D[apEq, t] /. (apEq /. Equal -> Rule) // ExpandAll) /. 
        (ϵ^(n_ /; n > 2)) -> 0

phippEq = 
      (D[phipEq, t] /. (apEq /. Equal -> Rule) // ExpandAll) /. 
        (ϵ^(n_ /; n > 2)) -> 0

Thus, this method can be viewed as a multiple scales procedure with a and ϕ being the scales.

Substituting the assumed expansion for u[t], expressed in function form, and the expansions for the derivatives of a[t] and ϕ[t] into eq31g, expanding the result for small ϵ, and keeping terms up to order (ϵ^3), we have

eq37a = 
      (eq31g /. u -> Function[t, uSol[[2]] // Evaluate] /. 
              ({apEq, phipEq, appEq, phippEq} /. Equal -> Rule) // 
            ExpandAll) /. (ϵ^(n_ /; n > 3)) -> 0

Equating coefficients of like powers of ϵ in eq37a yields

eqEps = 
      Table[Coefficient[eq37a[[1]], (ϵ^i)] == 0, {i, 2, 3}] /. 
        f_[t] -> f

Next, we use trigonometric identities to expand the nonhomogeneous terms in eqEps[[1]] in a Fourier series and obtain

order2Eq = TrigReduce /@ eqEps[[1]]

Eliminating the terms that produce secular terms in (Subscript[u, 2]) from order2Eq demands that

rule37a = 
      Solve[Coefficient[order2Eq[[1]], {Cos[ϕ], Sin[ϕ]}] == 
              0 // Thread, {(Subscript[Φ, 1])[a], (Subscript[A, 1])[a]}][[1]]

Then, the solution of the resulting order2Eq can be written as

rule37b = 
      DSolve[order2Eq /. rule37a, (Subscript[u, 2])[a, ϕ], {a, ϕ}][
            [1]] /. C[_] -> (0 &) // Simplify

which can be expressed in function form as

rule37c = {
        (Subscript[u, 2]) -> Function[{a, ϕ}, rule37b[[1, 2]] // Evaluate]}

Substituting rule37a and rule37c into eqEps[[2]], we obtain

eq37b = eqEps[[2]] /. rule37a /. rule37c

Eliminating the terms that produce secular terms in (Subscript[u, 3]) from eq37b demands that

eq37c = 
      Coefficient[
            eq37b[[1]] // TrigReduce, {Cos[ϕ], Sin[ϕ]}] == 0 // 
        Thread

or

rule37d = Solve[eq37c, {(Subscript[A, 2])[a], (Subscript[Φ, 2])[a]}][[1]]

Substituting rule37a and rule37d into apEq and phipEq yields

apEq /. (rule37a~Join~rule37d /. a -> a[t])

phipEq /. (rule37a~Join~rule37d /. a -> a[t]) // ExpandAll

which are in full agreement with those obtained by using the method of multiple scales and the generalized method of averaging.

3.8 The Method of Normal Forms

In this section, we use the method of normal forms to determine a second-order uniform expansion of the solution of eq31g. To this end, we start with the corresponding first-order complex-valued equation

eq343a=ζ'[t]==\(I\ ζ[t] + 1\/2\ I\ α\_2\ ζ[t]\^2 + 
            1\/2\ I\ α\_3\ ζ[t]\^3 + 
            I\ α\_2\ ζ[t]\ ζ\&_[t] + 
            3\/2\ I\ α\_3\ ζ[t]\^2\ ζ\&_[t] + 
            1\/2\ I\ α\_2\ ζ\&_[t]\^2 + 
            3\/2\ I\ α\_3\ ζ[t]\ ζ\&_[t]\^2 + 
            1\/2\ I\ α\_3\ ζ\&_[t]\^3\);

According to the method of normal forms, we introduce the near-identity transformation

\(\(basicTerms = {η[t], η\&_[t]}; \)\)

\(zetaRule = {

        ζ -> 
          Function[t, 
            ϵ\ η[t] + 
                Sum[ϵ\^j\ h\_j\ @@\ basicTerms, {j, 2, 3}] // 
              Evaluate], 

        ζ\&_ -> 
          Function[t, 
            ϵ\ η\&_[t] + 
                Sum[ϵ\^j\ h\&_\_j\ @@\ basicTerms, {j, 2, 3}] // 
              Evaluate]}\)

that results in the simplest possible equation

etaRule = {
        (η^′)[t] -> 
          I η[t] + Sum[(ϵ^j) (Subscript[g, 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 zetaRule, into eq343a, using etaRule, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3),  we have

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

Equating coefficients of like powers of ϵ in eq38a yields

eqEps = 
      Table[Coefficient[Subtract @@ eq38a, (ϵ^i)] == 0, {i, 2, 3}]

Second-Order Solution

We choose the (Subscript[h, i]) to eliminate as many terms from eqEps, thereby reducing them into their simplest possible form. It turns out that we can eliminate all nonresonance terms. To determine the resonance terms in eqEps[[1]], 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

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

Hence, the possible resonance terms are given by

ResonantQTerm = 
      ((E^(-I t)) possibleQTerms /. form /. 
                (E^(_ t)) -> 0) possibleQTerms // Union // Rest

Consequently, there are no resonance terms in the second-order problem and the nonresonance terms can be defined as

NRQT = Complement[possibleQTerms, ResonantQTerm]

We associate with each of them a coefficient according to

coeffsQ = Table[(Subscript[Γ, j]), {j, Length[NRQT]}]

Therefore, (Subscript[h, 2]) and its complex conjugate have the form

\(hFormQ = {
        h\_2 -> \((
            Evaluate[coeffsQ . NRQT /. Thread[basicTerms -> {#1, #2}]]\ &)\), 
        h\&_\_2 -> 
          \((Evaluate[
                \(coeffsQ . NRQT /. conjugateRule\) /. 
                  Thread[basicTerms -> {#1, #2}]]\ &)\)}\)

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

coeffsQRule = 
      Solve[Coefficient[eqEps[[1, 1]] /. hFormQ, NRQT] == 0 // Thread, 
          coeffsQ][[1]]

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

gRuleQ = 
      (Subscript[g, 1])[t] -> 
        - Coefficient[eqEps[[1, 1]], ResonantQTerm] . ResonantQTerm

In this case, there are on resonance terms and hence (Subscript[g, 1]) = 0.

Third-Order Equations

Substituting the second-order results into the third-order equation yields

order3expr = 
      eqEps[[2, 1]] /. hFormQ /. coeffsQRule /. 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 = 
      ((E^(-I t)) possibleCTerms /. form /. 
                (E^(_ t)) -> 0) possibleCTerms // Union // Rest

Then, the nonresonance terms are the complement of these resonance terms; that is,

NRCT = Complement[possibleCTerms, ResonantCTerm]

Again, we associate coefficients with these terms as follows:

coeffsC = Table[(Subscript[Λ, j]), {j, Length[NRCT]}]

Consequently, (Subscript[h, 3]) and its complex conjugate have the forms

\(hFormC = {
        h\_3 -> \((
            Evaluate[coeffsC . NRCT /. Thread[basicTerms -> {#1, #2}]]\ &)\), 
        h\&_\_3 -> 
          \((Evaluate[
                \(coeffsC . NRCT /. conjugateRule\) /. 
                  Thread[basicTerms -> {#1, #2}]]\ &)\)}\)

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

coeffsCRule = 
      Solve[Coefficient[order3expr /. hFormC, NRCT] == 0 // Thread, coeffsC]
          [[1]]

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

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

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

\(moduEq = 
      2\ I\ \(E\^\(\(-I\)\ t\)\) 
              \((\(η\^′\)[t] - 
                  \((\(\(\(η\^′\)[t] /. etaRule\) /. gRuleQ\) /. 
                      gRuleC)\))\) == 0 /. {
            η -> \((A[#]\ E\^\(I\ #\)&)\), 
            η\&_ -> \((A\&_[#]\ E\^\(\(-I\)\ #\)&)\)} // ExpandAll\)

which is in full agreement with that obtained by using the other techniques.