Chapter 4
Forced Oscillations of the Duffing Equation

Off[General::spell1]

In contrast with the preceding chapters, which deal with free oscillations, this chapter and the next two deal with forced oscillations. We consider

\(\(FDuffingEq = 
      \(u\^″\)[t] + ω\^2\ u[t] + 
          2  μ\&^\ \(u\^′\)[t] + δ\ u[t]\^2 + 
          α\ u[t]\^3 == F\ Cos[Ω\ t]; \)\)

where \(TraditionalForm\`μ\&^\) is a positive constant. In this chapter, we determine second-order solutions to FDuffingEq beginning with the straightforward expansion in the next section. We investigate this straightforward expansion and determine under what conditions it breaks down. This leads to the so-called resonance values of Ω. In Section 4.2, we use the method of multiple scales to determine second-order uniform expansions of the solutions of FDuffingEq for all resonance cases, including the effect of light viscous damping. In Sections 4.3 and 4.4, we use the generalized method of averaging and the method of normal forms, respectively, to determine second-order uniform expansions for the case of subharmonic resonance of order one-half; that is, Ω ≈ 2  ω.

4.1 Straightforward Expansion

We seek a straightforward expansion for the solutions of FDuffingEq in the form

uSol = u -> (Sum[(ϵ^i) (Subscript[u, i])[#], {i, 3}] &);

Substituting uSol into FDuffingEq , letting \(TraditionalForm
      \`F = \(ϵ\ f\ \ and\ \ \(\ μ\)\&^ = 
          ϵ\ μ\)\), expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we have

\(eq41a = 
      \((\(FDuffingEq /. uSol\) /. {F -> ϵ\ f, 
                μ\&^ -> ϵ\ μ} // ExpandAll)\) /. 
        ϵ\^\(n_ /; n > 3\) -> 0\)

Equating coefficients of like powers of ϵ on both sides yields

eqEpsa = 
      CoefficientList[Subtract @@ eq41a, ϵ] == 0 // Thread // 
        Rest

The particular solution of the first-order equation, linear and undamped problem, can be expressed as

sol1p = 
      DSolve[eqEpsa[[1]], (Subscript[u, 1])[t], t][[1]] /. C[_] -> 0 // 
          TrigToExp // Simplify

We note that sol1p consists of a term whose denominator is very small when Ω ≈ ω. Such a term is called a small-divisor term. Moreover, u[t] tends to infinity as Ω → ω, and the excitation is referred to as a resonance excitation. Because the small-divisor term appears in the first-order problem, we speak of a primary or main resonance. When the small-divisor terms appear in the higher-order problems, we speak of secondary resonances. In the case of primary resonance, the scaling F → ϵ f is not valid. Physically, as soon as the motion becomes large, the damping and nonlinearity are activated to counter the effect of the resonance. Consequently, to obtain a uniform expansion in this case, we rescale F and \(TraditionalForm\`μ\&^\), as discussed in Section 4.2, so that the influence of the damping and nonlinearity balances the influence of the primary resonance.

When Ω is away from ω, we add the homogeneous solution to sol1p and obtain the general solution of the first-order problem as

fRule = {
        f -> 2  Λ ((Ω^2) - (ω^2))};

\(sol1 = 
      u\_1 -> Function[t, 
          \(A\ E\^\(I\ ω\ t\) + A\&_\ E\^\(\(-I\)\ ω\ t\) + 
                  sol1p[\([1, 2]\)] /. fRule // Expand\) // Evaluate]\)

where A is a complex-valued constant. Substituting sol1 into the second-order equation, eqEpsa[[2]], yields

order2Eqa = 
      ((Subscript[u, 2])^″)[t] + (ω^2) (Subscript[u, 2])[t] - # & /@ 
              eqEpsa[[2]] /. sol1 // ExpandAll // Reverse

Solving for the particular solution of order2Eqa, we obtain

sol2p = 
      ((Subscript[u, 2])[t] /. 
                    DSolve[order2Eqa[[1]] == #, (Subscript[u, 2])[t], t][[1]]
                       /. C[_] -> 0 // TrigToExp // Simplify) & /@ 
        order2Eqa[[2]]

Clearly, sol2p breaks down because it contains secular terms proportional to the damping coefficient μ. Moreover, sol2p breaks down when Ω ≈ 2  ω, Ω ≈ ω/2, or Ω ≈ 0  because they produce small-divisor terms and hence secondary resonances. As discussed in Section 4.2, the case Ω ≈ 2  ω is called subharmonic resonance of order one-half and the case Ω ≈ ω/2 is called superharmonic resonance of order two.

When Ω is away from ω, 2 ω, ω/2, or 0, sol2p is valid only if we rescale the damping term so that it first appears in the third-order equation; that is, \(TraditionalForm\`μ\&^ -> \(ϵ\^2\) μ\). Using this scaling, we rewrite eq41a as

eq41b = 
      eq41a /. μ -> ϵ μ /. (ϵ^4) -> 0

eqEpsb = 
      CoefficientList[Subtract @@ eq41b, ϵ] == 0 // Thread // 
        Rest

Comparing eqEpsb to eqEpsa, we immediately obtain its solution from sol2p as

sol2pNew = sol2p /. μ -> 0

Hence, the general solution of eqEpsb[[2]] can be written as

\(sol2New = 
      u\_2 -> Function[t, 
          B\ E\^\(I\ ω\ t\) + B\&_\ E\^\(\(-I\)\ ω\ t\) + 
              sol2pNew // Evaluate]\)

where B is a complex-valued constant.

Substituting sol1 and sol2New into eqEpsb[[3]] yields

order3Eqb = 
      ((Subscript[u, 3])^″)[t] + (ω^2) (Subscript[u, 3])[t] - # & /@ 
                eqEpsb[[3]] /. sol1 /. sol2New // ExpandAll // Reverse

The particular solution of order3Eqb is given by

(sol3p = 
          ((Subscript[u, 3])[t] /. 
                        DSolve[order3Eqb[[1]] == #, (Subscript[u, 3])[t], t][
                          [1]] /. C[_] -> 0 // TrigToExp // Simplify)
                 & /@ order3Eqb[[2]]) // Timing

Clearly, sol3p breaks down because it contains secular terms and small-divisor terms when Ω ≈ ω, Ω ≈ 0, Ω ≈ 2  ω, Ω ≈ ω/2, Ω ≈ 3  ω, and Ω ≈ ω/3. As shown in the next section, the small-divisor terms arising from Ω ≈ 3  ω produce a subharmonic reosonance of order one-third and those arising from Ω ≈ ω/3 produce a superharmonic resonance of order three. Carrying out the expansion to higher order, one finds that other resonances may occur.

We note that the resonances that occur depend on the order of the nonlinearity. Quadratic nonlinearities produce (a) subharmonic resonances of order one-half and superharmonic resonances of order two at second order and (b) subharmonic resonances of order one-third and superharmonic resonances of order three at third order. On the other hand, cubic nonlinearities produce subharmonic resonances of order one-third and superharmonic resonances of order three at third order. For a given system and order of approximation, the resonances produced can be easily identified by carrying out a straightforward expansion as done above.

In the next three sections, we use the method of multiple scales, generalized method of averaging, and method of normal forms to determine second-order uniform expansions for the solutions of FDuffingEq that do not contain secular or small-divisor terms.

4.2 The Method of Multiple Scales

4.2.1 Preliminaries

For a uniform second-order expansion, we need three time scales, which we symbolize and list as follows:

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

The maximum order of the expansion is related, in general, to the number of time scales by

maxOrder = Length[timeScales] - 1;

In what follows, we need the complex conjugates of A, ζ, and Γ, which we define by

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

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

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

\(\(displayRule = {
        \(\(Derivative[a__]\)[ζ_\_i_]\)[__] :> 
          SequenceForm[
            Times\ @@\ MapIndexed[D\_\(#2[\([1]\)] - 1\)\%#1\ &, {a}], 
            ζ\_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\&_], ζ_\_i_[__] -> ζ\_i, A[__] -> A, 
        A\&_[__] -> A\&_}; \)\)

Using the time scales (Subscript[T, 0]), (Subscript[T, 1]),  and (Subscript[T, 2]), we transform FDuffingEq from an ordinary-differential equation into a partial-differential equation according to the rule

multiScalesRule = {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])]], 
        t -> (Subscript[T, 0])};

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

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

For a uniform expansion, we need to eliminate the secular and small-divisor terms. To accomplish this, we need to distinguish between primary and secondary resonances. They are treated separately beginning with the case of primary resonance.

4.2.2 Primary Resonance

We scale the damping \(TraditionalForm\`μ\&^\) and forcing F so that the damping term and resonance terms appear at the same order at which the secular and / or small-divisor terms generated by the nonlinearities appear. In this case, they appear at order (ϵ^3),  and hence we scale the forcing and damping as

\(\(scaling = {F -> ϵ\^3\ f, μ\&^ -> ϵ\^2\ μ}; 
    \)\)

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

\(eq422a = 
      \((\(\(FDuffingEq /. scaling\) /. multiScalesRule\) /. solRule // 
            ExpandAll)\) /. ϵ\^\(n_ /; n > maxOrder + 1\) -> 0; 

    eq422a /. displayRule\)

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

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

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

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

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

Array[eqOrder, 3] /. displayRule

The general solution of eqOrder[1] can be expressed as

\(\(sol1 = 
      u\_1 -> Function[{T\_0, T\_1, T\_2}, 
          A[T\_1, T\_2]\ Exp[I\ ω\ T\_0] + 
            A\&_[T\_1, T\_2]\ Exp[\(-I\)\ ω\ T\_0]]; \)\)

Using sol1, we rewrite eqOrder[2] as

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

    order2Eq /. displayRule\)

Eliminating the terms that produce secular terms in order2Eq demands that

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

Solving SCond1 for Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])] yields

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

Hence, A = A[(Subscript[T, 2])].

Substituting the solvability condition into order2Eq, we have

\(order2Eqm = 
      \(order2Eq /. SCond1Rule\) /. \((SCond1Rule /. conjugateRule)\); 

    order2Eqm /. displayRule\)

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

\(sol2p = 
      \(\(DSolve[order2Eqm, u\_2\ @@\ timeScales, timeScales]\)[\([1]\)] /. 
            C[_] -> \((0\ &)\) // TrigToExp\) // ExpandAll; 

    sol2p /. displayRule\)

Next, we use sol2p[[1,2]] 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])}, sol2p[[1, 2]] // Evaluate];

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

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

    order3Eq /. displayRule\)

To express the nearness of Ω to ω, we introduce a detuning parameter σ = O(1) defined by

OmgRule = {Ω -> ω + (ϵ^2) σ};

Using the OmgRule, we convert the small-divisor term arising from the excitation into a secular term. Then, eliminating the secular terms from (Subscript[u, 3]) demands that

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

\(SCond2 = 
      Coefficient[order3Eq[\([2]\)] /. expRule1, Exp[I\ ω\ T\_0]] == 
        0; 
SCond2 /. displayRule\)

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

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

4.2.3 Secondary Resonances Due to Cubic Nonlinearities

In this case Ω is away from ω, 2  ω, and ω/2 and small-divisor terms first appear at O((ϵ^3)). We scale the damping \(TraditionalForm\`μ\&^\) and forcing F so that the damping term and resonance terms appear at the same order as the cubic nonlinearity according to

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

\(eq423a = 
      \((\(\(FDuffingEq /. scaling\) /. multiScalesRule\) /. solRule // 
            ExpandAll)\) /. ϵ\^\(n_ /; n > maxOrder + 1\) -> 0; 

    eq423a /. displayRule\)

Equating coefficients of like powers of ϵ, we obtain

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

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

eqOrder[i_] := 
      (eqEps[[1, 1]] /. f -> 0 /. (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

The particular solution of eqOrder[1] can be expressed as

sol1p = 
      DSolve[eqOrder[1], (Subscript[u, 1]) @@ timeScales, timeScales][[1]] /. 
          C[_] -> (0&) // Simplify

The general solution of eqOrder[1] can be written as

fRule = {
        f -> 
          -2 Λ 
            (-(ω^2) + (Ω^2))};

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

Substituting sol1 into eqOrder[2] yields

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

    order2Eq /. displayRule\)

The Case Ω ≈ 3 ω

To express the nearness of Ω to 3  ω, we introduce a detuning parameter σ = O(1) defined by

OmgRule = {Ω -> 3  ω + (ϵ^2) σ};

expRule1 = {
        Exp[arg_] :> 
          Exp[Expand[arg /. OmgRule] /. (ϵ^2) (Subscript[T, 0]) -> (Subscript[T, 2])]};

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

SCond1 = 
      Coefficient[order2Eq[[2]] /. expRule1, Exp[I ω (Subscript[T, 0])]] == 0

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

whose complex conjugate is

ccSCond1Rule = SCond1Rule /. conjugateRule

Substituting the solvability condition, SCond1Rule, and its complex conjugate, into order2Eq, we have

(order2Eqm = order2Eq /. SCond1Rule /. ccSCond1Rule) /. 
      displayRule

In order to efficiently use DSolve to determine the particular solution of order2Eqm, we first transform the partial-differential equation to an ordinary-differential equation. Then we solve for the particular solution of the resulting equation and obtain

\(sol2p = 
      \(\((\(\(u\_2[T\_0] /. 
                    \(DSolve[
                        \((order2Eqm[\([1]\)] /. u\_2 -> \((u\_2[#1]\ &)\))
                            \) == #, u\_2[T\_0], T\_0]\)[\([1]\)]\) /. 
                  C[_] -> 0 // TrigToExp\) // Expand)\)\ &\)\ /@\ 
        order2Eqm[\([2]\)]; 
sol2p /. displayRule\)

which 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])}, sol2p // Evaluate];

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

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

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

\(SCond2 = 
      Coefficient[
            \(order3Eq[\([2]\)] /. expRule1\) /. 
              Ω -> 3  ω, Exp[I\ ω\ T\_0]] == 0 /. 
        D[SCond1Rule, T\_1]; 
SCond2 /. displayRule\)

The Case Ω ≈ (1/3)ω

To express the nearness of Ω to (1/3) ω, we introduce a detuning parameter σ = O(1) defined by

OmgRule = {
        Ω -> 
          (1/3) (ω + (ϵ^2) σ)};

expRule1 = {
        Exp[arg_] :> 
          Exp[Expand[arg /. OmgRule] /. (ϵ^2) (Subscript[T, 0]) -> (Subscript[T, 2])]};

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

SCond1 = 
      Coefficient[order2Eq[[2]] /. expRule1, Exp[I ω (Subscript[T, 0])]] == 0

or

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

whose complex conjugate is

ccSCond1Rule = SCond1Rule /. conjugateRule;

Substituting the solvability condition, SCond1Rule, and its complex conjugate, into order2Eq, we have

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

    order2Eqm /. displayRule\)

Again, we first transform the partial-differential equation into an ordinary-differential equation and then use DSolve to obtain the particular solution of order2Eqm as

(sol2p = 
          ((Subscript[u, 2])[(Subscript[T, 0])] /. 
                      DSolve[
                          (order2Eqm[[1]] /. (Subscript[u, 2]) -> ((Subscript[u, 2])[#1] &))
                               == #, (Subscript[u, 2])[(Subscript[T, 0])], (Subscript[T, 0])][[1]] /. 
                    C[_] -> 0 // TrigToExp) & /@ order2Eqm[[2]])
         /. displayRule

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

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

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

Converting the terms that produce small-divisor terms into terms that produce secular terms and then eliminating the secular terms from (Subscript[u, 3]) demands that

\(SCond2 = 
      Coefficient[
            \(order3Eq[\([2]\)] /. expRule1\) /. 
              Ω -> \(1\/3\) ω, 
            Exp[I\ ω\ T\_0]] == 0 /. D[SCond1Rule, T\_1]; 

    SCond2 /. displayRule\)

4.2.4 Secondary Resonances Due to Quadratic Nonlinearities

In this case Ω is away from ω and small-divisor terms first appear at O((ϵ^2)). We scale the damping \(TraditionalForm\`μ\&^\) and forcing F so that the damping term and resonance terms appear at the same order as the quadratic nonlinearity; that is,

\(\(scaling = {F -> ϵ\ f, μ\&^ -> ϵ\ μ}; \)\)

\(eq424a = 
      \((\(\(FDuffingEq /. scaling\) /. multiScalesRule\) /. solRule // 
            ExpandAll)\) /. ϵ\^\(n_ /; n > maxOrder + 1\) -> 0; 

    eq424a /. displayRule\)

Equating coefficients of like powers of ϵ, we obtain

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

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

eqOrder[i_] := 
      (eqEps[[1, 1]] /. f -> 0 /. (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

The particular solution of eqOrder[1] can be expressed as

sol1p = 
      DSolve[eqOrder[1], (Subscript[u, 1]) @@ timeScales, timeScales][[1]] /. 
          C[_] -> (0&) // Simplify

fRule = 
      f -> -2 (-(ω^2) + (Ω^2)) 
          Λ;

The general solution of eqOrder[1] can be written in pur function form as

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

Substituting sol1 into eqOrder[2] yields

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

    order2Eq /. displayRule\)

The Case Ω ≈ 2 ω

To express the nearness of Ω to 2  ω, we introduce a detuning parameter σ = O(1) defined by

OmgRule = {Ω -> 2 ω + ϵ σ};

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

expRule2 = {
        Exp[arg_] :> 
          Exp[Expand[arg /. OmgRule] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]};

Converting the terms that produce small-divisor terms into terms that produce secular terms and then eliminating the secular terms from (Subscript[u, 2]) demands that

\(SCond1 = 
      Coefficient[order2Eq[\([2]\)] /. expRule2, Exp[I\ ω\ T\_0]] == 
        0; 
SCond1 /. displayRule\)

or

SCond1Rule=(Solve[SCond1,Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll);
SCond1Rule /. displayRule

Substituting sigRule into SCond1Rule yields

\(SCond1Rulem = 
      SCond1Rule /. 
        Exp[a_] :> 
          Exp[\(a /. sigRule\) /. T\_1 -> ϵ\ T\_0 // Expand]; 

    SCond1Rulem /. displayRule\)

whose complex conjugate is

ccSCond1Rulem = SCond1Rulem /. conjugateRule;

Substituting these conditions into order2Eq yields

\(order2Eqm = \(order2Eq /. SCond1Rulem\) /. ccSCond1Rulem // ExpandAll; 
    
order2Eqm /. displayRule\)

Transforming order2Eqm into an ordinary-differential equation, we can obtain the particular solution as

\(sol2p = 
      \(\((\(u\_2[T\_0] /. 
                  \(DSolve[
                      \((order2Eqm[\([1]\)] /. u\_2 -> \((u\_2[#1]\ &)\))
                          \) == #, u\_2[T\_0], T\_0]\)[\([1]\)]\) /. 
                C[_] -> 0 // TrigToExp)\)\ &\)\ /@\ order2Eqm[\([2]\)]; 

    sol2p /. displayRule\)

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

Then, using sol1 and sol2 in eqOrder[3], we have

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

Converting the terms that produce small-divisor terms into terms that produce secular terms and then eliminating the secular terms from (Subscript[u, 3]) demands that

\(SCond2 = 
      \(\(Coefficient[
                  \(order3Eq[\([2]\)] /. expRule2\) /. 
                    Ω -> 2  ω, 
                  Exp[I\ ω\ T\_0]] == 0 /. D[SCond1Rule, T\_1]\) /. 
            SCond1Rule\) /. \((SCond1Rule /. conjugateRule)\) // ExpandAll; 

    SCond2 /. displayRule\)

or

SCond2Rule=Solve[SCond2,Derivative[0, 1][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll;
SCond2Rule /. displayRule

The two partial-differential equations, SCond1Rule and SCond2Rule, can be reconstituted to obtain an ordinary-differential equation governing A. The result is

(2 I ω (A^′) == 
          (2 I ω dt[1][A[(Subscript[T, 1]), (Subscript[T, 2])]] /. SCond1Rule /. 
                  SCond2Rule // Expand // Collect[#, ϵ] &))
         /. displayRule

4.2.5 First-Order Real-Valued System

As in Section 4.2.4 , we determine a second-order uniform expansion of the solution of the FDuffingEq for the case of subharmonic resonance of order one-half. However, instead of treating the second-order form of this equation, we treat its corresponding first-order form; that is, we first transform it into a system of two real-valued first-order equations. To this end, we introduce the transformation

vRule = (u^′)[t] -> v[t];

In order that the influence of the nonlinearity and damping balance the subharmonic resonance, we scale the damping coefficient \(TraditionalForm\`μ\&^\) and forcing amplitude F as

\(\(scaling = {F -> ϵ\ f, μ\&^ -> ϵ\ μ}; \)\)

Using the scaling and vRule, we transform the FDuffingEq into the following system of two real-valued first-order equations:

eq425a = {vRule /. Rule -> Equal, 
        FDuffingEq /. scaling /. vRule /. D[vRule, t]}

To determine a second-order uniform expansion of the solution of eq425a using the method of multiple scales, we first transform these equations from ordinary-differential equations into partial-differential equations in terms of the three time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]) according to

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])]],t -> (Subscript[T, 0])};

Next, we 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 eq425a, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

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

Equating coefficients of like powers of ϵ in eq425b yields

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

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

\(linearSys = \(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f -> 0; 

    linearSys /. displayRule\)

\(eqOrder[i_] := 
      \((linearSys /. u_\_1 -> u\_i)\) == 
          \((linearSys /. u_\_1 -> u\_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: Linear System

We seek a solution of the 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}]

whose 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 = {c1, 1} /. 
        Solve[(hermitian[coefMat] . {c1, 1})[[1]] == 0, c1][
          [1]]

The complex conjugate of leftVec is

ccleftVec = leftVec /. conjugateRule

Thefore, the homogeneous solution of eqOrder[1] is given by

sol1h = 
      # + (# /. conjugateRule) & /@ 
        (rightVec A[(Subscript[T, 1]), (Subscript[T, 2])] (E^(I ω (Subscript[T, 0]))))

Transforming eqOrder[1] into ordinary-differential equations and using DSolve, we obtain the particular solution of eqOrder[1] as

sol1p = 
      DSolve[
                eqOrder[1] /. {(Subscript[u, 1]) -> ((Subscript[u, 1])[#1] &), 
                    (Subscript[v, 1]) -> ((Subscript[v, 1])[#1] &)}, {(Subscript[u, 1])[(Subscript[T, 0])], (Subscript[v, 1])[(Subscript[T, 0])]}, 
                (Subscript[T, 0])][[1]] /. C[_] -> 0 // TrigToExp // Simplify

Then, the general solution of the first-order equations can be expressed in pure function form as

fRule = 
      f -> -2 (-(ω^2) + (Ω^2)) 
          Λ;

sol1 = {
        (Subscript[u, 1]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol1h[[1]] + 
                ((Subscript[u, 1])[(Subscript[T, 0])] /. sol1p /. fRule // Expand) // 
              Evaluate], 
        (Subscript[v, 1]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol1h[[2]] + 
                ((Subscript[v, 1])[(Subscript[T, 0])] /. sol1p /. fRule // Expand) // Evaluate]}

Second-Order Equations

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

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

    order2Eq /. displayRule\)

To express quantitatively the nearness of the subharmonic resonance of order one-half, we introduce the detuning parameter σ defined by

OmgRule = {Ω -> 2  ω + ϵ σ};

To convert the terms that produce small-divisor terms into terms that produce secular terms, we use the rule

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

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

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

Then, the solvability condition of order2Eq demands that STerms1 be orthogonal to the solution of the adjoint; that is, ccleftVec. The result is

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

whose complex conjugate is

ccSCond1 = SCond1 /. conjugateRule

Next, we use the solvability condition to eliminate (Subscript[D, 1]) A and its complex conjugate from order2Eq and then find their particular solution. To simplify the resulting equations using Mathematica, we need to replace (Subscript[T, 1]) with ϵ (Subscript[T, 0]). To this end, we express the detuning parameter σ in terms of ω and Ω as

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

and define the rule

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

This rule enables us to rewrite SCond1 and ccSCond1 as

\(SCond1f = Join[SCond1, ccSCond1] /. expRule2; 

    SCond1f /. displayRule\)

With this form of the solvability condition, order2Eq becomes

\(order2Eqm = order2Eq /. SCond1f // ExpandAll; 

    order2Eqm /. displayRule\)

Next, we use the method of undetermined coefficients to determine the particular solution of order2Eqm. To accomplish this, we first identify the form of their nonhomogeneous terms, which we will refer to as possible terms, and then seek the solution as a linear combination of them. Substituting this solution into order2Eqm and equating the coefficient of each possible term to zero, we obtain a set of pairs of algebraic equations that need to be solved for the unknown coefficients. The solutions of these pairs are unique except for the coefficients of (E^(I ω (Subscript[T, 0]))) and (E^(-I ω (Subscript[T, 0]))). To determine a unique solution corresponding to these terms, we use two different approaches.

Approach 1:

In the first approach, we replace ω with another symbol, say (Subscript[ω, 0]), find the solution of the system of algebraic equations, and then take the limit as (Subscript[ω, 0]) tends to ω. To accomplish this, we define

\(\(collectForm = {E\^\(I\ T\_0\ ω\)\ A[T\_1, T\_2], 
        E\^\(\(-I\)\ T\_0\ ω\)\ A\&_[T\_1, T\_2], 
        E\^\(I\ T\_0\ Ω\)\ Λ, 
        E\^\(\(-I\)\ T\_0\ Ω\)\ Λ}; \)\)

Then, the form of the possible terms on the right-hand sides of order2Eqm can be determined as

possibleTerms = 
      (collectForm /. ω -> (Subscript[ω, 0]))~Join~
        (Outer[Times, collectForm, collectForm] // Flatten // Union)

Next, we associate with each possible term an undetermined coefficient in (Subscript[u, 2]) by

symbolList1 = Table[(Subscript[Γ, i]), {i, Length[possibleTerms]}]

and an undetermined coefficient in (Subscript[v, 2]) by

symbolList2 = symbolList1 /. Γ -> Λ

Then, we seek a particular solution of order2Eqm in the form

sol2Rule = {
        (Subscript[u, 2]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            symbolList1 . possibleTerms // Evaluate], 
        (Subscript[v, 2]) -> Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            symbolList2 . possibleTerms // Evaluate]};

Substituting sol2Rule into order2Eqm and replacing (E^(±I ω (Subscript[T, 0]))) on the right-hand sides of order2Eqm by (E^(±I (Subscript[ω, 0]) (Subscript[T, 0]))), we have

order2Eqf = 
      Subtract @@ # & /@ order2Eqm /. sol2Rule /. {
            Exp[I ω (Subscript[T, 0])] -> Exp[I (Subscript[ω, 0]) (Subscript[T, 0])], 
            Exp[-I ω (Subscript[T, 0])] -> Exp[-I (Subscript[ω, 0]) (Subscript[T, 0])]} // 
        Expand;

Applying directly the function Coefficient to collect the coefficients of all possibleTerms could result in extra terms that are functions of (Subscript[T, 0]) , which we eliminate. The result is

algEqs = 
      Flatten[Coefficient[order2Eqf, #]&/@ possibleTerms] == 0 /. 
          Exp[_ (Subscript[T, 0])] -> 0 // Thread

Using the function Solve, we directly solve the system of algebraic equations algEqs for these undetermined coefficients and obtain

var = Join[symbolList1, symbolList2];

coefs = Solve[algEqs, var][[1]];

Substituting coefs into sol2Rule and changing (Subscript[ω, 0]) back to ω, we obtain the solution of the second-order equations in pure function form as

\(sol2a = 
      sol2Rule /. 
        Function[{T\_0, T\_1, T\_2}, b_] :> 
          Function[{T\_0, T\_1, T\_2} // Evaluate, 
            \(b /. coefs\) /. ω\_0 -> ω // Evaluate]; 

    sol2a /. displayRule\)

Approach 2:

As discussed above, solutions of the algebraic equations resulting from applying the method of undetermined coefficients are unique except those corresponding to (E^(I ω (Subscript[T, 0]))) and (E^(-I ω (Subscript[T, 0]))). To render them unique, we demand that they be orthogonal to solutions of the corresponding adjoint problems; that is, we damand that the coefficient vector of (E^(I ω (Subscript[T, 0]))) be {(Subscript[u, 2 ]), (Subscript[v, 2])} be orthogonal to ccleftVec. To accomplish this, we first determine the form of all possible terms on the right-hand sides of order2Eqm and then identify those among them corresponding to the resonance terms and those corresponding to nonresonance terms.

possibleTerms = 
      collectForm~Join~
        (Outer[Times, collectForm, collectForm] // Flatten // Union)

Next, we identify the resonance terms among them by

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

Hence, the nonresonance part of possibleTerms is the complement of RT; that is,

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

Next, we associate with each possible resonance term an undetermined coefficient by using the rule

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

and we associate with each possible nonresonance term an undetermined coefficient by using the rule

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

Hence, the coefficients of RT in the nonhomogeneous terms in order2Eqm are

coefRT = 
      Coefficient[#[[2]], RT] & /@ order2Eqm /. 
        Exp[_ (Subscript[T, 0]) + _. ] -> 0

It follows from

ccleftVec . coefRT === {0, 0}

that coefRT is orthogonal to the solution of the adjoint, ccleftVec, as it should.

Imposing the condition that {(Subscript[u, 2]) , (Subscript[v, 2])} is orthogonal to ccleftVec, we seek (Subscript[u, 2]) and (Subscript[v, 2]) in pure function form as

\(\(sol2Form = {
		
        u\_2 -> Function[{T\_0, T\_1, T\_2}, 
            RTsymbolList . RT + \((RTsymbolList . RT /. conjugateRule)\) + 
                NRTsymbolList[1] . NRT // Evaluate], 
        v\_2 -> Function[{T\_0, T\_1, T\_2}, 
            \(-I\)\ ω\ RTsymbolList . RT + 
                \((\(-I\)\ ω\ RTsymbolList . RT /. conjugateRule)\) + 
                NRTsymbolList[2] . NRT // Evaluate]}; \)\)

Substituting sol2Form into order2Eqm and collecting coefficients of NRT, we obtain the following set of algebraic equations:

algEqs1 = 
      (Coefficient[Subtract @@ #, NRT] & /@ 
                  (order2Eqm /. sol2Form) /. Exp[_ (Subscript[T, 0])] -> 0 // Flatten)
             == 0 // Thread

Solving algEqs1 for the undetermined coefficients NRTsymbolList yields

coef1 = 
      Solve[algEqs1, Array[NRTsymbolList, 2] // Flatten][[1]]

Substituting sol2Form into either of the two order2Eqm and collecting the coefficients of RT, we obtain the following set of algebraic equations:

algEqs2 = 
      Coefficient[Subtract @@ order2Eqm[[1]] /. sol2Form, RT] == 0 /. 
          Exp[_ (Subscript[T, 0])] -> 0 // Thread

Solving algEqs2 for the undetermined coefficients RTsymbolList yields

coef2a = Solve[algEqs2, RTsymbolList][[1]]

coef2 = Join[coef2a, coef2a /. conjugateRule]

Substituting coef1 and coef2 into sol2Form, we obtain the solution of the second-order equations in pure function form as

\(sol2b = 
      sol2Form /. 
        Function[{T\_0, T\_1, T\_2}, b_] :> 
          Function[{T\_0, T\_1, T\_2} // Evaluate, 
            \(\(b /. coef1\) /. coef2 // Expand\) // Evaluate]; 

    sol2b /. displayRule\)

These two approaches generate the same results because

Subtract @@ ((Subscript[u, 2])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])] /. {sol2a, sol2b}) // Simplify

Subtract @@ ((Subscript[v, 2])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])] /. {sol2a, sol2b}) // Simplify

The second approach is more convenient to use for a general system of equations. Hence, we will adopt this approach in the remainder of this book.

Third-Order Equations

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

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

Using the expRule1 to convert the terms that produce small-divisor terms into terms that produce secular terms in the right-hand sides of order3Eq and collecting the terms that could produce secular terms, we have

\(STerms2 = 
      \(Coefficient[
            \(#[\([2]\)] /. expRule1\) /. Ω -> 2\ ω, 
            Exp[I\ ω\ T\_0]]\ &\)\ /@\ order3Eq; 

    STerms2 /. displayRule\)

Then, the solvability condition demands that STerms2 be orthogonal to the solution of the adjoint ccleftVec. The result is

SCond2=Solve[ccleftVec . STerms2 == 0,Derivative[0, 1][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll;
SCond2 /. displayRule

The two partial-differential equations, SCond1 and SCond2, can be reconstituted to obtain an ordinary-differential equation governing A. The result is

(2 I ω (A^′) == 
          (2 I ω dt[1][A[(Subscript[T, 1]), (Subscript[T, 2])]] /. SCond1 /. 
                  SCond2 // Expand // Collect[#, ϵ] &)) /. 
      displayRule

4.2.6 First-Order Complex-Valued System

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

\(\(transformRule = {u -> \((ζ + ζ\&_)\), 
        u\^′ -> I\ ω\ \((ζ - ζ\&_)\)}; \)\)

gRule = 
      g -> - 2 ϵ μ (u^′) - δ (u^2) - 
          α (u^3) + ϵ f Cos[Ω t];

then FDuffingEq becomes

eq426a = 
      (ζ^′) == I ω ζ - (I/(2 ω)) g /. 
            gRule /. transformRule // ExpandAll

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

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

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

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

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

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

Equating coefficients of like powers of ϵ yields

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

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

eqOrder[i_] := 
      (eqEps[[1, 1]] /. f -> 0 /. (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 pure function form as

sol1Form = 
      (DSolve[eqOrder[1], (Subscript[ζ, 1])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], timeScales][
                  [1, 1, 2]] /. C[1] -> A // Expand) /. 
          Exp[a_] :> Exp[a // Expand] //. 
        c1_ Exp[x_] + c2_ Exp[x_] :> Factor[c1 + c2] Exp[x]

\(sol1 = {
        ζ\_1 -> Function[{T\_0, T\_1, T\_2}, sol1Form // Evaluate], 
        ζ\&_\_1 -> 
          Function[{T\_0, T\_1, T\_2}, sol1Form /. conjugateRule // Evaluate]}
        \)

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

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

    order2Eq /. displayRule\)

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

SCond1 = 
      Coefficient[order2Eq[[2]] /. expRule1, Exp[I ω (Subscript[T, 0])]] == 0

or

SCond1Rule=(Solve[SCond1,Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll)/.c1_ μ + c2_ μ :> Factor[c1 + c2] μ

With this solvability condition, we rewrite order2Eq as

order2Eqm = 
      Expand /@ 
        (order2Eq /. SCond1Rule /. 
            Exp[arg_] :> 
              Exp[arg /. (Subscript[T, 1]) -> ϵ (Subscript[T, 0]) /. sigRule // Expand])
          ;

There are at least two approaches we can use to solve for the particular solution of order2Eqm. In the first approach, we transform order2Eqm into an ordinary-differential equation and then use DSolve and Map to determine the following particular solution:

(sol2a = 
          ((Subscript[ζ, 2])[(Subscript[T, 0])] /. 
                        DSolve[
                            (order2Eqm[[1]] /. 
                                  (Subscript[ζ, 2]) -> ((Subscript[ζ, 2])[#1] &))
                                 == #, (Subscript[ζ, 2])[(Subscript[T, 0])], (Subscript[T, 0])][[1]]
                           /. C[_] -> 0) & /@ order2Eqm[[2]] /. 
              Exp[a_] :> Exp[Expand[a]] // Expand) // Timing

In the second approach, we multiply order2Eqm by an integrating factor and then use the function Integrate. To reduce the computation time, we define first the rule

intRule = int[a_ + b_, arg_] -> int[a, arg] + int[b, arg];

Then, using this rule and an integrating factor, we find a particular solution of order2Eqm as

(sol2b = 
          ((int[order2Eqm[[2]]*Exp[-I (Subscript[T, 0]) ω] // 
                            Expand, (Subscript[T, 0])] //. intRule /. int -> Integrate)
                    *Exp[I (Subscript[T, 0]) ω] // Expand) /. 
            Exp[a_] :> Exp[Expand[a]]) // Timing

The CPU time is decreased by a factor of

%%%[[1]]/%[[1]]

sol2a === sol2b

Then, the particular solution of order2Eqm can be expressed in pure function form as

\(\(sol2 = {
        ζ\_2 -> Function[{T\_0, T\_1, T\_2}, sol2b // Evaluate], 
        ζ\&_\_2 -> 
          Function[{T\_0, T\_1, T\_2}, sol2b /. conjugateRule // Evaluate]}; 
    \)\)

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

order3Eq = Expand /@ (eqOrder[3] /. sol1 /. sol2);

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

SCond2Rule=(Solve[Coefficient[order3Eq[[2]] /. expRule1, 
                            Exp[I ω (Subscript[T, 0])]] == 0,Derivative[0, 1][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]/.SCond1Rule/.(SCond1Rule /. conjugateRule)//ExpandAll)//.{c1_ μ + c2_ μ :> Factor[c1 + c2] μ, 
            c1_ (f^2) (δ^2) + c2_ (f^2) (δ^2) :> 
              Factor[c1 + c2] (f^2) (δ^2), 
            c1_ (f^2) α + c2_ (f^2) α :> 
              Factor[c1 + c2] (f^2) α};
SCond2Rule /. displayRule

or

fRule = {
        f -> 2 ((ω^2) - (Ω^2)) Λ};

SCond2Rule /. fRule /. Ω -> 2  ω /. 
      displayRule

which is in agreement with that obtained in the preceding section.

4.2.7 The Function MMS1

Collecting the steps described in Section 4.2.5 for a system of two real-valued first-order equations, we can build a function named MMS1 (Method of Multiple Scales for 1DOF system) specifically for FDuffingEq. A more general function (a Package) can be similarly created by considering as arguments the governing equation, symbols for the dependent variable, independent variable, excitation amplitudes and frequencies, and all other related quantities which allow the program to identify their respective meanings. We then use MMS1 to solve for different resonance cases.

MMS1[scaling_List, ResonanceCond : {__Equal}]:=
   Module[{},
		vRule = (u^′)[t] -> v[t];
     OmgRule = Solve[ResonanceCond, Ω][[1]];
		eqa = {vRule /. Rule -> Equal, 
                FDuffingEq /. scaling /. vRule /. D[vRule, t]};
		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])]],t -> (Subscript[T, 0])};
		solRule = {
                u -> (Sum[(ϵ^j) (Subscript[u, j])[#1, #2, #3], {j, 3}] &), 
                v -> (Sum[(ϵ^j) (Subscript[v, j])[#1, #2, #3], {j, 3}] &)}
                ;
		eqb = (
                  eqa /. multiScales /. solRule // TrigToExp // 
                    ExpandAll) /. (ϵ^(n_ /; n > 3)) -> 0;
		eqEps = 
              Thread[CoefficientList[Subtract @@ #, ϵ] == 0] &
                       /@ eqb // Transpose // Rest;
		
	    		
		linearSys = #[[1]] & /@ eqEps[[1]] /. f -> 0;
		\(eqOrder[i_] := 
              \((linearSys /. u_\_1 -> u\_i)\) == 
                  \((linearSys /. u_\_1 -> u\_i)\) - 
                    \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([i]\)])\) // Thread\);
		coefList = 
              (E^(-I ω (Subscript[T, 0]))) linearSys /. {
                    (Subscript[u, 1]) -> (P (E^(I ω #))&), 
                    (Subscript[v, 1]) -> (Q (E^(I ω #))&)} // Expand;
		coefMat = Outer[Coefficient, coefList, {P, Q}];
		hermitian[mat_?MatrixQ] := 
              mat /. conjugateRule // Transpose;
		rightVec = {1, c1} /. 
                Solve[(coefMat . {1, c1})[[1]] == 0, c1][[1]]
                  ;
		leftVec = {c1, 1} /. 
                Solve[(hermitian[coefMat] . {c1, 1})[[1]] == 0, c1]
                    [[1]];
		ccleftVec = leftVec /. conjugateRule;
		sol1h = 
              # + (# /. conjugateRule) & /@ 
                (rightVec A[(Subscript[T, 1]), (Subscript[T, 2])] (E^(I ω (Subscript[T, 0]))));
		sol1p = 
              DSolve[
                        eqOrder[1] /. {(Subscript[u, 1]) -> ((Subscript[u, 1])[#1] &), 
                            (Subscript[v, 1]) -> ((Subscript[v, 1])[#1] &)}, {(Subscript[u, 1])[(Subscript[T, 0])], 
                          (Subscript[v, 1])[(Subscript[T, 0])]}, (Subscript[T, 0])][[1]] /. C[_] -> 0 // 
                  TrigToExp // Simplify;
		fRule = 
              f -> -2 (-(ω^2) + (Ω^2)) 
                  Λ;
		sol1 = {
                (Subscript[u, 1]) -> 
                  Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
                    sol1h[[1]] + 
                        ((Subscript[u, 1])[(Subscript[T, 0])] /. sol1p /. fRule // Expand) // 
                      Evaluate], 
                (Subscript[v, 1]) -> 
                  Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
                    sol1h[[2]] + 
                        ((Subscript[v, 1])[(Subscript[T, 0])] /. sol1p /. fRule // Expand) // 
                      Evaluate]};
		
	    		
		order2Eq = eqOrder[2] /. sol1 // ExpandAll;
		expRule1 = 
              Exp[arg_] :> 
                Exp[Expand[arg /. OmgRule] /. 
                    (ϵ^n)_.  (Subscript[T, 0]) :> timeScales[[n + 1]]];
		STerms1 = 
              Coefficient[#[[2]] /. expRule1, Exp[I ω (Subscript[T, 0])]] &
                   /@ order2Eq;
		SCond1=Solve[ccleftVec . STerms1 == 0,Derivative[1, 0][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]];
		ccSCond1 = SCond1 /. conjugateRule;
		sigRule = 
              Solve[OmgRule /. Rule -> Equal, σ][[1]];
		expRule2 = 
              Exp[a_] :> 
                Exp[a /. {(Subscript[T, 1]) -> ϵ (Subscript[T, 0]), 
                          (Subscript[T, 2]) -> (ϵ^2) (Subscript[T, 0])} /. sigRule // Expand]
                  ;
		SCond1f = Join[SCond1, ccSCond1] /. expRule2;
		order2Eqm = order2Eq /. SCond1f // ExpandAll;
	  \(collectForm = 
              Join[{E\^\(I\ T\_0\ ω\)\ A[T\_1, T\_2], 
                  E\^\(\(-I\)\ T\_0\ ω\)\ A\&_[T\_1, T\_2]}, 
                If[\((F/f /. scaling)\) === ϵ, {
                    E\^\(I\ T\_0\ Ω\)\ Λ, 
                    E\^\(\(-I\)\ T\_0\ Ω\)\ 
                      Λ}, {}]]\);
		possibleTerms = 
              collectForm~Join~
                (Outer[Times, collectForm, collectForm] // Flatten // 
                    Union);
		RT = (
                      # /. {a_ /; a =!= 0 -> 1} & /@ 
                        ((E^(-I ω (Subscript[T, 0]))) possibleTerms /. 
                              expRule1 /. Exp[_ (Subscript[T, 0]) + _. ] -> 0)) 
                    possibleTerms // Union // Rest;
		NRT = Complement[possibleTerms, Join[RT, RT /. conjugateRule]]
                ;
		RTsymbolList = Table[(Subscript[Γ, j]), {j, Length[RT]}];
		NRTsymbolList[i_] = 
              Table[(Subscript[Λ, i, j]), {j, Length[NRT]}];
		\(sol2Form = {
\ \ \ \ \ 		
                u\_2 -> 
                  Function[{T\_0, T\_1, T\_2}, 
                    RTsymbolList . RT + 
                        \((RTsymbolList . RT /. conjugateRule)\) + 
                        NRTsymbolList[1] . NRT // Evaluate], 
                v\_2 -> 
                  Function[{T\_0, T\_1, T\_2}, 
                    \(-I\)\ ω\ RTsymbolList . RT + 
                        \((\(-I\)\ ω\ RTsymbolList . RT /. 
                            conjugateRule)\) + NRTsymbolList[2] . NRT // 
                      Evaluate]}\);
		algEqs1 = 
              (Coefficient[Subtract @@ #, NRT] & /@ 
                          (order2Eqm /. sol2Form) /. Exp[_ (Subscript[T, 0])] -> 0 // 
                      Flatten) == 0 // Thread;
		coef1 = 
              Solve[algEqs1, Array[NRTsymbolList, 2] // Flatten][[1]]
                ;
		algEqs2 = 
              Coefficient[Subtract @@ order2Eqm[[1]] /. sol2Form, 
                      RT] == 0 /. Exp[_ (Subscript[T, 0])] -> 0 // Thread;
		coef2a = Solve[algEqs2, RTsymbolList][[1]];
		coef2 = Join[coef2a, coef2a /. conjugateRule];
		sol2 = 
              sol2Form /. 
                Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, b_] :> 
                  Function[{(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])} // Evaluate, 
                    b /. coef1 /. coef2 // Expand // Evaluate];
		
	    
		order3Eq = eqOrder[3] /. sol1 /. sol2 // ExpandAll;
		STerms2 = 
              Coefficient[
                    #[[2]] /. expRule1 /. 
                      (OmgRule /. ϵ -> 0), 
                    Exp[I ω (Subscript[T, 0])]] & /@ order3Eq;
		SCond2=Solve[ccleftVec . STerms2 == 0,Derivative[0, 1][A][(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll;
		
	    
		moduEq = 
              2 I ω (A^′) == 
                (2 I ω dt[1][ A[(Subscript[T, 1]), (Subscript[T, 2])]] /. SCond1
                         /. SCond2 // Collect[#, ϵ] &);
		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[
              u[t] == (
                  u @@ timeScales /. solRule /. 
                          ϵ^3 -> 0 /. sol1 /. sol2 /. 
                    displayRule)];
		If[(F/f /. scaling) === ϵ,Print["\!\(\*
StyleBox[\(where\),
FontColor->RGBColor[0, 0, 1]]\)\n Λ==f/(2((ω^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]
  ]

Primary Resonance:  Ω ≈ ω

\(\(scaling1 = {F -> ϵ\^3\ f, μ\&^ -> ϵ\^2\ μ}; 
    \)\)

ResonanceCond1 = {
        Ω == ω + (ϵ^2) σ};

MMS1[scaling1, ResonanceCond1] // Timing

Subharmonic Resonance of Order 1/2:  Ω ≈ 2ω

\(\(scaling2 = {F -> ϵ\ f, μ\&^ -> ϵ\ μ}; 
    \)\)

ResonanceCond2 = {
        Ω == 2  ω + ϵ σ};

MMS1[scaling2, ResonanceCond2] // Timing

Subharmonic Resonance of Order 1/3:  Ω ≈ 3ω

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

ResonanceCond3 = {
        Ω == 3  ω + (ϵ^2) σ};

MMS1[scaling3, ResonanceCond3] // Timing

Superharmonic Resonance of Order 2:  2Ω ≈ ω

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

ResonanceCond4 = {
        2  Ω == ω + ϵ σ};

MMS1[scaling4, ResonanceCond4] // Timing

Superharmonic Resonance of Order 3:  3Ω ≈ ω

\(\(scaling5 = {F -> ϵ\ f, μ\&^ -> ϵ\^2\ μ}; 
    \)\)

ResonanceCond5 = {
        3  Ω == ω + (ϵ^2) σ};

MMS1[scaling5, ResonanceCond5] // Timing

Slowly-Modulated Load:  Ω ≈ 0

\(\(scaling6 = {F -> ϵ\ f, μ\&^ -> ϵ\ μ}; 
    \)\)

ResonanceCond6 = {Ω == ϵ σ};

MMS1[scaling6, ResonanceCond6] // Timing

4.3 The Generalized Method of Averaging

To apply either the method of averaging or the generalized method of averaging to FDuffingEq, we need first to use the method of variation of parameters to transform it into a system of two first-order equations governing the amplitude and phase. To this end, we introduce the following transformation:

states = {a[t], ϕ[t]};

transformRule = {
        u -> a[t] Cos[ϕ[t]] + Λ Cos[Ψ[t]], 
        (u^′) -> 
          - ω a[t] Sin[ϕ[t]] - 
            Ω Λ Sin[Ψ[t]]};

gRule=g->-2 ϵ μ u'-ϵ δ (u^2)-(ϵ^2) α (u^3);

After some algebraic manipulations (Nayfeh, 1973, 1981), we transform FDuffingEq into two first-order equations as

eq43a = 
      D[states, t] == 
          ({-(g/ω) Sin[ϕ[t]], 
                        ω - (g/(a[t] ω)) Cos[ϕ[t]]} /. 
                      gRule /. transformRule // Expand // TrigReduce // 
              Expand) // Thread;

We seek approximate solutions to eq43a in the form

basicTerms = {η[t], φ[t], Ψ[t]};

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

where η[t] and φ[t] are expanded in power series in ϵ as

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

The functions (Subscript[a, 1]), (Subscript[a, 2]),  (Subscript[ϕ, 1]),  and (Subscript[ϕ, 2]) are fast varying functions of φ, while it follows from basicDRule that η, 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 solRule and basicDRule, expand the result for small ϵ, discard terms of order higher than (ϵ^2), and obtain

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

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

eq43bRHS = 
      CoefficientList[
            Series[#[[2]] /. solRule, {ϵ, 0, 2}] // Normal
                 // Expand, ϵ] & /@ eq43a;

Equating coefficients of like powers of ϵ in eq43bLHS and eq43bRHS, we obtain

eqEps = 
      MapThread[Equal, {eq43bLHS, eq43bRHS}, 2] // Transpose // Rest;

Next, we use the method of separation of variables to separate fast and slowly varying terms in the first-order equations, eqEps[[1]]. We first introduce the rules

SVT[j_] := {(Subscript[A, j]) @@ basicTerms, (Subscript[Φ, j]) @@ basicTerms}solVar[j_] := {(Subscript[a, j]), (Subscript[ϕ, j])}

psiRule = {φ[t] -> ω t + β[t], 
        Ψ[t] -> Ω t};

betaRule = 
      Solve[psiRule /. Rule -> Equal, {β[t], Ω}][
        [1]]

For the case of subharmonic resonance of O((1/2)), we define

OmgRule = {Ω -> 2  ω + ϵ σ};

sigRule = Solve[OmgRule /. Rule -> Equal, σ][[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]];

Using these rules, we find that 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, 
              SVT[1][[i]]] [[1]], {i, 2}] // Flatten

whereas the fast varying parts are given by

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

To determine a particular solution corresponding to these fast varying terms, we use the method of undetermined coefficients. To accomplish this, we first determine the possible forms of the terms in FVT1 as follows:

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

Using the principle of superposition, we seek a particular solution corresponding to the fast-varying terms as a linear combination of these possible forms:

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 equations for the undetermined coefficients, and then substituting the result back into sol1Form, we obtain the solution

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

which can be expressed in pure function form as

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

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

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

Since we are seeking an expansion valid up 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

SVT2Rule = 
      Table[Solve[
              order2Expr[[i]] == 0 /. 
                        Thread[solVar[2] -> (0 &)] /. expRule1 /. 
                    f_[_ t + _. ] -> 0 /. expRule2 /. 
                (OmgRule /. ϵ -> 0), SVT[2][[i]]], {i, 
              2}] // Flatten // ExpandAll

Hence, to the second approximation, we find  that

u[t] == 
      (TrigReduce /@ 
            (Normal[
                    Series[u /. transformRule /. solRule, {ϵ, 0, 
                        1}]] /. sol1Rule // Expand) // 
          Collect[#, ϵ]&)

where the equations governing the amplitude and phase are

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

4.4 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 FDuffingEq. To this end, we start with the corresponding first-order complex-valued equation, which was derived in Section 4.2.6 and a summary of the derivation is given below:

\(\(transformRule = {u[t] -> \((ζ[t] + ζ\&_[t])\), 
        \(u\^′\)[t] -> I\ ω\ \((ζ[t] - ζ\&_[t])\), 
        f\ Cos[Ω\ t] -> \(1\/2\) \((z[t] + z\&_[t])\)}; \)\)

gRule = 
      g -> -2 ϵ μ (u^′)[t] - 
          ϵ δ u[t]^2 - 
          (ϵ^2) α u[t]^3 + f Cos[Ω t];

eq44a = {
        (ζ^′)[t] == 
                I ω ζ[t] - I g/(2 ω) /. gRule
               /. transformRule // ExpandAll, 
        (z^′)[t] == I Ω z[t]}

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

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

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

that results in the simplest possible equation

etaRule = {
        (η^′)[t] -> 
          I ω η[t] + Sum[(ϵ^j) (Subscript[g, j])[t], {j, 2}], 
        (z^′)[t] -> I Ω z[t]}

where the overbar denotes the complex conjugate

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

Substituting the expansion for ζ, the zetaRule, into eq44a, using the etaRule, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^2),  we have

eq44b = 
      (eq44a[[1]] /. zetaRule /. etaRule /. 
              (etaRule /. conjugateRule) // ExpandAll) /. 
        (ϵ^(n_ /; n > 2)) -> 0;

Equating coefficients of like powers of ϵ in eq44b yields

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

First-Order Solution

Using the method of undetermined coefficients, we assume

coeffsL = Table[(Subscript[Δ, i]), {i, Length[basicTerms]}]

\(\(hFormL = {
        h\_0 -> \((
            Evaluate[
                coeffsL . basicTerms /. 
                  Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\), 
        h\&_\_0 -> 
          \((Evaluate[
                \(coeffsL . basicTerms /. conjugateRule\) /. 
                  Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\)}; \)\)

Substituting the hFormL into eqEps[[1]], collecting the coefficients of the basicTerms, and solving for the coeffsL, we obtain

coeffsLRule = 
      Solve[Coefficient[eqEps[[1, 1]] /. hFormL, basicTerms] == 0 // 
            Thread, coeffsL][[1]]

Substituting the coeffsLRule into the hFormL yields

hSolL = 
      (Subscript[h, 0]) @@ basicTerms /. hFormL /. coeffsLRule /. 
        (Subscript[Δ, i])_ -> 0

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

\(hRuleL = {
        h\_0 -> \((
            Evaluate[hSolL /. Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\), 
        h\&_\_0 -> 
          \((Evaluate[
                \(hSolL /. conjugateRule\) /. 
                  Thread[basicTerms -> {#1, #2, #3, #4}]]\ &)\)}\)

Second-Order Solution

Substituting the hRuleL into the left-hand side of the second-order equation, eqEps[[2]], yields

order2expr = eqEps[[2, 1]] /. hRuleL // Expand

We choose (Subscript[h, 1]) to eliminate as many terms from order2expr, thereby reducing it to its simplest possible form. It turns out that we can eliminate all nonresonance terms. To determine the resonance terms in order2expr, we first determine all of the possible forms of its nonhomogeneous terms:

possibleQTerms = 
      basicTerms~Join~
        (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\), 
        z[t] -> E\^\(I\ Ω\ t\), 
        z\&_[t] -> E\^\(\(-I\)\ Ω\ t\)}; \)\)

Hence, the possible resonance terms are given by

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

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

gRuleQ = 
      (Subscript[g, 1])[t] -> 
        -(Coefficient[order2expr /. (Subscript[h, 1]) -> (0&), ResonantQTerm] /. 
                Thread[basicTerms -> 0]) . ResonantQTerm

The nonresonance terms are the complement of the resonance terms; that is,

NRQT = Complement[possibleQTerms, ResonantQTerm];

We associate with each of them an undetermined coefficient according to

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

Therefore, (Subscript[h, 1]) has the form

hFormQ = {
        (Subscript[h, 1]) -> (
            Evaluate[coeffsQ . NRQT /. Thread[basicTerms -> {#1, #2, #3, #4}]]
                 &)}

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

coeffsQRule = 
      Solve[Coefficient[order2expr /. gRuleQ /. hFormQ // Expand, 
                  NRQT] == 0 /. Thread[basicTerms -> 0] // Thread, coeffsQ][
        [1]]

Substituting the coeffsQRule into the hFormQ yields

hSolQ = (Subscript[h, 1]) @@ basicTerms /. hFormQ /. coeffsQRule

Hence, we can write the second-order solution in pure function form as

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

Third-Order Equations

Substituting the hRuleL, hRuleQ, and gRuleQ into the left-hand side of the third-order equation, eqEps[[3]], yields

order3expr = 
      eqEps[[3, 1]] /. hRuleL /. hRuleQ /. gRuleQ /. 
          (gRuleQ /. conjugateRule) // Expand;

The nonhomogeneous terms in order3expr are proportional to

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

possibleTerms = Join[possibleQTerms, possibleCTerms] // Union

Next, we determine the resonance terms according to

ResonantCTerm = 
      ((E^(-I ω t)) possibleTerms /. form /. 
                  Ω -> 2 ω /. (E^(_ t)) -> 0) 
            possibleTerms // Union // Rest

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

gRuleC = 
      (Subscript[g, 2])[t] -> 
        -Factor[
              Coefficient[order3expr /. (Subscript[h, 2]) -> (0&), ResonantCTerm] /. 
                Thread[basicTerms -> 0]] . ResonantCTerm

The nonresonance terms are the complement of the resonance terms; that is,

NRCT = Complement[possibleTerms, ResonantCTerm];

Again, we associate an undetermined coefficient with each of these terms as follows:

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

Consequently, (Subscript[h, 2]) has the form

hFormC = {
        (Subscript[h, 2]) -> (
            Evaluate[coeffsC . NRCT /. Thread[basicTerms -> {#1, #2, #3, #4}]]
                 &)};

Substituting for (Subscript[h, 2]) 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 /. gRuleC /. hFormC // Expand, 
                    NRCT] == 0 /. Thread[basicTerms -> 0] // Thread, coeffsC]
            [[1]] // ExpandAll;

Substituting the coeffsCRule into the hFormC yields

hSolC = (Subscript[h, 2]) @@ basicTerms /. hFormC /. coeffsCRule;

Combining the etaRule, gRuleQ, and gRuleC and letting

\(\(solForm = {η -> \((A[#]\ E\^\(I\ ω\ #\)&)\), 
        η\&_ -> \((A\&_[#]\ E\^\(\(-I\)\ ω\ #\)&)\), 
        z -> \((f\ E\^\(I\ Ω\ #\)&)\), 
        z\&_ -> \((f\ E\^\(\(-I\)\ Ω\ #\)&)\)}; \)\)

we obtain the modulation equation

moduEq = 
      (2 I ω (E^(-I ω t)) 
                  ((η^′)[t] - 
                      ((η^′)[t] /. etaRule /. gRuleQ
                             /. gRuleC)) /. solForm // Expand // 
            Collect[#, ϵ]&) == 0

or

moduEq /. 
        f -> 2  Λ ((Ω^2) - (ω^2)) /. 
      Ω -> 2  ω

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