Chapter 7
Continuous Systems with Cubic Nonlinearities

Elastic systems, such as beams, plates, and shells, are usually modeled by partial-differential equations with specified boundary conditions. For small oscillations, the responses of such deformable bodies, continuous or distributed-parameter systems, can be adequately described by linear equations and boundary conditions. However, as the amplitude of oscillations increases, nonlinear effects in the governing equations, boundary conditions, or both come into play. The sources of nonlinearities may be geometric, inertial, material, or damping in nature.

Geometric nonlinearities may be due to one or more of the following: nonlinear relationships among the strains and the displacements, large rotations, free surfaces in fluids, time-dependent constraints, mid-plane stretching, and large curvatures. Nonlinear stretching of the midplane of a deformable body accompanies its transverse vibrations if it is supported in such a way as to restrict the movement of its ends and / or edges. If large-amplitude vibrations are accompanied by large changes in the curvature, it is necessary to employ a nonlinear relationship between the curvature and the displacement. Inertial nonlinearities are caused by concentrated and / or distributed masses, convective accelerations, and Corriolis and centripetal accelerations. Material nonlinearities occur whenever the constitutive relations are nonlinear, such as the case when the stresses are nonlinear functions of the strains; the resistive, inductive, and capacitive circuit elements are nonlinear; and the feedback control forces and moments in servomechanisms are nonlinear. Damping nonlinearities occur due to form drag and hysteresis.

Since exact solutions are, in general, not available for determining the dynamic responses of nonlinear continuous systems to external or parametric excitations, recourse has been made to approximate analyses by using either purely numerical techniques, or purely analytical techniques, or a combination of numerical and analytical techniques. Application of purely numerical techniques to such problems may be costly in terms of computation time and may not reveal some of their intricate and complicated responses. With purely analytical methods, it may be difficult to treat systems with inhomogeneities or complicated geometries. With a combination of numerical and analytical techniques, one can determine some of the complicated responses of systems with inhomogeneities and complicated shapes.

The numerical-analytic approaches can be divided into two groups: discretization and direct methods.  In the discretization methods, one postulates the solution in the form

w[x,t]=(∑_(m = 1)^M)(Subscript[ϕ, m])(x) (Subscript[q, m])(t)

where M is a finite integer. Then, one assumes the spatial functions (Subscript[ϕ, m])(x), space discretization, or the temporal functions (Subscript[q, m])(t), time discretization. With this discretization, the (Subscript[q, m])(t) are usually taken to be harmonic and the method of harmonic balance is used to obtain a set of nonlinear boundary-value problems for the (Subscript[ϕ, m])(x).

With space discretization, the (Subscript[ϕ, m])(x) (and thus the spatial dependence) are assumed a priori. If the boundary conditions are homogeneous, the (Subscript[ϕ, m])(x) are usually taken to be the eigenfunctions of the linearized problem. The method of weighted residuals or variational principles can then be used to determine a set of ordinary-differential equations governing the (Subscript[q, m])(t). The obtained set of ordinary-differential equations can be studied using any of a great number of methods developed for discrete systems.

The most common implementation of weighted residuals is the Galerkin method. The great majority of recent studies of forced vibrations assume that the response can be expressed in terms of only the linear modes that are directly or indirectly excited (Nayfeh and Mook, 1979). For example, if a system is driven near the natural frequency of a linear mode and that mode is not involved in an internal resonance with other modes, the response is assumed to consist of only that mode. Such an approach is usually referred to as a single-mode approximation.

In the direct approach, a reduction method, such as the method of multiple scales, is applied directly to the governing partial-differential equations and associated boundary conditions and no assumptions are made a priori regarding the spatial or temporal dependence of the response. The principal advantage of this approach is in the treatment of the boundary conditions at higher orders. This approach has been used recently by Nayfeh and Nayfeh (1979), Nayfeh (1975, 1996), Nayfeh and Asfar (1986), Nayfeh and Bouguerra (1990), Pai and Nayfeh (1990), Raouf and Nayfeh (1990), Nayfeh, Nayfeh, and Mook (1992), Nayfeh and Nayfeh (1993, 1995), Nayfeh, Nayfeh, and Pakdemirli (1995), Pakdemirli, Nayfeh, and Nayfeh (1995), Chin and Nayfeh (1996), Nayfeh and Lacarbonara (1997, 1998), Lacarbonara, Nayfeh, and Kreider (1998), Nayfeh, Lacarbonara, and Chin (1999), and Rega, et al. (1999).

Some of the aforementioned studies show that the discretization and direct approaches yield the same results for systems with cubic nonlinearities provided that first-order approximations are sought. For systems with quadratic and cubic nonlinearities, the discretization approach might produce quantitative, and in some case qualitative, errors unless many modes are included in the discretized model. In other words, one has to include as many terms in the discretized model as needed for convergence. Since we are dealing with systems with cubic nonlinearities in this chapter, we will obtain first-order results using both the discretization and direct approaches.

Preliminaries

Off[General::spell1, Integrate::gener]

Needs["Utilities`Notation`"]

To use the method of multiple scales, we introduce different time scales (Subscript[T, 0]) =  t ,  (Subscript[T, 1])= ϵ t, and (Subscript[T, 2])= (ϵ^2)t, symbolize them as

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

and form a list of them as follows:

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

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

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

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

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

To manipulate some complicated integrals without Mathematica being choked, we define the following rules:

\(intRule1 = {int[fun_, arg2__] :> int[Expand[fun], arg2]}; 

    intRule2 = {int[a_ + b_, arg2__] :> int[a, arg2] + int[b, arg2], 
        int[ϵ\^n_. \ fun_, arg2__] :> ϵ\^n\ int[fun, arg2], 
        int[a_\ fun_, a1_, b1___] :> 
          a\ int[fun, a1, b1] /; FreeQ[a, First[a1]], 
        int[int[a1_, a2_]\ fun_. , a3__] :> 
          int[a1, a2]\ int[fun, a3] /; FreeQ[Rest[a2], First[a2]]}; \)

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

\(\(displayRule = {
        \(\(Derivative[a_, b__]\)[w_\_i_]\)[x, __] :> 
          SequenceForm[
            If[arg1 = 
                Times\ @@\ MapIndexed[D\_\(#2[\([1]\)] - 1\)\%#1\ &, {b}]; 
              arg1 =!= 1, arg1, ""], 
            w\_i\^SequenceForm\ @@\ Table["'", {a}]], 
        \(\(Derivative[a__]\)[A\_i_]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\_i], \(\(Derivative[a__]\)[A\&_\_i_]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\&_\_i], w_\_i_[x, __] -> w\_i, A\_i_[__] -> A\_i, 
        A\&_\_i_[__] -> A\&_\_i, 
        Exp[a_.  + b_. \ Complex[0, m_]\ T\_0 + 
              c_. \ Complex[0, n_]\ T\_0] :> 
          Exp[a + \((m*b + n*c)\) I\ T\_0], int -> Integrate}; \)\)

7.1 Solvability Conditions and the Concept of Adjoint

In directly attacking continuous systems, one often encounters nonhomogeneous boundary-value problems whose homogeneous parts have nontrivial solutions. Consequently, the nonhomogeneous problems have solutions only if solvability or consistency conditions are satisfied (Nayfeh, 1981). In this section, we describe how to determine such solvability conditions.

We start the discussion with the simple problem

\(eq1 = 
      \(y\^″\)[x] + π\^2\ y[x] == π\ Sin[π\ x]; 

    bc1 = {y[0] == β\_1, y[1] == β\_2}; \)

The corresponding homogeneous problem has the nontrivial solution

DSolve[{eq1[[1]] == 0}~Join~bc1 /. (Subscript[β, i])_ -> 0, y[x], x][
      [1]]

Hence, the nonhomogeneous problem will not have a solution unless a solvability condition is satisfied. To determine this solvability condition, we use two approaches.

In the first approach, we find the general solution of eq1 using the function DSolve as

yRule = DSolve[eq1, y[x], x][[1]] // Simplify

where C[1] and C[2] are arbitrary constants. Imposing the boundary conditions bc1, we have

bc2 = bc1 /. Flatten[yRule /. {{x -> 0}, {x -> 1}}]

These equations are inconsistent unless

Equal @@ Plus @@ (List @@ # & /@ bc2)

which is the desired solvability condition.

In the second approach, instead of determining the general solution of the homogeneous differential equations and then enforcing the boundary conditions to determine the solvability conditions, we use the concept of adjoint as described next. This approach is attractive if one is not interested in determining the solution of the nonhomogeneous problem but interested only in determining the solvability conditions, as is the case in many applications of perturbation methods.

To determine the adjoint and then the solvability conditions, we need to perform integration by parts of some products. To accomplish this with Mathematica, we define a function named intByParts according to

\(intByParts[expr_ /; Head[expr] =!= Equal, depVar_] := 
\ \ \ 
      Module[{intRule1, intRule2, intRule3}, 
		\ \ 
        intRule1 = f_[terms_, x_] :> int[Expand[terms], x]; 
\ \ \ \ \ \ 
        intRule2 = int[term1_ + term2_, x_] :> int[term1, x] + int[term2, x]; 
        
\ \ \ \ \ \ 
        intRule3 = 
          c1_. \ int[term1_, x_] + c2_. \ int[term2_, x_] :> 
            int[c1\ term1 + c2\ term2, x]; 
\ \ \ \ \ \ 
        IBP[u_, v_, n_ /; n >= 1, range : {x_, a_, b_}] := 
          \((u\ \(\(Derivative[n - 1]\)[v]\)[x] /. x -> b)\) - 
            \((u\ \(\(Derivative[n - 1]\)[v]\)[x] /. x -> a)\) - 
            IBP[D[u, x], v, n - 1, range]; 
\ \ \ \ \ \ 
        IBP[u_, v_, 0, range : {x_, a_, b_}] := int[u\ v[x], range]; 
		

        \ \ \ \ \ \ 
        \(\(\(\(expr /. intRule1\) //. intRule2\) /. 
              int[u_\ \(\(Derivative[n_]\)[depVar]\)[x_], arg2_] :> 
                IBP[u, depVar, n, arg2]\) //. intRule3\) /. int -> Integrate

        \ \ ]\)

To determine the solvability condition of eq1 and bc1 with this approach, we multiply eq1 by u[x], integrate the result by parts from x = 0 to x = 1, and obtain

eq1a = intByParts[int[u[x] #, {x, 0, 1}], y] & /@ eq1

To determine the adjoint, we set the coefficient of y[x] in the integrand on the left-hand side of eq1a equal to zero; that is,

adjointEq1 = 
      Cases[eq1a[[1]], Integrate[expr_, _] :> Coefficient[expr, y[x]]]
            [[1]] == 0

To determine the adjoint boundary conditions, we consider the homogeneous problem (i.e., f[x] = 0, (Subscript[β, 1]) = 0, and (Subscript[β, 2]) = 0), use the adjointEq1, and obtain from eq1a that

bc1a = 
      eq1a /. Integrate -> (0 &) /. 
        (bc1 /. Equal -> Rule /. (Subscript[β, i])_ -> 0)

To determine the adjoint boundary conditions, we set each of the coefficients of (y^′)[0] and (y^′)[1] in bc1a equal to zero and obtain

adjointBC1=Solve[Coefficient[bc1a[[1]],{y'[0],y'[1]}]==0,{u[0], u[1]}][[1]]/.Rule -> Equal

Therefore, the adjoint u is defined by the adjoint system consisting of adjointEq1 and adjointBC1. Since they are the same as the homogeneous parts of eq1 and bc1, the problem is said to be self-adjoint. Hence, u[x] = Sin[x].

adjoint = {u ->  Function[x,  Sin[π x]]}

Once the adjoint problem has been defined, we return to the nonhomogeneous problem to determine the solvability condition. Substituting for the adjoint and the boundary conditions bc1 into eq1a yields the solvability condition

SolvCond1 = 
      eq1a[[2]] == 
        (eq1a[[1]] /. Integrate -> (0 &) /. 
              (bc1 /. Equal -> Rule) /. 
            (adjointBC1 /. Equal -> Rule))

SolvCondf = SolvCond1 /. adjoint

which is the same as the solvability condition obtained with the other method.

7.1.1 Hinged-Clamped Beam

In treating a hinged-clamped beam, we need to determine the adjoint of the boundary-value problem

eq1a=-(ω^2) ϕ[x]-2 P ϕ''[x]+Derivative[4][ϕ][x]==f[x];

bc1a={ϕ[0] == 0,ϕ''[0]==0,ϕ[1] == 0,ϕ'[1]==0};

Multiplying eq1a by u[x] and integrating the result by parts from x = 0 to x = 1, we obtain

eq1b = intByParts[int[u[x] #, {x, 0, 1}], ϕ] & /@ eq1a

To determine the adjoint equation, we set the coefficient of ϕ[x] in the integrand on the left-hand side of eq1b equal to zero; that is,

adjointEq1 = 
      Cases[eq1b[[1]], 
            Integrate[expr_, _] :> Coefficient[expr, ϕ[x]]][
          [1]] == 0

To determine the adjoint boundary conditions, we consider the homogeneous problem, use the adjointEq1, and obtain from eq1b that

bc1b = 
      eq1b /. Integrate -> (0 &) /. (bc1a /. Equal -> Rule)

To determine the adjoint boundary conditions, we set each of the coefficients of (ϕ^′)[0] , (ϕ^((3)))[0] , 
      (ϕ^″)[1], and (ϕ^((3)))[1] in bc1b equal to zero and obtain

adjointBC1=Solve[Coefficient[bc1b[[1]],{ϕ'[0],Derivative[3][ϕ][0],ϕ''[1],Derivative[3][ϕ][1]}]==0,{u[0], (u^″)[0], u[1], (u^′)[1]}
                ][[1]]/.Rule -> Equal

Therefore, the adjoint u is defined by the adjoint system consisting of adjointEq1 and adjointBC1. Since they are the same as the homogeneous parts of eq1a and bc1a, the problem is said to be self-adjoint.

Once the adjoint problem has been defined, we return to the nonhomogeneous problem to determine the solvability condition. Substituting for the adjoint and the boundary conditions bc1a into eq1b yields the solvability condition

SolvCond1 = 
      eq1b[[2]] == 
        (eq1b[[1]] /. Integrate -> (0 &) /. 
              (bc1a /. Equal -> Rule) /. 
            (adjointBC1 /. Equal -> Rule))

7.1.2 Cantilever Beam

The boundary-value problem for a cantilever beam can be written as

eq2a=-(ω^2) ϕ[x]+Derivative[4][ϕ][x]==f[x];

bc2a={ϕ[0] == 0,(ϕ^′)[0] == 0,(ϕ^″)[1] == 0,Derivative[3][ϕ][1]==0};

To determine the adjoint of this problem, we multiply eq2a by u[x], integrate the result by parts from x = 0 to x = 1, and obtain

eq2b = intByParts[int[u[x] #, {x, 0, 1}], ϕ] & /@ eq2a

We set the coefficient of ϕ[x] in the integrand on the left-hand side of eq2b equal to zero and obtain the adjoint equation as

adjointEq2 = 
      Cases[eq2b[[1]], 
            Integrate[expr_, _] :> Coefficient[expr, ϕ[x]]][
          [1]] == 0

To determine the adjoint boundary conditions, we consider the homogeneous problem, use the adjointEq2, and obtain from eq2b that

bc2b = 
      eq2b /. Integrate -> (0 &) /. (bc2a /. Equal -> Rule)

To determine the adjoint boundary conditions, we set each of the coefficients of (ϕ^″)[0] , 
      (ϕ^((3)))[0] , ϕ[1], and (ϕ^′)[1] in bc2b equal to zero and obtain

adjointBC2=Solve[Coefficient[bc2b[[1]],{(ϕ^″)[0],Derivative[3][ϕ][0],ϕ[1],(ϕ^′)[1]}]==0,{u[0],(u^′)[0],(u^″)[1],Derivative[3][u][1]}][[1]]/.Rule -> Equal

Therefore, the adjoint u is defined by the adjoint system consisting of adjointEq2 and adjointBC2. Since they are the same as the homogeneous parts of eq2a and bc2a, the problem is said to be self-adjoint.

Once the adjoint problem has been defined, we return to the nonhomogeneous problem to determine the solvability condition. Substituting for the adjoint and the boundary conditions bc2a into eq2b yields the solvability condition

SolvCond2 = 
      eq2b[[2]] == 
        (eq2b[[1]] /. Integrate -> (0 &) /. 
              (bc2a /. Equal -> Rule) /. 
            (adjointBC2 /. Equal -> Rule))

7.2 Hinged-Clamped Beam

7.2.1 EOM and BC's

We consider the nonlinear planar response of a hinged-clamped uniform prismatic beam to a harmonic axial load. The beam is subjected to a static axial load and one of its ends is restrained by a linear spring. We assume that the curvature and inertia nonlinearities are much smaller than the nonlinearity caused by the mid-plane stretching. The transverse deflection w[x, t] of the beam at the position x and time t is governed by the nondimensional integral-partial-differential equation (Nayfeh and Mook, 1979)

EOM = 
      (∂_(t, t))w[x, t] - 2  P (∂_(x, x))w[x, t] + 
          (∂_(x, x, x, x))w[x, t] == 
        -2 ϵ μ[x] (∂_t) w[x, t] + 
          4  ϵ F (∂_(x, x))w[x, t] 
            Cos[Ω t] + 
          ϵ α (∂_(x, x))w[x, t] 
            (∫_0^1)((((∂_x) w[x, t]))^2) 
                \[DifferentialD]x

and homogeneous boundary conditions

BC = {w[0, t] == 0, (∂_(x, x))w[x, t] == 0 /. x -> 0, 
        w[1, t] == 0, (∂_x) w[x, t] == 0 /. x -> 1}

We assume that the static axial load is such that the lowest two natural frequencies of the beam are in the ratio of three-to-one; that is, there is a three-to-one internal resonance between the second and first modes. Moreover, we assume that neither of these two modes is involved in an internal resonance with any other mode. We consider three parametric resonances, namely, principal parametric resonance of the first mode, principal parametric resonance of the second mode, and combination parametric resonance of the first and second modes.

7.2.2 Direct Attack of the Continuous Problem

In this section, we directly attack the integral-partial-differential equation EOM and associated boundary conditions BC and seek a first-order uniform expansion of their solution in the form

solRule = 
      w -> (Evaluate[Sum[(ϵ^j) (Subscript[w, j])[#1, #2, #3], {j, 0, 1}]] &)

where the first independent variable stands for x and the last two independent variables stand for the two time scales (Subscript[T, 0]) and (Subscript[T, 1]). Substituting the solRule into EOM, transforming the total time derivatives into partial derivatives in terms of (Subscript[T, 0]) and (Subscript[T, 1]), expanding the result for small ϵ, discarding terms of order higher than ϵ, and using the intRule2 to simplify the expansions of the integrands, we obtain

eq722a = 
      ((Join[{EOM}, BC] /. Integrate -> int /. {
                        w[x_, t] -> w[x, (Subscript[T, 0]), (Subscript[T, 1])], 
                        Derivative[m_, n_][w][x_, t] -> 
                          dt[n][D[w[x, (Subscript[T, 0]), (Subscript[T, 1])], {x, m}]], t -> (Subscript[T, 0])}
                         /. solRule // ExpandAll) //. intRule2 // 
            ExpandAll) /. (ϵ^(n_ /; n > 1)) -> 0;

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

\(eqEps = 
      \(Thread[CoefficientList[Subtract\ @@\ #, ϵ] == 0]\ &\)\ /@\ 
          eq722a // Transpose; 
eqEps /. displayRule\)

Because in the presence of damping, all modes that are not directly excited by the forcing or indirectly excited by the internal resonance will decay with time (Nayfeh and Mook, 1979), the solution of eqEps[[1]] can be expressed in terms of the lowest two linear free-vibration modes; that is,

\(\(sol0 = 
      w\_0 -> Function[{x, T\_0, T\_1}, 
          Sum[A\_i[T\_1]\ ϕ\_i[x]\ Exp[I\ ω\_i\ T\_0] + 
              A\&_\_i[T\_1]\ ϕ\_i[x]\ Exp[\(-I\)\ ω\_i\ T\_0], {
              i, 2}]]; \)\)

where (Subscript[ω, 1 ])and (Subscript[ω, 2]) are the natural frequencies of these modes. For later use, we define the list

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

One can easily show that the mode shapes (Subscript[ϕ, i])[x] are orthogonal. We assume that these modes are normalized so that

(∫_0^1) (Subscript[ϕ, i])[x] (Subscript[ϕ, j])[x] dx -> 
      (Subscript[δ, ij])

where (Subscript[δ, ij]) is the Kronecker delta function.

Substituting sol0 into the first-order equation, eqEps[[2,1]], and using intRule1 and intRule2, we obtain

\(order1Eq = 
      \((eqEps[\([1, 1, 1]\)] /. w\_0 -> w\_1)\) == 
        \((\(\((eqEps[\([1, 1, 1]\)] /. w\_0 -> w\_1)\) - 
                  \((\(Subtract\ @@\ eqEps[\([2, 1]\)] /. sol0 // TrigToExp
                        \) // Expand)\) /. intRule1\) //. intRule2 // Expand)
          \); 
order1Eq /. displayRule\)

It follows from eqEps[[2]] that the first-order boundary conditions are

order1BC = eqEps[[2]] // Rest

Principal Parametric Resonance of the First Mode

In this case, Ω ≈ 2  (Subscript[ω, 1]). To describe the nearness of the internal and principal parametric resonances, we introduce the two detuning parameters (Subscript[σ, 1]) and (Subscript[σ, 2]) defined by

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

and define the following rules:

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand side of order1Eq and obtain the source of secular terms as

\(ST = 
      Table[Coefficient[order1Eq[\([2]\)] /. expRule[i], 
          Exp[I\ ω\_i\ T\_0]], {i, 2}]; 
ST /. displayRule\)

Because the homogeneous part of order1Eq and order1BC has a nontrivial solution, the corresponding nonhomogeneous problem has a solution only if solvability conditions are satisfied. It follows from Section 7.1.1 that the homogeneous problem is self-adjoint and that the solvability conditions demand that ST be orthogonal to (Subscript[ϕ, 1])[x] and (Subscript[ϕ, 2])[x], respectively. Imposing these conditions, we have

\(SCond = 
      Table[\(\(int[ϕ\_j[x]\ ST[\([j]\)], {x, 0, 1}] /. intRule1\) //. 
                intRule2\) /. int -> Integrate, {j, 2}] == 0 // Thread; 

    SCond /. displayRule\)

To simplify the notation in the solvability conditions, we use the orthonormality of the mode shapes and define the following parameters:

notationRule={(∫_0^1) (Subscript[ϕ, i])_[x]^2 \[DifferentialD]x -> 1,(∫_0^1) μ[x] (Subscript[ϕ, i])_[x]^2 
                  \[DifferentialD]x -> (Subscript[μ, i]),F(∫_0^1)(Subscript[ϕ, i])_[x] (Subscript[ϕ, j_]^″)[x]\[DifferentialD]x->-(Subscript[ω, i]) (Subscript[s, i, j]) f};

To identify the forms of the nonlinear terms in the solvability conditions, we first identify all of the possible forms of the nonhomogeneous terms in the first-order problem. To this end, we let

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

Then, all of the possible forms of the nonlinear terms in the first-order problem are given by

\(cubicTerms = 
      \(Nest[Outer[Times, basicTerms, #]\ &, basicTerms, 2] // Flatten\) // 
        Union; 
cubicTerms /. displayRule\)

Out of these terms, only the terms that may lead to secular terms appear in the solvability conditions, which can be identified according to

\(secularTerms = 
      \(\((\(\(\(E\^\(\(-I\)\ ω\_#\ T\_0\)\) cubicTerms /. 
                    expRule[#]\) /. Exp[_\ T\_0 + _. ] -> 0 // Union\) // 
              Rest)\)\ &\)\ /@\ {1, 2}\)

Next, we define the following parameters:

symbolList = {
        -8 (Subscript[ω, 1]) {(Subscript[γ, 11]), (Subscript[δ, 1]), (Subscript[γ, 12])}, 
        -8 (Subscript[ω, 2]) {(Subscript[δ, 2]), (Subscript[γ, 21]), (Subscript[γ, 22])}}

and express SCond in a more concise form as

\(eqMod = 
      \(Expand[\(-\(1\/ω\_#\)\) 
                \((\((\(SCond[\([#, 1]\)] /. notationRule\) /. 
                        Thread[secularTerms[\([#]\)] -> 0])\) + 
                    symbolList[\([#]\)] . secularTerms[\([#]\)])\)] == 0\ &
          \)\ /@\ {1, 2}\)

where

-symbolList[[1]] -> 
          (-Coefficient[SCond[[1, 1]], #] & /@ 
              secularTerms[[1]]) // Thread // TableForm

-symbolList[[2]] -> 
          (-Coefficient[SCond[[2, 1]], #] & /@ 
              secularTerms[[2]]) // Thread // TableForm

Modulation Equations in Polar Form

The complex-valued solvability conditions can be expressed in real-valued form by introducing the polar transformation

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

into eqMod and obtain

\(expr1 = 
      \(Expand[eqMod[\([#, 1]\)]\ Exp[\(-I\)\ θ\_#[T\_1]] /. ruleA]\ &
          \)\ /@\ {1, 2}\)

Next, we separate the real and imaginary parts of expr1 to obtain the equations governing the modulation of the amplitudes (Subscript[a, i]) and phases (Subscript[θ, i]). To accomplish this, we define the rule

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

Then, the equations governing the amplitudes of motion correspond to the imaginary parts of expr1; that is,

ampEq=Solve[ComplexExpand[Im[#]] == 0 /. realRule & /@ expr1,{(Subscript[a, 1]^′)[(Subscript[T, 1])],(Subscript[a, 2]^′)[(Subscript[T, 1])]}][[1]]/.Rule -> Equal//ExpandAll

Moreover, the equations governing the phases of motion correspond to the real parts of expr1; that is,

phaseEq = 
      Solve[ComplexExpand[Re[#]] == 0 /. realRule & /@ expr1, {
                ((Subscript[θ, 1])^′)[(Subscript[T, 1])], 
                ((Subscript[θ, 2])^′)[(Subscript[T, 1])]}][[1]] /. 
          Rule -> Equal // ExpandAll

These modulation equations are nonautonomous because they depend explicitly on (Subscript[T, 1]).

To determine an autonomous set of modulation equations, we start with expr1 and identify the independent phase arguments in it. To accomplish this, we first identify all possible phase arguments; that is,

expTerms = I Cases[expr1, Exp[arg_] -> arg, Infinity] // Expand

Out of these arguments, only two are independent because

Outer[D, expTerms, {(Subscript[θ, 1])[(Subscript[T, 1])], (Subscript[θ, 2])[(Subscript[T, 1])]}] // RowReduce

Denoting these two independent arguments by (Subscript[γ, 1]) and (Subscript[γ, 2]), we have

gammaList = {2  (Subscript[γ, 1])[(Subscript[T, 1])], (Subscript[γ, 2])[(Subscript[T, 1])]} == 
          expTerms[[{1, 2}]] // Thread

which can be solved for (Subscript[θ, 1]) and (Subscript[θ, 2]) to obtain

thetaRule = 
      Solve[gammaList, {(Subscript[θ, 1])[(Subscript[T, 1])], (Subscript[θ, 2])[(Subscript[T, 1])]}][
          [1]] // ExpandAll

Substituting for the (Subscript[θ, i]) in expr1 yields

expr2 = 
      expr1 /. 
            Table[(Subscript[θ, i]) -> 
                (Evaluate[thetaRule[[i, 2]] /. (Subscript[T, 1]) -> #] &), {i, 2}]
               /. Exp[arg_] :> Exp[Expand[arg]] // Expand

Then, the equations governing the amplitudes of motion correspond to the imaginary parts of expr2; that is,

(ampEq=Solve[ComplexExpand[Im[#]] == 0 /. realRule & /@ expr2
                      ,{(Subscript[a, 1]^′)[(Subscript[T, 1])],(Subscript[a, 2]^′)[(Subscript[T, 1])]}][[1]]/.Rule -> Equal//ExpandAll)/.f_[(Subscript[T, 1])] -> f

Moreover, the equations governing the phases of motion correspond to the real parts of expr2; that is,

(phaseEq = 
          Solve[ComplexExpand[Re[#]] == 0 /. realRule & /@ expr2, {
                    ((Subscript[γ, 1])^′)[(Subscript[T, 1])], 
                    ((Subscript[γ, 2])^′)[(Subscript[T, 1])]}][[1]] /. 
              Rule -> Equal // ExpandAll) /. f_[(Subscript[T, 1])] -> f

The modulation equations are autonomous because they are independent of the independent variable (Subscript[T, 1]).

Modulation Equations in Cartesian Form and Symmetry Property

In performing bifurcation analyses of the dynamics of the system, one might find it more convenient, in some cases, to represent the modulation equations in Cartesian rather than polar form. As a byproduct, one can also obtain the symmetries of the system.To this end, we introduce the Cartesian transformation

\(cartRule = {
        A\_k_ -> 
          \((1\/2\ \((p\_k[#] - I\ q\_k[#])\)\ Exp[I\ λ\_k[#]]&)\), 
        A\&_\_k_ -> 
          \((1\/2\ \((p\_k[#] + I\ q\_k[#])\)\ Exp[\(-I\)\ λ\_k[#]]&)
            \)}; 
realRule = {Im[x_] -> 0, Re[x_] -> x, 
        \(λ\_k_\^′\)[T\_1] -> ν\_k}; \)

where the (Subscript[λ, i]) are chosen to render the resulting modulation equations autonomous. Substituting the cartRule into eqMod yields

eqModCart[k_] := 
      Exp[-I (Subscript[λ, k])[(Subscript[T, 1])]] Subtract @@ eqMod[[k]] /. 
          cartRule // Expand;

Next, we choose the (Subscript[λ, i]) to render eqModCart autonomous. To accomplish this, we identify its independent arguments by identifying first all of its arguments. They are given by

list1 = 
      -I (
            Cases[#, Exp[x_] -> x, Infinity] & /@ 
                  Array[eqModCart, 2] // Flatten // Union) // Expand

To determine the number of independent relations in list1, we use the function RowReduce and obtain

Outer[D, list1, {(Subscript[λ, 1])[(Subscript[T, 1])], (Subscript[λ, 2])[(Subscript[T, 1])]}] // RowReduce

Hence, there are only two linearly independent arguments in list1. We choose the first two to determine the (Subscript[λ, i]). For an invariant set of modulation equations, we set each of these relations to be an even multiple of π. Solving the resulting relations for the (Subscript[λ, i]), we obtain

list1 = list1[[{1, 2}]];

lambdaRule = 
      Solve[list1 + {(Subscript[m, 1]), (Subscript[m, 2])} 2  π == 0 // Thread, 
          Table[(Subscript[λ, i])[(Subscript[T, 1])], {i, 2}]][[1]]

where the (Subscript[m, i]) are positive or negative integers. Letting (Subscript[ν, i]) = ((Subscript[λ, i])^′) yields

D[lambdaRule, (Subscript[T, 1])] /. realRule

Then, separating the real and imaginary parts of eqModCart, we obtain the following Cartesian form of the modulation equations:

(eqs = 
            Solve[Flatten[
                        Table[ComplexExpand[{Im[#], Re[#]}] & /@ 
                              (eqModCart[k] /. Exp[_] -> 1) //. 
                            realRule, {k, 2}]] == 0 // Thread, 
                  Table[{((Subscript[p, k])^′)[(Subscript[T, 1])], 
                        ((Subscript[q, k])^′)[(Subscript[T, 1])]}, {k, 2}] // Flatten][
                [1]] /. Rule -> Equal) /. f_[(Subscript[T, 1])] -> f // TableForm

The symmetry property of the system can also be obtained as follows:

phase = 
      Flatten[Table[
            Mod[#[[2]] & /@ lambdaRule /. (Subscript[T, 1]) -> 0, 2  π], {
              (Subscript[m, 1]), 0, 10}, {(Subscript[m, 2]), 0, 10}], 1] // Union

rotMat[a_] = {{Cos[a], Sin[a]}, {-Sin[a], Cos[a]}};

symmetryList = 
      Join @@ Table[rotMat[#[[i]]] . {(Subscript[p, i]), (Subscript[q, i])}, {i, 2}] & /@ 
        phase

Hence, if {(Subscript[p, 1]), (Subscript[q, 1]), (Subscript[p, 2]), (Subscript[q, 2])} is a solution of the modulation equations, then {-(Subscript[p, 1]), -(Subscript[q, 1]), -(Subscript[p, 2]), -(Subscript[q, 2])} is also another solution of these equations. If these two sets of solutions are the same, then the system response is symmetric; otherwise, it will be asymmetric.

The Functions PolarForm and CartesianForm

Following the procedures described in the preceding section, we build two functions, PolarForm and CartesianForm, to automate the process provided that the complex dependent variables in the modulation equations are written in the form (Subscript[S, i]), i = 1,  ... , neq and the corresponding complex conjugates are \(TraditionalForm\`S\&_\_i\). The functions can be slightly modified to allow a more general form for the modulation equations. To obtain the polar form or Cartesian form of the complex modulation equations, we simply provide the list of modulation equations and S as two required inputs in PolarForm or CartesianForm, respectively. If S = A, the list of modulation equations is the only required input.

\(PolarForm[eqs_List, S_Symbol:  A] := 
	
      Module[{neq = Length[eqs], t}, 
		\ \ \ 
        t = \(Cases[eqs, S\_i_[a_] -> a, ∞]\)[\([1]\)]; 
		\ \ \ 
        polarRule = {
            S\_k_ -> \((\(1\/2\) a\_k[#]\ Exp[I\ β\_k[#]]\ &)\), 
            S\&_\_k_ -> \((\(1\/2\) a\_k[#]\ Exp[\(-I\)\ β\_k[#]]&)\)}; 
        
\ \ \ 		realRule = {Im[x_] -> 0, Re[x_] -> x}; 
\ \ \ 		
        eqModPolar[k_] := 
          Exp[\(-I\)\ β\_k[t]]\ Subtract\ @@\ eqs[\([k]\)] /. 
              polarRule // Expand; 
		\ \ \ 
        list1 = \(I\ 
                Cases[Array[eqModPolar, neq], Exp[a_] -> a, ∞] // 
              Expand\) // Union; 
\ \ \ 		
        newLists[list_, n_] := 
			\ \ \ 
          With[{v = Table[Unique[], {n}]}, 
				\ \ 
            ReplaceList[list, 
              Append[Flatten[\({___, Pattern[#, _]}\ &\)\ /@\ v], ___] -> 
                v]]; 
	\ \ \ 	betaList = Table[β\_i[t], {i, neq}]; 

        \ \ \ 		list2 = {}; 
		\ \ \ 
        Scan[If[NullSpace[Outer[Coefficient, #, betaList]] === {}, 
              list2 = \(-#\); Return[]]\ &, newLists[list1, neq]]; 
		\ \ 
        If[list2 === {}, 
          Print["Autonomous system is not possible!!!
The phase list is:"]; Return[list1]
		\ \ ]; 
\ \ \ \ \ \ \ 
        betaRule1 = 
          \(Solve[list3 = Table[γ\_i[t], {i, neq}] == list2 // Thread, 
              betaList]\)[\([1]\)]; 
\ \ 		
        betaRule2 = 
          Table[β\_i -> 
              Function[t // Evaluate, 
                β\_i[t] /. betaRule1 // Evaluate], {i, neq}]; 
\ \ 		
        Do[eq0[i] = 
            \(eqModPolar[i] /. betaRule2\) /. Exp[a_] :> Exp[Expand[a]]; 

          \ \ \ \ \ \ \ \ \ 
          eq[i] = \(Solve[ComplexExpand[Im[eq0[i]]] == 0 /. realRule, 
                  \(a\_i\^′\)[t]]\)[\([1, 1]\)] /. Rule -> Equal; 

          \ \ \ \ \ \ \ \ \ 
          eqRe[i] = ComplexExpand[Re[eq0[i]]] == 0 /. realRule, {i, neq}

          		\ \ ]; 
\ \ \ \ \ \ \ 
        gammapList = Table[\(γ\_i\^′\)[t], {i, neq}]; 
		\ \ 
        eqList = 
          \(Solve[Array[eqRe, neq], gammapList]\)[\([1]\)] /. 
              Rule -> Equal // ExpandAll; 
\ \ \ \ \ \ \ 
        Do[eq[i + neq] = eqList[\([i]\)], {i, neq}]; 

        \ \ \ \ \ \ \ {Array[eq, 2*neq], list3}
\ ]\)

\(CartesianForm[eqs_List, S_Symbol:  A] := 
	
      Module[{neq = Length[eqs], t}, 
\ \ 		\ 
        t = \(Cases[eqs, S\_i_[a_] -> a, ∞]\)[\([1]\)]; 
		\ \ \ 
        cartRule = {
            S\_k_ -> 
              \((1/2 \((p\_k[#] - I\ q\_k[#])\)\ Exp[I\ λ\_k[#]]\ &)
                \), S\&_\_k_ -> 
              \((1/2 \((p\_k[#] + I\ q\_k[#])\)\ 
                    Exp[\(-I\)\ λ\_k[#]]\ &)\)}; 
\ \ \ 		
        realRule = 
          \ {Im[x_] -> 0, Re[x_] -> x, 
            \(λ\_k_\^′\)[t] -> ν\_k}; 
\ 	\ \ 	
        eqModCart[k_] := 
          Exp[\(-I\)\ λ\_k[t]]\ Subtract\ @@\ eqs[\([k]\)] /. 
              cartRule // Expand; 
\ \ \ 		
        eqs1 = \({ComplexExpand[Im[#]] == 0, ComplexExpand[Re[#]] == 0}\ &
                  \)\ /@\ \((Array[eqModCart, neq] /. E^_ -> 1)\) //. 
              realRule // Flatten; 
\ \ \ 	\ 
        eqs2 = \(Solve[eqs1, 
                Table[{\(p\_k\^′\)[t], \(q\_k\^′\)[t]}, {k, 
                      neq}] // Flatten]\)[\([1]\)] /. Rule -> Equal; 		

        		\ \ list1 = 
          \(I\ Cases[Array[eqModCart, neq], Exp[a_] -> a, ∞] // 
              Expand\) // Union; 
		\ \ 
        newLists[list_, n_] := 
			\ \ \ 
          With[{v = Table[Unique[], {n}]}, 
				\ \ 
            ReplaceList[list, 
              Append[Flatten[\({___, Pattern[#, _]}\ &\)\ /@\ v], ___] -> 
                v]]; 
	\ \ \ \ 
        lambdaList = Table[λ\_i[t], {i, neq}]; 
\ \ \ 	\ 
        list2 = {}; 
		\ \ 
        Scan[If[NullSpace[Outer[Coefficient, #, lambdaList]] === {}, 
              list2 = \(-#\); Return[]]\ &, newLists[list1, neq]]; 
		\ \ 
        If[list2 === {}, 
          Print["Autonomous system is not possible!!!
The phase list is:"]; Return[list1]
		\ \ ]; 
\ \ 		
        lambdaRule = 
          \(Solve[list2 + 2\ Table[m\_i, {i, neq}]\ π == 0 // Thread, 
              lambdaList]\)[\([1]\)]; 
		\ \ 
        iterList = Table[{m\_i, 0, 10}, {i, neq}]; 
	\ \ 	
        phase = Flatten[
              Table[Mod[
                    \(Select[#, \(! FreeQ[#, π]\)\ &]\ &\)\ /@\ 
                      \((lambdaList /. lambdaRule // Expand)\), 2  π] // 
                  Evaluate, Sequence\ @@\ iterList // Evaluate], neq - 1] // 
            Union; 
\ \ 		
        rotMat[θ_] := {{Cos[θ], \(-Sin[θ]\)}, {
              Sin[θ], Cos[θ]}}; 
		\ \ 
        symmetryList = 
          \(Flatten[Table[rotMat[#[\([k]\)]] . {p\_k, q\_k}, {k, 2}]]\ &\)\ /@
            \ phase; 

        \ 	\ \ \ {eqs2, D[lambdaRule, t] /. realRule, symmetryList}
\ ]\)

Principal Parametric Resonance of the Second Mode

In this case, Ω ≈ (Subscript[ω, 2]). To describe the nearness of the internal and principal parametric resonances, we introduce the two detuning parameters (Subscript[σ, 1]) and (Subscript[σ, 2]) defined by

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

and define the following rules:

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand side of order1Eq and obtain the source of secular terms as

\(ST = 
      Table[Coefficient[order1Eq[\([2]\)] /. expRule[i], 
          Exp[I\ ω\_i\ T\_0]], {i, 2}]; 
ST /. displayRule\)

As discussed in the case of principal parametric resonance of the first mode, the solvability conditions of order1Eq and order1BC demand that ST be orthogonal to (Subscript[ϕ, 1])[x] and (Subscript[ϕ, 2])[x], respectively. Imposing these conditions, we have

\(SCond = 
      Table[\(\(int[ϕ\_j[x]\ ST[\([j]\)], {x, 0, 1}] /. intRule1\) //. 
                intRule2\) /. int -> Integrate, {j, 2}] == 0 // Thread; 

    SCond /. displayRule\)

Using the notations introduced in the preceding section, we can express SCond in a more concise form as

\(eqMod = 
      \(Expand[\(-\(1\/ω\_#\)\) 
                \((\((\(SCond[\([#, 1]\)] /. notationRule\) /. 
                        Thread[secularTerms[\([#]\)] -> 0])\) + 
                    symbolList[\([#]\)] . secularTerms[\([#]\)])\)] == 0\ &
          \)\ /@\ {1, 2}\)

where

-symbolList[[1]] -> 
          (-Coefficient[SCond[[1, 1]], #] & /@ 
              secularTerms[[1]]) // Thread // TableForm

-symbolList[[2]] -> 
          (-Coefficient[SCond[[2, 1]], #] & /@ 
              secularTerms[[2]]) // Thread // TableForm

Modulation Equations in Polar Form

Using eqMod and the function PolarForm defined in the preceding section, we obtain the modulation equations in polar form and the definitions for (Subscript[γ, i]) as

PolarForm[eqMod]

Modulation Equations in Cartesian Form and Symmetry Property

Using eqMod and the function CartesianForm defined in the preceding section, we obtain the modulation equations in Cartesian form, the definitions for (Subscript[ν, i]), and the symmetry property as

CartesianForm[eqMod]

Combination Parametric Resonance of the Two Modes

In this case, \(TraditionalForm
      \`Ω ≈ ω\_1 + ω\_2\). To describe the nearness of the internal and combination parametric resonances, we introduce the two detuning parameters (Subscript[σ, 1]) and (Subscript[σ, 2]) defined by

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

and define the following rules:

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand side of order1Eq and obtain the source of secular terms as

\(ST = 
      Table[Coefficient[order1Eq[\([2]\)] /. expRule[i], 
          Exp[I\ ω\_i\ T\_0]], {i, 2}]; 
ST /. displayRule\)

The solvability conditions demand that ST be orthogonal to solutions of the adjoint homogeneous problem. Since the problem is self-adjoint, we have

\(SCond = 
      Table[\(\(int[ϕ\_j[x]\ ST[\([j]\)], {x, 0, 1}] /. intRule1\) //. 
                intRule2\) /. int -> Integrate, {j, 2}] == 0 // Thread; 

    SCond /. displayRule\)

Using the notations introduced in the preceding section, we can express SCond in a more concise form as

\(eqMod = 
      \(Expand[\(-\(1\/ω\_#\)\) 
                \((\((\(SCond[\([#, 1]\)] /. notationRule\) /. 
                        Thread[secularTerms[\([#]\)] -> 0])\) + 
                    symbolList[\([#]\)] . secularTerms[\([#]\)])\)] == 0\ &
          \)\ /@\ {1, 2}\)

where

-symbolList[[1]] -> 
          (-Coefficient[SCond[[1, 1]], #] & /@ 
              secularTerms[[1]]) // Thread // TableForm

-symbolList[[2]] -> 
          (-Coefficient[SCond[[2, 1]], #] & /@ 
              secularTerms[[2]]) // Thread // TableForm

Modulation Equations in Polar Form

Using eqMod and the function PolarForm, we obtain the modulation equations in polar form and the definitions for (Subscript[γ, i]) as

PolarForm[eqMod]

Modulation Equations in Cartesian Form and Symmetry Property

Using eqMod and the function CartesianForm, we obtain the modulation equations in Cartesian form, the definitions for (Subscript[ν, i]), and the symmetry property as

CartesianForm[eqMod]

7.2.3 Discretization of the Continuous Problem

As an alternative, we apply the method of multiple scales to the discretized system of EOM and BC. To determine the discretized form, we expand w[x, t] in terms of the linear mode shapes (Subscript[ϕ, m])[x] as

wRule1 = 
      w -> Function[{x, t}, 
          (∑_(m = 1)^∞) (Subscript[u, m])[t] (Subscript[ϕ, m])[x]];

To simplify the computation, we temporarily drop the ∑ sign in wRule1, as long as we know the repeated m represents a summation index, and rewrite it as

wRule2 = w -> Function[{x, t}, (Subscript[u, m])[t] (Subscript[ϕ, m])[x]];

where the (Subscript[u, m])[t] are the generalized coordinates. Substituting wRule2 into EOM, multiplying the result with (Subscript[ϕ, n])[x], and integrating the outcome from x = 0 to x = 1 using intRule1 and intRule2, we obtain the following discretized form of the equations describing the response of the beam:

eq723a = 
      int[(Subscript[ϕ, n])[x] (# /. Integrate -> int /. wRule2), {x, 
                    0, 1}] /. intRule1 //. intRule2 /. int -> Integrate &
           /@ EOM

where the damping is assumed to be modal. The nonlinear term, in general, should be rewritten as

NT=(Coefficient[eq723a[[2]], (Subscript[u, m])[t]^3]/.(((Subscript[ϕ, m]^′)[x])^2)->((Subscript[ϕ, k])^′)[x] 
                ((Subscript[ϕ, ℓ])^′)[x]) (Subscript[u, m])[t] (Subscript[u, k])[t] (Subscript[u, ℓ])[t]

Using the orthonormality of the mode shapes (Subscript[ϕ, i])[x] and the results from the corresponding eigenvalue problem, we define the following rules:

notationRule1={(∫_0^1) (Subscript[ϕ, m])[x] (Subscript[ϕ, n])[x] \[DifferentialD]x -> 
              (Subscript[δ, m, n]),(∫_0^1)(Subscript[ϕ, n])[x] (Subscript[ϕ, m]^((4)))[x]\[DifferentialD]x->2P (∫_0^1)(Subscript[ϕ, n])[x] (Subscript[ϕ, m]^″)[x]\[DifferentialD]x+(Subscript[ω, m]^2) (Subscript[δ, m, n])};

notationRule2={(∫_0^1) μ[x] (Subscript[ϕ, m])[x] (Subscript[ϕ, n])[x] 
                  \[DifferentialD]x -> (Subscript[μ, m]) (Subscript[δ, m, n]),(∫_0^1)(Subscript[ϕ, n])[x] (Subscript[ϕ, m]^″)[x]\[DifferentialD]x->-(Subscript[γ, n, m]),((∫_0^1)(Subscript[ϕ, k]^′)[x] (Subscript[ϕ, ℓ]^′)[x]\[DifferentialD]x) ((∫_0^1)(Subscript[ϕ, n])[x] (Subscript[ϕ, m]^″)[x]\[DifferentialD]x)->-(Subscript[γ, k, ℓ]) (Subscript[γ, n, m]),a_ (Subscript[δ, m, n]) :> (a /. m -> n)};

and then rewrite eq723a as

EOM1 = 
      Expand[eq723a[[1]] /. notationRule1] == 
          (eq723a[[2]] /. (Subscript[u, m])[t]^3 -> 0) + NT //. notationRule2

where m, k,  and ℓ represent summation indicies.

Using the method of multiple scales, we seek a first-order uniform expansion in the form

solRule = 
      (Subscript[u, n])_ -> ((∑_(j = 0)^1) (ϵ^j) (Subscript[u, n, j])[#1, #2] &)
          ;

Transforming the total time derivatives in EOM1 into partial derivatives in terms of (Subscript[T, 0]) and (Subscript[T, 1]), substituting the solRule into EOM1, expanding the result for small ϵ, and discarding terms of order higher than ϵ, we obtain

eq723b = 
      (EOM1 /. {(Subscript[u, s])_[t] -> (Subscript[u, s])[(Subscript[T, 0]), (Subscript[T, 1])], 
                  Derivative[n_][(Subscript[u, s])_][t] -> 
                    dt[n][(Subscript[u, s])[(Subscript[T, 0]), (Subscript[T, 1])]], t -> (Subscript[T, 0])} /. solRule // 
            ExpandAll) /. (ϵ^(n_ /; n > 1)) -> 0

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

eqEps = CoefficientList[#, ϵ] & /@ eq723b // Thread

Because in the presence of damping, all modes that are not directly or indirectly excited decay with time, all of the (Subscript[u, n, 0]) decay with time except (Subscript[u, 1, 0]) and (Subscript[u, 2, 0]). Then it follows from eqEps[[1]] that

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

Taking account of the summation indicies, we define

sumRule = {
        (Subscript[γ, k, ℓ]) (Subscript[γ, n, m]) (Subscript[u, k, 0])[a__] 
            (Subscript[u, m, 0])[a__] (Subscript[u, ℓ, 0])[a__] -> 
          (∑_(k = 1)^2)
              (∑_(m = 1)^2)
                (∑_(ℓ = 1)^2) (Subscript[γ, k, ℓ]) 
                  (Subscript[γ, n, m]) (Subscript[u, k, 0])[a] (Subscript[u, m, 0])[a] 
                  (Subscript[u, ℓ, 0])[a], 
        (Subscript[γ, n, m]) (Subscript[u, m, 0])[a__] -> 
          (∑_(m = 1)^2) (Subscript[γ, n, m]) (Subscript[u, m, 0])[a]};

Hence, we rewrite eqEps[[2]] as

Clear[order1Eq]

order1Eq[n_] = 
      (eqEps[[1, 1]] /. (Subscript[u, n, 0]) -> (Subscript[u, n, 1])) == 
          (eqEps[[1, 1]] /. (Subscript[u, n, 0]) -> (Subscript[u, n, 1])) - 
            (Subtract @@ eqEps[[2]]) /. sumRule

Substituting sol0 into order1Eq[1] and order1Eq[2] yields

eqOrder1 = Array[order1Eq, 2] /. sol0 // TrigToExp // ExpandAll;

Principal Parametric Resonance of the First Mode

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

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand sides of eqOrder1 and obtain the solvability condition as

\(SCond = 
      Table[Coefficient[\(-eqOrder1[\([i, 2]\)]\) /. expRule[i], 
            Exp[I\ ω\_i\ T\_0]] == 0, {i, 2}]; 

    SCond /. displayRule\)

which is in agreement with that obtained by direct approach.

Principal Parametric Resonance of the Second Mode

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

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand sides of eqOrder1 and obtain the solvability condition as

\(SCond = 
      Table[Coefficient[\(-eqOrder1[\([i, 2]\)]\) /. expRule[i], 
            Exp[I\ ω\_i\ T\_0]] == 0, {i, 2}]; 

    SCond /. displayRule\)

which is in agreement with that obtained by direct approach.

Combination Parametric Resonance of the Two Modes

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

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand sides of eqOrder1 and obtain the solvability condition as

\(SCond = 
      Table[Coefficient[\(-eqOrder1[\([i, 2]\)]\) /. expRule[i], 
            Exp[I\ ω\_i\ T\_0]] == 0, {i, 2}]; 

    SCond /. displayRule\)

which is in agreement with that obtained by direct approach.

7.2.4 Method of Time-Averaged Lagrangian

As a second alternative, we derive the modulation equations by using the method of time-averaged Lagrangian. The nondimensional Lagrangian of the beam can be expressed as (Nayfeh, 1998)

Lagr1 = 
      (1/2) (∫_0^1)((((∂_t) w[x, t]))^2) 
              \[DifferentialD]x - 
        (1/2) 
          (∫_0^1)((((∂_(x, x))w[x, t]))^2) 
              \[DifferentialD]x - 
        (P + 2  ϵ F Cos[Ω t]) 
          (∫_0^1)((((∂_x) w[x, t]))^2) 
              \[DifferentialD]x - 
        (1/4) ϵ α 
          ((((∫_0^1)((((∂_x) w[x, t]))^2) 
                  \[DifferentialD]x))^2) + C[t];

where C[t] is independent of w. Transforming the total time derivatives into partial derivatives in terms of (Subscript[T, 0]) and (Subscript[T, 1]), we modify Lagr1 as

Lagr2 = 
      Lagr1 /. Integrate -> int /. {w[x_, t] -> w[x, (Subscript[T, 0]), (Subscript[T, 1])], 
          Derivative[m_, n_][w][x_, t] -> 
            dt[n][D[w[x, (Subscript[T, 0]), (Subscript[T, 1])], {x, m}]], t -> (Subscript[T, 0])}

To this end, we let

\(solRule = 
      w -> \((Evaluate[
              Sum[A\_i[#3]\ ϕ\_i[#1]\ E\^\(I\ ω\_i\ #2\) + 
                  A\&_\_i[#3]\ ϕ\_i[#1]\ 
                    E\^\(\(-I\)\ ω\_i\ #2\), {i, 2}]]\ &)\)\)

where the first independent variable (#1) stands for x and the last two independent variables (#2 and #3) stand for the two time scales (Subscript[T, 0]) and (Subscript[T, 1]). Substituting solRule into Lagr2, using intRule1 and intRule2 to simplify the expansions of the integrands, and collecting the coefficients of ϵ, we obtain

order1Lagr = 
      Lagr2 /. solRule /. intRule1 //. intRule2 // TrigToExp // 
          Expand // Coefficient[#, ϵ] &;

Principal Parametric Resonance of the First Mode

ResonanceConds = {
          (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
          Ω == 2  (Subscript[ω, 1]) + ϵ (Subscript[σ, 2])} /. 
        Equal -> Rule

Using ResonanceConds, we obtain the slowly varying terms from order1Lagr as

TAL = 
      order1Lagr /. 
            Exp[a_] :> 
              Exp[Expand[a /. ResonanceConds] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]
                 /. Exp[_ (Subscript[T, 0]) + _. ] -> 0 /. int -> Integrate

The Euler-Lagrange equations corresponding to the TAL can be written as

eqMod1=D[D[TAL,(Subscript[\(A\&_\), 1]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_1[T\_1]]\)==0

and

eqMod2=D[D[TAL,(Subscript[\(A\&_\), 2]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_2[T\_1]]\)==0

Adding linear viscous damping to eqMod1 and eqMod2 and performing integration by parts yields the same modulation equations as obtained by directly attacking the partial-differential system.

Principal Parametric Resonance of the Second Mode

ResonanceConds = 
      Solve[{(Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
            Ω == 2  (Subscript[ω, 2]) + ϵ (Subscript[σ, 2])}, {
            (Subscript[ω, 2]), Ω}][[1]]

Using ResonanceConds, we obtain the slowly varying terms from order1Lagr as

TAL = 
      order1Lagr /. 
            Exp[a_] :> 
              Exp[Expand[a /. ResonanceConds] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]
                 /. Exp[_ (Subscript[T, 0]) + _. ] -> 0 /. int -> Integrate

The Euler-Lagrange equations corresponding to the TAL can be written as

eqMod1=D[D[TAL,(Subscript[\(A\&_\), 1]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_1[T\_1]]\)==0

and

eqMod2=D[D[TAL,(Subscript[\(A\&_\), 2]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_2[T\_1]]\)==0

Adding linear viscous damping to eqMod1 and eqMod2 and performing integration by parts yields the same modulation equations as obtained by directly attacking the partial-differential system.

Combination Parametric Resonance of the Two Modes

ResonanceConds = 
      Solve[{(Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + ϵ (Subscript[σ, 1]), 
            Ω == 
              (Subscript[ω, 1]) + (Subscript[ω, 2]) + ϵ (Subscript[σ, 2])}, {
            (Subscript[ω, 2]), Ω}][[1]]

Using ResonanceConds, we obtain the slowly varying terms from order1Lagr as

TAL = 
      order1Lagr /. 
            Exp[a_] :> 
              Exp[Expand[a /. ResonanceConds] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]
                 /. Exp[_ (Subscript[T, 0]) + _. ] -> 0 /. int -> Integrate

The Euler-Lagrange equations corresponding to the TAL can be written as

eqMod1=D[D[TAL,(Subscript[\(A\&_\), 1]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_1[T\_1]]\)==0

and

eqMod2=D[D[TAL,(Subscript[\(A\&_\), 2]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_2[T\_1]]\)==0

Adding linear viscous damping to eqMod1 and eqMod2 and performing integration by parts yields the same modulation equations as obtained by directly attacking the partial-differential system.

7.3 Cantilever Beam

7.3.1 EOM and BC's

We consider the nonlinear nonplanar response of an inextensional cantilever beam to an external excitation of one of its flexural modes. The lowest torsional frequencies of the beam considered are much higher than the frequencies of the excited flexural modes so that the torsional inertia can be neglected. We assume that the beam is uniform and homogeneous. The transverse deflections v[x, t] and w[x, t] of the beam at the position x and time t are governed by the nondimensional integral-partial-differential equation (Crespo da Silva and Glynn, 1978)

EOM1 = {
        (∂_(t, t))v[x, t] + 
            (Subscript[β, y]) (∂_(x, x, x, x))v[x, t] == 
          -2 ϵ μ[x] (∂_t) v[x, t] + 
            ϵ (1 - (Subscript[β, y])) 
              (∂_x)(
                  (∂_(x, x))w[x, t] 
                      (∫_1^x)(∂_(x, x))v[x, t] 
                          (∂_(x, x))w[x, t] \[DifferentialD]x - 
                    (∂_(x, x, x))w[x, t] 
                      (∫_0^x)(∂_(x, x))v[x, t] 
                          (∂_x) w[x, t] \[DifferentialD]x) - 
            ϵ ((((1 - (Subscript[β, y])))^2)/(Subscript[β, γ])) 
              (∂_(x, x))(
                  (∂_(x, x))w[x, t] 
                    (∫_0^x)
                        (∫_1^x)(∂_(x, x))v[x, t] 
                          (∂_(x, x))w[x, t] \[DifferentialD]x 
                          \[DifferentialD]x) - 
            ϵ (Subscript[β, y]) 
              (∂_x)(
                  (∂_x) v[x, t] 
                    (∂_x)(
                        (∂_x) v[x, t] 
                            (∂_(x, x))v[x, t] + 
                          (∂_x) w[x, t] 
                            (∂_(x, x))w[x, t])) - 
            ϵ (1/2) 
              (∂_x)(
                  (∂_x) v[x, t] 
                    (∫_1^x)
                        (∂_(t, t))(
                            (∫_0^x)
                                (((((∂_x) v[x, t]))^2) + 
                                    ((((∂_x) w[x, t]))^2)) 
                                \[DifferentialD]x) \[DifferentialD]x) - 
            ϵ F[x] Cos[Ω t], 
        (∂_(t, t))w[x, t] + (∂_(x, x, x, x))w[x, t] == 
          -2 ϵ μ[x] (∂_t) w[x, t] - 
            ϵ (1 - (Subscript[β, y])) 
              (∂_x)(
                  (∂_(x, x))v[x, t] 
                      (∫_1^x)(∂_(x, x))v[x, t] 
                          (∂_(x, x))w[x, t] \[DifferentialD]x - 
                    (∂_(x, x, x))v[x, t] 
                      (∫_0^x)(∂_(x, x))w[x, t] 
                          (∂_x) v[x, t] \[DifferentialD]x) - 
            ϵ ((((1 - (Subscript[β, y])))^2)/(Subscript[β, γ])) 
              (∂_(x, x))(
                  (∂_(x, x))v[x, t] 
                    (∫_0^x)
                        (∫_1^x)(∂_(x, x))v[x, t] 
                          (∂_(x, x))w[x, t] \[DifferentialD]x 
                          \[DifferentialD]x) - 
            ϵ 
              (∂_x)(
                  (∂_x) w[x, t] 
                    (∂_x)(
                        (∂_x) v[x, t] 
                            (∂_(x, x))v[x, t] + 
                          (∂_x) w[x, t] 
                            (∂_(x, x))w[x, t])) - 
            ϵ (1/2) 
              (∂_x)(
                  (∂_x) w[x, t] 
                    (∫_1^x)
                        (∂_(t, t))(
                            (∫_0^x)
                                (((((∂_x) v[x, t]))^2) + 
                                    ((((∂_x) w[x, t]))^2)) 
                                \[DifferentialD]x) \[DifferentialD]x)};

along with the homogeneous boundary conditions

BC1 = {v[0, t] == 0, (∂_x) v[x, t] == 0 /. x -> 0, 
        (∂_(x, x))v[x, t] == 0 /. x -> 1, 
        (∂_(x, x, x))v[x, t] == 0 /. x -> 1}

BC = Join[BC1, BC1 /. v -> w]

The spatial derivatives outside the integrals on the right-hand sides of EOM1 result in lengthy expressions. To simplify the symbolic computation, we first define some operators as

op1[v_, w_] = 
      (∂_(x, x))w[x, t] 
            (∫_1^x)(∂_(x, x))v[x, t] 
                (∂_(x, x))w[x, t] \[DifferentialD]x - 
          (∂_(x, x, x))w[x, t] 
            (∫_0^x)(∂_(x, x))v[x, t] 
                (∂_x) w[x, t] \[DifferentialD]x /. 
        Integrate -> int[1]

op2[v_, w_] = 
      Hold[(∂_(x, x))w[x, t] 
              (∫_0^x)
                  (∫_1^x)(∂_(x, x))v[x, t] 
                    (∂_(x, x))w[x, t] \[DifferentialD]x 
                    \[DifferentialD]x] /. Integrate -> int[2] // 
        ReleaseHold

op3[v_, w_] = 
      (∂_x) v[x, t] 
          (∫_1^x)
              (∂_(t, t))(
                  (∫_0^x)
                      (((((∂_x) v[x, t]))^2) + 
                          ((((∂_x) w[x, t]))^2)) 
                      \[DifferentialD]x) \[DifferentialD]x /. 
        Integrate -> int[1]

where int[i] indicates the whole expression, including the integrand and the term multiplied by int[i] [expr], will later be wrapped by the ith spatial derivative in the solvability condition.

Then we rewrite EOM1 as

EOM = {
        (∂_(t, t))v[x, t] + 
            (Subscript[β, y]) (∂_(x, x, x, x))v[x, t] == 
          -2 ϵ μ[x] (∂_t) v[x, t] + 
            ϵ (1 - (Subscript[β, y])) op1[v, w] - 
            ϵ ((((1 - (Subscript[β, y])))^2)/(Subscript[β, γ])) 
              op2[v, w] - 
            ϵ (Subscript[β, y]) 
              (∂_x)(
                  (∂_x) v[x, t] 
                    (∂_x)(
                        (∂_x) v[x, t] 
                            (∂_(x, x))v[x, t] + 
                          (∂_x) w[x, t] 
                            (∂_(x, x))w[x, t])) - 
            ϵ (1/2) op3[v, w] - 
            ϵ F[x] Cos[Ω t], 
        (∂_(t, t))w[x, t] + (∂_(x, x, x, x))w[x, t] == 
          -2 ϵ μ[x] (∂_t) w[x, t] - 
            ϵ (1 - (Subscript[β, y])) op1[w, v] - 
            ϵ ((((1 - (Subscript[β, y])))^2)/(Subscript[β, γ])) 
              op2[w, v] - 
            ϵ 
              (∂_x)(
                  (∂_x) w[x, t] 
                    (∂_x)(
                        (∂_x) v[x, t] 
                            (∂_(x, x))v[x, t] + 
                          (∂_x) w[x, t] 
                            (∂_(x, x))w[x, t])) - 
            ϵ (1/2) op3[w, v]};

It follows from EOM1 and BC that the linear undamped natural frequencies and their corresponding mode shapes are given by

\(\({ω\&^\_\(1, m\) == λ\_m\%2\ \@β\_y, 
      ω\_\(2, n\) == λ\_n\%2, 

      ϕ\_i[x] == 
        Cosh[λ\_i\ x] - Cos[λ\_i\ x] + 
          \(Cos[λ\_i] + Cosh[λ\_i]\)\/\(Sin[λ\_i] + 
                Sinh[λ\_i]\)\ 
            \((Sin[λ\_i\ x] - Sinh[λ\_i\ x])\)\  (*\ i = m, n
        \ *) }; \)\)

where \(TraditionalForm\`ω\&^\_\(1, \ m\)\) and (Subscript[ω, 2,  n]) are the natural frequencies in the y and z directions and the (Subscript[λ, i]) are the roots of

1 + Cos[(Subscript[λ, i])] Cosh[(Subscript[λ, i])] == 0

We consider the case of one-to-one internal resonance between the mth mode in the y or v direction and the nth mode in the z or w direction; that is, \(TraditionalForm
      \`ω\&^\_\(1, \ m\) ≈ ω\_\(2, \ n\)\). To express the nearness of these frequencies quantitatively, we let

betayRule = {
        (Subscript[β, y]) -> 1 + (Subscript[δ, 0]) + ϵ (Subscript[δ, 1])};

so that \(TraditionalForm
      \`ω\_\(1, \ m\) = 
        \(\(λ\_m\%2\) \@\(1 + δ\_0\) = 
          \(λ\_n\%2 = ω\_\(2, \ n\)\)\)\). The beam has a near-square cross-section when (Subscript[δ, 0]) = 0. We assume neither of these two modes is involved in an internal resonance with any other mode. Moreover, we consider a primary resonance of the flexural mode in the y direction.

To manipulate some complicated integrals, we define the following rules:

\(intRule3 = {
        \(int[i_]\)[fun_, arg2__] :> \(int[i]\)[Expand[fun], arg2]}; 

    intRule4 = {
        \(int[i_]\)[a_ + b_, arg2__] :> 
          \(int[i]\)[a, arg2] + \(int[i]\)[b, arg2], 
        \(int[i_]\)[ϵ\^n_. \ fun_, arg2__] :> 
          ϵ\^n\ \(int[i]\)[fun, arg2], 
        \(int[i_]\)[a_\ fun_, a1_, b1___] :> 
          a\ \(int[i]\)[fun, a1, b1] /; FreeQ[a, First[a1]], 
        \(int[i_]\)[\(int[i_]\)[a1_, a2_]\ fun_. , a3__] :> 
          \(int[i]\)[a1, a2]\ \(int[i]\)[fun, a3] /; 
            FreeQ[Rest[a2], First[a2]], 
        \(int[i_]\)[\(int[i_]\)[a1_, a2_], a3__] :> \(int[i]\)[a1, a3, a2]}; 
    \)

7.3.2 Direct Attack of the Continuous Problem

In this section, we directly attack the integral-partial-differential equation EOM and associated boundary conditions BC and seek a first-order uniform expansion of their solution in the form

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

where the first independent variable stands for x and the last two independent variables stand for the two time scales (Subscript[T, 0]) and (Subscript[T, 1]). Transforming the total time derivatives in EOM into partial derivatives in terms of (Subscript[T, 0]) and (Subscript[T, 1]), substituting the solRule and betayRule into EOM and BC, expanding the result for small ϵ, discarding terms of order higher than ϵ, and using intRule3 and intRule4 to simplify the expansions of the integrands, we obtain

(eq732a = 
          ((Join[EOM, BC] /. {v[x_, t] -> v[x, (Subscript[T, 0]), (Subscript[T, 1])], 
                                w[x_, t] -> w[x, (Subscript[T, 0]), (Subscript[T, 1])], 
                                Derivative[m_, n_][w_][x_, t] -> 
                                  dt[n][D[w[x, (Subscript[T, 0]), (Subscript[T, 1])], {x, m}]], 
                                t -> (Subscript[T, 0])} /. solRule /. betayRule // 
                        ExpandAll) /. intRule3 //. intRule4 // ExpandAll)
               /. (ϵ^(n_ /; n > 1)) -> 0; ) // Timing

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

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

eqEps[[1]] /. displayRule

In the presence of damping, all modes that are not directly excited by the forcing or indirectly excited by the internal resonance will decay with time. Hence, the solution of eqEps[[1]] can be expressed in terms of the two excited linear free-vibration modes; that is,

\(\(sol0 = {
        v\_0 -> Function[{x, T\_0, T\_1}, 
            A\_1[T\_1]\ ϕ\_m[x]\ Exp[I\ ω\_\(1, m\)\ T\_0] + 
              A\&_\_1[T\_1]\ ϕ\_m[x]\ 
                Exp[\(-I\)\ ω\_\(1, m\)\ T\_0]], 
        w\_0 -> Function[{x, T\_0, T\_1}, 
            A\_2[T\_1]\ ϕ\_n[x]\ Exp[I\ ω\_\(2, n\)\ T\_0] + 
              A\&_\_2[T\_1]\ ϕ\_n[x]\ 
                Exp[\(-I\)\ ω\_\(2, n\)\ T\_0]]}; \)\)

One can easily show that the mode shapes (Subscript[ϕ, i])[x] are orthogonal. We assume that these modes are normalized so that

(∫_0^1) (Subscript[ϕ, i])[x] (Subscript[ϕ, j])[x] dx -> 
      (Subscript[δ, ij])

where (Subscript[δ, ij]) is the Kronecker delta function.

Substituting sol0 into the first-order equations, eqEps[[2,1]] and eqEps[[2,2]], and using intRule3 and intRule4, we obtain

Clear[order1Eq]

order1Eq[1] = 
      (eqEps[[1, 1, 1]] /. (Subscript[v, 0]) -> (Subscript[v, 1])) == 
        ((eqEps[[1, 1, 1]] /. (Subscript[v, 0]) -> (Subscript[v, 1])) - 
                  (Subtract @@ eqEps[[2, 1]] /. sol0 // TrigToExp
                         // Expand) /. intRule3 //. intRule4 // Expand)
          ;

order1Eq[2] = 
      (eqEps[[1, 2, 1]] /. (Subscript[w, 0]) -> (Subscript[w, 1])) == 
        ((eqEps[[1, 2, 1]] /. (Subscript[w, 0]) -> (Subscript[w, 1])) - 
                  (Subtract @@ eqEps[[2, 2]] /. sol0 // TrigToExp
                         // Expand) /. intRule3 //. intRule4 // Expand)
          ;

It follows from eqEps[[2]] that the first-order boundary conditions are

order1BC = Drop[eqEps[[2]], 2]

We consider the case of primary resonance of the flexural mode in the y direction, \(TraditionalForm
      \`Ω ≈ ω\_\(1, \ m\)\). To describe the nearness of the primary resonance, we introduce the detuning parameter σ by

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

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

and define the following rules:

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand side of order1Eq[1] and order1Eq[2] and obtain the source of secular terms as

ST = 
      Table[Coefficient[order1Eq[i][[2]] /. expRule[i], 
          Exp[I omgList[[i]] (Subscript[T, 0])]], {i, 2}];

The solvability conditions of order1Eq[i] and order1BC demand that ST[[i]] be orthogonal to (Subscript[ϕ, i])[x]. Imposing these conditions, we have

SCond1 = 
      int[(Subscript[ϕ, m])[x] ST[[1]], {x, 0, 1}] == 0 /. intRule1 //. 
        intRule2;

SCond2 = 
      int[(Subscript[ϕ, n])[x] ST[[2]], {x, 0, 1}] == 0 /. intRule1 //. 
        intRule2;

Recall that the ith spatial derivative should be recovered back from int[i] as stated in the preceding section. Hence, we define

intRule5[mode_] = 
      int[mode a_ int[i_][b_, c__], d_] :> 
        int[mode HoldForm[D[a int[b, c], {x, i}]], d];

and rewrite SCond1 and SCond2 as

SCond = {SCond1 /. intRule5[(Subscript[ϕ, m])[x]], 
        SCond2 /. intRule5[(Subscript[ϕ, n])[x]]};

To simplify the notation in the solvability conditions, we use the orthonormality of the mode shapes and define the following parameters:

notationRule={int[(Subscript[ϕ, i])_[x]^2, {x, 0, 1}] -> 1,int[(Subscript[ϕ, i])_[x] (Subscript[ϕ, i_]^((4)))[x],{x, 0, 1}]->(Subscript[λ, i]^4),int[μ[x] (Subscript[ϕ, i])_[x]^2, {x, 0, 1}] -> (Subscript[μ, i]),int[F[x] (Subscript[ϕ, i])_[x], {x, 0, 1}] -> 
              (Subscript[ω, 1, i]) f};

To identify the forms of the nonlinear terms in the solvability conditions, we first identify all of the possible forms of the nonhomogeneous terms in the first-order problem. To this end, we let

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

Then, all of the possible forms of the nonlinear terms in the first-order problem are given by

\(cubicTerms = 
      \(Nest[Outer[Times, basicTerms, #]\ &, basicTerms, 2] // Flatten\) // 
        Union; 
cubicTerms /. displayRule\)

Out of these terms, only the terms that may lead to secular terms appear in the solvability conditions, which can be identified according to

secularTerms = 
      ((E^(-I omgList[[#]] (Subscript[T, 0]))) cubicTerms /. 
                    expRule[#] /. Exp[_ (Subscript[T, 0]) + _. ] -> 0 // Union // 
              Rest) & /@ {1, 2}

Next, we define the following parameters:

coef = 
      Table[Coefficient[SCond[[i, 1]], secularTerms[[i]]], {i, 2}] // 
        Collect[#, (Subscript[δ, 0])] &;

Clear[symbolList]

symbolList[1] = 
      MapIndexed[
        If[#1 =!= 0, 8  (Subscript[ω, 1, m]) (Subscript[α, 1, #2[[1]]]), 0]
             &, coef[[1]]]

symbolList[2] = 
      MapIndexed[
        If[#1 =!= 0, 8  (Subscript[ω, 2, n]) (Subscript[α, 2, #2[[1]]]), 0]
             &, coef[[2]]]

and express SCond in a more concise form as

eqMod = 
      Expand[-(1/omgList)[[#]] 
                ((SCond[[#, 1]] /. notationRule /. 
                        Thread[secularTerms[[#]] -> 0]) + 
                    symbolList[#] . secularTerms[[#]])] == 0 & /@
         {1, 2}

where

(symbolList[1] -> coef[[1]] // Thread // Union // Rest)
           /. int -> Integrate // Timing

(symbolList[2] -> coef[[2]] // Thread // Union // Rest)
           /. int -> Integrate // Timing

Modulation Equations in Polar Form

Using eqMod and the function PolarForm defined in the preceding section, we obtain the modulation equations in polar form and the definitions for (Subscript[γ, i]) as

PolarForm[eqMod]

Modulation Equations in Cartesian Form and Symmetry Property

Using eqMod and the function CartesianForm defined in the preceding section, we obtain the modulation equations in Cartesian form, the definitions for (Subscript[ν, i]), and the symmetry property as

CartesianForm[eqMod]

7.3.3 Discretization of the Continuous Problem

As an alternative, we apply the method of multiple scales to the discretized system of EOM and BC. To determine the discretized form, we express v[x, t] and w[x, t] in terms of the linear mode shapes (Subscript[ϕ, m])[x] and (Subscript[ϕ, n])[x], respectively, as

discretRule = {v -> Function[{x, t}, (Subscript[q, m])[t] (Subscript[ϕ, m])[x]], 
        w -> Function[{x, t}, (Subscript[q, n])[t] (Subscript[ϕ, n])[x]]};

modes = {(Subscript[ϕ, m])[x], (Subscript[ϕ, n])[x]};

where the (Subscript[q, i])[t] are the generalized coordinates. We consider only (Subscript[q, m])[t] and (Subscript[q, n])[t] because in the presence of damping, all other modes that are not directly or indirectly excited decay with time. Substituting discretRule into each of EOM, multiplying the results with (Subscript[ϕ, m])[x] and (Subscript[ϕ, n])[x], respectively, and integrating the outcome from x = 0 to x = 1 using intRule3, intRule4, intRule1, intRule2, and intRule5, we obtain the following discretized form of the equations describing the response of the beam:

eq733a = 
      Table[int[modes[[k]] # /. discretRule, {x, 0, 1}] /. 
                        intRule3 //. intRule4 /. intRule1 //. intRule2
                   /. intRule5[modes[[k]]] //. 
              int[a_, b__] -> HoldForm[Integrate[a, b]] & /@ 
          EOM[[k]], {k, 2}]

where the damping is assumed to be modal. Using the orthonormality of the mode shapes (Subscript[ϕ, i])[x], we define

notationRule={int[(Subscript[ϕ, i])_[x]^2, {x, 0, 1}] -> 1,int[(Subscript[ϕ, i])_[x] (Subscript[ϕ, i_]^((4)))[x],{x, 0, 1}]->(Subscript[λ, i]^4),int[μ[x] (Subscript[ϕ, i])_[x]^2, {x, 0, 1}] -> (Subscript[μ, i]),int[F[x] (Subscript[ϕ, i])_[x], {x, 0, 1}] -> 
              (Subscript[ω, 1, i]) f}/.int[a_, b__] -> HoldForm[Integrate[a, b]]

To simplify the expression for eq733a, we first find possible nonlinear terms as follows:

basicTerms = Table[D[(Subscript[q, m])[t], {t, i}], {i, 0, 2}]

linearTerms = Join[basicTerms, basicTerms /. m -> n]

cubicTerms = 
      ϵ Nest[Outer[Times, linearTerms, #] &, linearTerms, 2] // 
          Flatten // Union;

The coefficients corresponding to cubicTerms in eq733a are obtained as

coefList = 
      Coefficient[eq733a[[#, 2]], cubicTerms] & /@ {1, 2};

We replace these coefficients by (Subscript[α, i,  j]) as

symbolList = 
      MapIndexed[If[#1 =!= 0, (Subscript[α, Sequence @@ #2]), 0] &, 
        coefList, {2}];

Using notationRule, symbolList, and cubicTerms, we can rewrite eq733a as

EOM1 = 
      Table[eq733a[[k, 1]] == 
            (eq733a[[k, 2]] /. Thread[cubicTerms -> 0]) + 
              symbolList[[k]] . cubicTerms /. notationRule, {k, 2}]

Using the method of multiple scales, we seek a first-order uniform expansion in the form

solRule = 
      (Subscript[q, n])_ -> ((∑_(j = 0)^1) (ϵ^j) (Subscript[q, n, j])[#1, #2] &)
          ;

Transforming the total time derivatives in EOM1 into partial derivatives in terms of (Subscript[T, 0]) and (Subscript[T, 1]), substituting the solRule and betayRule into EOM1, expanding the result for small ϵ, and discarding terms of order higher than ϵ, we obtain

eq733b = 
      (EOM1 /. {(Subscript[q, s])_[t] -> (Subscript[q, s])[(Subscript[T, 0]), (Subscript[T, 1])], 
                    Derivative[n_][(Subscript[q, s])_][t] -> 
                      dt[n][(Subscript[q, s])[(Subscript[T, 0]), (Subscript[T, 1])]], t -> (Subscript[T, 0])} /. solRule
                 /. betayRule // ExpandAll) /. 
        (ϵ^(n_ /; n > 1)) -> 0;

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

eqEps = 
      Thread[CoefficientList[Subtract @@ #, ϵ] == 0] & /@ 
          eq733b // Transpose

Then it follows from eqEps[[1]] that

\(\(sol0 = {
        q\_\(m, 0\) -> 
          Function[{T\_0, T\_1}, 
            A\_1[T\_1]\ Exp[I\ ω\_\(1, m\)\ T\_0] + 
              A\&_\_1[T\_1]\ Exp[\(-I\)\ ω\_\(1, m\)\ T\_0]], 
        q\_\(n, 0\) -> 
          Function[{T\_0, T\_1}, 
            A\_2[T\_1]\ Exp[I\ ω\_\(2, n\)\ T\_0] + 
              A\&_\_2[T\_1]\ Exp[\(-I\)\ ω\_\(2, n\)\ T\_0]]}; \)\)

Substituting sol0 into eqEps[[2]] yields

order1Eq = 
      (#[[1]] & /@ eqEps[[1]] /. 
                    (Subscript[q, k_, 0]) -> (Subscript[q, k, 1])) == 
                (#[[1]] & /@ eqEps[[1]] /. 
                      (Subscript[q, k_, 0]) -> (Subscript[q, k, 1])) - 
                  (#[[1]] & /@ eqEps[[2]]) /. sol0 // 
            TrigToExp // ExpandAll // Thread;

We consider the case of primary resonance of the flexural mode in the y direction, \(TraditionalForm
      \`Ω ≈ ω\_\(1, \ m\)\). To describe the nearness of the primary resonance, we introduce the detuning parameter σ and define

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

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

OmgRule = 
      Solve[ResonanceConds, Drop[omgList, {#}]~Join~{Ω}][
            [1]] & /@ {1, 2}

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

We substitute ResonanceConds into the right-hand sides of order1Eq and obtain the solvability condition as

\(SCond = 
      Table[Coefficient[\(-order1Eq[\([i, 2]\)]\) /. expRule[i], 
            Exp[I\ omgList[\([i]\)]\ T\_0]] == 0, {i, 2}]; 

    SCond /. displayRule\)

where

# -> 
          (# /. (Subscript[α, i_, j_]) :> coefList[[i, j]] /. 
                    betayRule /. ϵ -> 0 // Expand // 
              Collect[#, (Subscript[δ, 0])] &) & /@ 
      (symbolList // Flatten // Union // Rest)

which is in agreement with that obtained by direct approach.

7.3.4 Method of Time-Averaged Lagrangian

As a second alternative, we derive the modulation equations by using the method of time-averaged Lagrangian. The nondimensional Lagrangian of the beam can be expressed as (Crespo da Silva and Glynn, 1978; Arafat, Nayfeh, and Chin, 1998)

Lagr1 = 
      Hold[(∫_0^1)
                  ((1/2) 
                        (ϵ 
                              ((((∂_t)
                                      (∫_0^x) (1/2) 
                                        (((((∂_x) v[x, t]))^2) + 
                                        ((((∂_x) w[x, t]))^2)) 
                                        \[DifferentialD]x))^2) + 
                            ((((∂_t) v[x, t]))^2) + 
                            ((((∂_t) w[x, t]))^2)) - 
                      ϵ (1 - (Subscript[β, y])) 
                        ((∂_(x, x))v[x, t] 
                            (∂_(x, x))w[x, t] 
                            (∫_0^x)(∂_(x, x))v[x, t] 
                                (∂_x) w[x, t] \[DifferentialD]x)
                           - (1/2) 
                        (((((∂_(x, x))w[x, t]))^2) + 
                            ϵ ((((∂_x) w[x, t]))^2) 
                              ((((∂_(x, x))w[x, t]))^2)) - 
                      (1/2) (Subscript[β, y]) 
                        (((((∂_(x, x))v[x, t]))^2) + 
                            ϵ ((((∂_x) v[x, t]))^2) 
                              ((((∂_(x, x))v[x, t]))^2) + 
                            ϵ 2 (∂_x) v[x, t] 
                              (∂_(x, x))v[x, t] 
                              (∂_x) w[x, t] 
                              (∂_(x, x))w[x, t]) - 
                      ϵ 
                        ((((1 - (Subscript[β, y])))^2)/(2 (Subscript[β, γ]))) 
                        (((((∫_1^x)(∂_(x, x))v[x, t] 
                                    (∂_(x, x))w[x, t] 
                                    \[DifferentialD]x))^2) + 
                            2 (∂_(x, x))v[x, t] 
                              (∂_(x, x))w[x, t] 
                              (∫_0^x)
                                  (∫_1^x)
                                    (∂_(x, x))v[x, t] 
                                    (∂_(x, x))w[x, t] 
                                    \[DifferentialD]x \[DifferentialD]x)
                           - ϵ F[x] Cos[Ω t] 
                        v[x, t]) \[DifferentialD]x] /. 
            Integrate -> int /. 
          HoldPattern[(∂_t) int[a_, b_]] :> 
            int[(∂_t) a, b] // ReleaseHold;

where Hold and int are used to keep Integrate from being evaluated. Transforming the total time derivatives into partial derivatives in terms of (Subscript[T, 0]) and (Subscript[T, 1]), we modify Lagr1 as

Lagr2 = 
      Lagr1 /. {v[x_, t] -> v[x, (Subscript[T, 0]), (Subscript[T, 1])], w[x_, t] -> w[x, (Subscript[T, 0]), (Subscript[T, 1])], 
          Derivative[m_, n_][w_][x_, t] -> 
            dt[n][D[w[x, (Subscript[T, 0]), (Subscript[T, 1])], {x, m}]], t -> (Subscript[T, 0])};

To this end, we let

\(\(solRule = {
        v -> Function[{x, T\_0, T\_1}, 
            A\_1[T\_1]\ ϕ\_m[x]\ Exp[I\ ω\_\(1, m\)\ T\_0] + 
              A\&_\_1[T\_1]\ ϕ\_m[x]\ 
                Exp[\(-I\)\ ω\_\(1, m\)\ T\_0]], 
        w -> Function[{x, T\_0, T\_1}, 
            A\_2[T\_1]\ ϕ\_n[x]\ Exp[I\ ω\_\(2, n\)\ T\_0] + 
              A\&_\_2[T\_1]\ ϕ\_n[x]\ 
                Exp[\(-I\)\ ω\_\(2, n\)\ T\_0]]}; \)\)

and define

intRule6 = {
        int[c1_.  int[a1__]^n_. , b1_] :> 
          c1 int[int[a1]^n /. intRule1 //. intRule2 /. 
                  ϵ -> 0, b1] /; FreeQ[c1, First[b1]], 
        int[c1_ int[a1__]^n_. , b1_] :> 
          int[c1 int[a1]^n /. intRule1 //. intRule2 /. 
              ϵ -> 0, b1]};

Substituting solRule and betayRule into Lagr2, using intRule1, intRule2, and intRule6 to simplify the expansions of the integrands, and collecting the coefficients of ϵ, we obtain

(order1Lagr = 
          Lagr2 /. solRule /. betayRule /. intRule1
                                 //. intRule2 //. intRule6 /. intRule1
                           //. intRule2 //. intRule6 /. intRule1 //. 
                  intRule2 // TrigToExp // Expand // 
            Coefficient[#, ϵ] &; ) // Timing

ResonanceConds = {(Subscript[ω, 2, n]) == (Subscript[ω, 1, m]), 
          Ω == (Subscript[ω, 1, m]) + ϵ σ} /. 
        Equal -> Rule

notationRule1={int[(Subscript[ϕ, i])_[x]^2, {x, 0, 1}] -> 1,int[F[x] (Subscript[ϕ, i])_[x], {x, 0, 1}] -> f (Subscript[ω, 1, i]),int[(((Subscript[ϕ, i_]^″)[x])^2),{x, 0, 1}]->(Subscript[λ, i]^4)};

Using ResonanceConds and notationRule1, we obtain the slowly varying terms from order1Lagr as

TAL = 
      order1Lagr /. 
            Exp[a_] :> 
              Exp[Expand[a /. ResonanceConds] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])]
                 /. Exp[_ (Subscript[T, 0]) + _. ] -> 0 /. notationRule1;

The Euler-Lagrange equations corresponding to the TAL can be written as

eqMod1=D[D[TAL,(Subscript[\(A\&_\), 1]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_1[T\_1]]\)==0//.int[a_, b__] -> HoldForm[Integrate[a, b]]

and

eqMod2=D[D[TAL,(Subscript[\(A\&_\), 2]^′)[(Subscript[T, 1])]],(Subscript[T, 1])]-\(D[TAL, A\&_\_2[T\_1]]\)==0//.int[a_, b__] -> HoldForm[Integrate[a, b]]

Adding linear viscous damping to eqMod1 and eqMod2 and performing integration by parts yields the same modulation equations as obtained by directly attacking the partial-differential system.