Chapter 2
The Duffing Equation

Off[General::spell1]

To symbolize some notations, we need to load the package

Needs["Utilities`Notation`"]

2.1 The Duffing Equation

The free oscillations of many conservative systems having a single degree of freedom are governed by an equation of the form

\(\(\(eq21a = \(x\^″\)[t] + f[x[t]] == 0; 
    \)
 (*\ An\ equation\ is\ set\ and\ then\ assigned\ to\ eq21a\ *) \)\)

where f is a nonlinear function of x[t]. Here, (x^″)[t] is the acceleration of the system, whereas f[x[t]] is the restoring force. Let x = (Subscript[x, 0]) be an equilibrium position of the system. Then (Subscript[x, 0]) is a constant and hence f[(Subscript[x, 0])] = 0. In order for Mathematica to interpret (Subscript[x, 0]) as a constant, we need to symbolize it. Otherwise,  if we replace x with an expression g[t], Mathematica will replace (Subscript[x, 0]) with g[t] sub 0 and not treat it as a constant. For example, replacing x with Sin[t] in the expression (x^2) + (Subscript[x, 0]) yields

(x^2) + (Subscript[x, 0]) /. x -> Sin[t]

Clearly, (Subscript[x, 0]) is treated as a variable x with subscript 0, which is wrong. Therefore, we need to symbolize  (Subscript[x, 0]) as  (enter and select (Subscript[x, 0]) followed by clicking Symbolize[◼] from NotationPalette)

Symbolize[(Subscript[x, 0])];

Then, Mathematica treats the symbolized expression (Subscript[x, 0]) as a constant. Consequently, replacing x with  Sin[t] in the expression (x^2) + (Subscript[x, 0]) yields

(x^2) +  (Subscript[x, 0]) /. x -> Sin[t]

which is correct.

Next, we assume that f is an analytic function of x at x = (Subscript[x, 0]) and expand it in a Taylor series around (Subscript[x, 0]) as

fexp = Series[f[x[t]], {x[t], (Subscript[x, 0]), 3}]

where terms up to cubic are retained. Using Normal to truncate the higher-order terms from fexp, we have

fpoly = fexp // Normal

Substituting fpoly for f[x[t]] and using the equilibrium condition f[(Subscript[x, 0])] = 0, we obtain

temp = eq21a /. {f[x[t]] -> fpoly} /. f[(Subscript[x, 0])] -> 0

Now, we introduce the transformation x[t] → (Subscript[x, 0]) + u[t] in temp to obtain the final form of the equation that will be used in the subsequent analysis. To accomplish this, we need to express this transformation in a pure function form as x → ((Subscript[x, 0]) + u[#] &); otherwise, Mathematica will not transform the derivatives of x[t] in terms of the derivatives of u[t] unless we explicitly define the rules for those derivatives; that is, it will not replace (x^′)[t] with (u^′)[t] and (x^″)[t] with (u^″)[t] by simply substituting x[t] → (Subscript[x, 0]) + u[t]. Moreover, it will not replace x[p] with (Subscript[x, 0]) + u[p] as desired. To see the difference between these two transformations, we compare the following two statements:

{x[t], x[p], (x^′)[t], (x^″)[t]} /. 
      x[t] -> (Subscript[x, 0]) + u[t]

{x[t], x[p], x[a, b], y[t], (x^′)[t], 
        (x^″)[t]} /. x -> ((Subscript[x, 0]) + u[#]&)

We note that the pure function form can be used for any expression with the head x, especially, for the derivatives. Consequently, we introduce this transformation in pure function form into temp and obtain

eq21b = temp /. x -> ((Subscript[x, 0]) + u[#]&)

To simplify the notation in eq21b, we denote the coefficient of u[t]^b by (Subscript[k, b]). The result is

eq21c = eq21b /. _.  u[t]^b_.  -> (Subscript[k, b]) u[t]^b

We relate the (Subscript[k, b]) to the derivatives of the function f[x] by using the built-in function Cases as

Cases[eq21b, a_.  u[t]^b_.  -> {(Subscript[k, b]) -> a}, {2}] // Flatten

Most of this chapter is devoted to the special case of cubic nonlinearities. Thus, we set the coefficient of the quadratic term in eq21c equal to zero and obtain

eq21d = eq21c /. (Subscript[k, 2]) -> 0

which is called the Duffing equation. We note that (Subscript[k, 1]) and (Subscript[k, 3]) may be positive or negative.

It is a good practice to nondimensionalize the governing equations before treating them with perturbation methods. To this end, we nondimensionalize u and t using a characteristic length U and a characteristic time T of the motion and denote the nondimensional variables by an asterisk. In nondimensional form, eq21d becomes

\(eq21e = 
      \(eq21d /. u -> \((U\ \(u\^*\)[#/T]\ &)\)\) /. t -> T\ \(t\^*\)\)

Next, we multiply the left-hand side of eq21e (using eq21e[[1]]) by (T^2)/U so that the coefficient of \(TraditionalForm\`\(\(u\^*\)\^″\)[\(t\^*\)]\) is unity. The result is

eq21f = (eq21e[[1]] (T^2)/U // Expand) == 0

We choose T so that the linear natural frequency of the system is unity; that is, we let (Subscript[k, 1]) (T^2) = 1, for ease of notation, we let ϵ = (Subscript[k, 3]) (T^2) (U^2), drop the asterisk, and rewrite eq21f as

\(DuffingEq = 
      \(eq21f /. {k\_1\ T\^2 -> 1, k\_3\ T\^2\ U\^2 -> ϵ}\) /. 
        \(anything_\^*\) -> anything\)

We note that ϵ is a nondimensional quantity, which is a measure of the strength of the nonlinearity.

2.2 Straightforward Expansion

To solve the DuffingEq, we augment it with the initial conditions

ic = {u[0] == u0, (u^′)[0] == v0};

The solution u of the initial-value problem (IVP), consisting of the DuffingEq and associated initial conditions ic, is a function of the independent variable t and the parameter ϵ. Next, we determine an approximation of the IVP for weak nonlinearity; that is, for small but finite ϵ.

First, we set ϵ equal to zero in the DuffingEq, augment the result with the initial conditions by using the Mathematica command Join, and obtain the following linear IVP:

linearIVP = Join[{DuffingEq /. ϵ -> 0}, ic]

The solution of this linear IVP can be obtained by using the Mathematica command DSolve. The result is

linearSol = DSolve[linearIVP, u[t], t][[1]]

When ϵ is small but different from zero, the solution of the nonlinear initial-value problem is no longer given by linearSol, but deviates from it. We try a correction in the form of a power series in ϵ; that is, we expand u[t] in a power series in ϵ as

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

where (Subscript[u, 0])[t] is the solution of the linear problem. For example, to first order, we have

u[t] /. uExpRule[1]

and to third order, we have

u[t] /. uExpRule[3]

and to fifth order, we have

u[t] /. uExpRule[5]

Next, we restrict our discussion to first order; that is, we keep two terms in the power series. Thus, we apply the rule uExpRule[1] to the DuffingEq, expand the result, and obtain

eq22atemp = DuffingEq /. uExpRule[1] // ExpandAll

Since we kept only terms up to O(ϵ) in the expansion of u[t], we need to neglect terms of order higher than ϵ in eq22atemp for consistency. To accomplish this, we use the rule (ϵ^(n /; n > k)) → 0 to discard terms with order higher than (ϵ^k) and obtain

eq22a = eq22atemp /. (ϵ^(n_ /; n > 1)) -> 0

Setting the coefficients of like powers of ϵ in eq22a equal to zero, we obtain the following hierarchy of linear equations:

eqEps = CoefficientList[eq22a[[1]], ϵ] == 0 // Thread

which needs to be solved in succession. That is, we first solve the first equation in eqEps to obtain (Subscript[u, 0])[t] and then substitute the result into the second equation in eqEps to obtain a linear nonhomogeneous equation, which can be solved for (Subscript[u, 1])[t].

The general solution of the zeroth-order problem (eqEps[[1]]) can be expressed in Cartesian form as

DSolve[eqEps[[1]], (Subscript[u, 0])[t], t][[1]]

Alternatively, we can express the solution of the zeroth-order problem in the polar form \(TraditionalForm
      \`u\_0[t] → \(a\_0\) Cos[t + β\_0]\) by adding two appropriate initial conditions to eqEps[[1]] as

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

where (Subscript[a, 0]) and (Subscript[b, 0]) are arbitrary constants. Because eqEps[[1]] represents a linear undamped oscillator, one can write down its general solution in either Cartesian or polar form without using the Mathematica command DSolve.

Substituting for (Subscript[u, 0])[t] in the first-order equation (eqEps[[2]]) and determining the general solution of the resulting equation for (Subscript[u, 1])[t], we obtain

sol1temp = 
      DSolve[eqEps[[2]] /. sol0, (Subscript[u, 1])[t], t][[1]] // Simplify

We note that the homogeneous part of sol1temp consists of the terms involving the C[i] (the usual case) and Cos[t + (Subscript[β, 0])] (obtained from DSolve in this case). Hence, we can replace them by (Subscript[a, 1]) Cos[t + (Subscript[β, 1])] and rewrite sol1temp as

sol1 = 
      (Subscript[u, 1])[t] -> 
        (sol1temp[[1, 2]] /. {C[_] -> 0, Cos[t + _] -> 0}) + 
          (Subscript[a, 1]) Cos[t + (Subscript[β, 1])]

Combining the solutions of the zeroth- and first-order problems, we obtain the general solution of the DuffingEq, to first order, as

sol = u[t] == (u[t] /. uExpRule[1] /. sol0 /. sol1)

where (Subscript[a, 0]), (Subscript[a, 1]), (Subscript[β, 0]), and (Subscript[β, 1]) are arbitrary constants.

We started with a second-order equation that can support two initial conditions but it appears that we ended up with four arbitrary constants. It turns out that only two of the constants (Subscript[a, 0]), (Subscript[a, 1]), (Subscript[β, 0]), and (Subscript[β, 1]) are arbitrary and that the two initial conditions in ic are sufficient to determine all of them. To see this, we first change the Equal form in ic to the Rule form as

icRule = ic /. Equal -> Rule

Then, we impose the initial conditions icRule by evaluating sol and its derivative at t = 0 according to

eq22b = {sol, D[sol, t]} /. t -> 0 /. icRule

We equate the coefficient of each power of ϵ in eq22b and obtain four algebraic equations for the constants (Subscript[a, 0]), (Subscript[a, 1]), (Subscript[β, 0]), and (Subscript[β, 1]). The result is

eq22c = 
      Table[Map[Coefficient[#, ϵ, i] &, eq22b, {2}], {i, 0, 1}]

Thus, once (Subscript[a, 0]) and (Subscript[β, 0]) are known from eq22c[[1]], (Subscript[a, 1]) and (Subscript[β, 1]) can be calculated from eq22c[[2]].

Alternatively, instead of including the homogeneous solution at each order and consider the arbitrary constants to be independent of ϵ as done above, one may disregard the homogeneous solution at all orders except the first and consider the arbitrary constants to depend on ϵ in imposing the initial conditions. With the latter approach, the solution to first order can be rewritten as

sol /. {(Subscript[a, 0]) -> a, (Subscript[β, 0]) -> β, (Subscript[a, 1]) -> 0} // ExpandAll

We note that the correction term \(TraditionalForm
      \`\(1\/32\) \(a\^3\) ϵ\ Cos[3  t + 3  β] - 
        \(3\/8\) \(a\^3\) t\ ϵ\ Sin[t + β]\) is small compared to the principal term a Cos[t + β], as it is supposed to be, only when ϵ t is small compared with unity. When ϵ t is O(1), the term that is supposed to be a small correction becomes the same order as the principal term. Moreover, when ϵ t > O(1), the "small-correction" term becomes larger than the principal term. Hence, the above obtained straightforward expansion is valid only for times such that ϵ t < O(1); that is, t < O((ϵ^(-1))). Consequently, we say such expansions are nonuniform or breakdown for long times and we call them pedestrian  or naive expansions. The reason for the breakdown of the above expansion is the presence of the term t Sin[t + β], a product of algebraic and circular terms. Such terms are called mixed-secular terms. In subsequent sections, we implement four methods that avoid secular terms, and hence yield uniform expansions.

2.3 The Lindstedt-Poincaré Technique

The breakdown in the straightforward expansion is due to its failure to account for the nonlinear dependence of the frequency of the system on the nonlinearity. A number of techniques that yield uniformly valid expansions have been developed. Four of these techniques are discussed in this chapter. We start with the Lindstedt-Poincaré technique in this section.

To account for the dependence of the frequency ω of the system on the nonlinearity, we explicitly exhibit ω in the governing differential equation. To this end, we introduce the transformation τ = ω t, where ω is a constant that depends on ϵ, and obtain

DuffingEq=u''[t]+u[t]+ϵ u[t]^3==0;
eq23a = 
        DuffingEq /. u -> (u[ω #] &) /. 
          t -> τ/ω

To determine a uniform expansion of the solution of eq23a, we expand both u[τ; ϵ] and ω[ϵ] in powers of ϵ using uExpRule[k] and

omgRule[maxOrder_] := 
      ω -> 
        (Subscript[ω, 0]) + Sum[(ϵ^i) (Subscript[ω, i]), {i, maxOrder}]

where (Subscript[ω, 0]) is the linear natural frequency of the system. In the present case, the linear natural frequency was normalized to unity. Hence, we set (Subscript[ω, 0]) = 1. The corrections to the linear frequency are determined in the course of the analysis by requiring the expansion of u to be uniform for all τ. Next, we show how one can determine (Subscript[ω, 1]).

Substituting the expansions of u and ω into eq23a, using the fact that (Subscript[ω, 0]) = 1, setting maxOrder = 1, expanding the result, and discarding terms of order higher than ϵ, we obtain

eq23b = 
      ExpandAll[eq23a /. uExpRule[1] /. omgRule[1]] /. 
          (Subscript[ω, 0]) -> 1 /. (ϵ^(n_ /; n > 1)) -> 0

Equating coefficients of like powers of ϵ on both sides of eq23b yields

eqEps = CoefficientList[eq23b[[1]], ϵ] == 0 // Thread

The general solution of eqEps[[1]], the zeroth-order equation, can be expressed in polar form as

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

where a and β are constants. We have expressed (Subscript[u, 0]) as a pure function so that we can evaluate its derivatives in the higher-order problems.

Substituting for (Subscript[u, 0])[τ] in eqEps[[2]], the first-order equation, yields

order1Eq = eqEps[[2]] /. sol0 // Simplify

whose particular solution is

u1pSol = 
      DSolve[order1Eq, (Subscript[u, 1])[τ], τ][[1]] /. C[_] -> 0 // 
        Simplify

It is clear that u1pSol contains a mixed-secular term, which makes the expansion nonuniform. In contrast with the straightforward expansion, where the secular term cannot be annihilated unless u is trivial, in this case, we can choose the parameter (Subscript[ω, 1]) to eliminate the secular term according to

omg1Rule = 
      Solve[Coefficient[u1pSol[[1, 2]], Sin[β + τ]] == 0, 
          (Subscript[ω, 1])][[1]]

Substituting for (Subscript[ω, 1]) into u1pSol, we have

sol1 = u1pSol /. omg1Rule // Simplify

We note that, to determine the condition for the elimination of the secular term from (Subscript[u, 1]), we do not need to determine the particular solution first as done above. Instead, we only need to inspect order1Eq and choose (Subscript[ω, 1]) so that the coefficient of Cos[β + τ], which produces secular terms in (Subscript[u, 1]), is equal to zero. To this end, we calculate the nonhomogeneous part of order1Eq and simplify the result using trigonometric identities as

expr23a = order1Eq[[1]] /. (Subscript[u, 1]) -> (0&) // TrigReduce

Next, we choose (Subscript[ω, 1]) to annihilate the term in expr23a that produces secular terms, namely Cos[β + τ], and obtain

Solve[Coefficient[expr23a, Cos[β + τ]] == 0, (Subscript[ω, 1])]
        [[1]]

Substituting sol0 and sol1 into uExpRule[1] yields

u[τ] = 
      (u[τ] /. uExpRule[1] /. sol0 /. sol1 // Simplify)

Replacing τ with ω t, substituting for ω, using the fact that (Subscript[ω, 0]) = 1, and expanding the arguments of the trigonometric functions, we obtain

solLP = 
      u[t] == (
          u[τ] /. τ -> ω t /. omgRule[1] /. 
                omg1Rule /. (Subscript[ω, 0]) -> 1 /. 
            f_[arg_] :> f[arg // Expand])

Clearly, the expansion solLP is free of secular terms and the correction term \(TraditionalForm
      \`\(1\/32\) a\^3\ ϵ\ 
        Cos[3  t + 3  β + \(9\/8\) ϵ\ a\^2\ t]\) is small compared with the principal term \(TraditionalForm
      \`a\ Cos[t + β + \(3\/8\) ϵ\ a\^2\ t]\) for all t and hence it is uniformly valid.

2.4 The Method of Multiple Scales

We note from solLP that the functional dependence of u on t and ϵ is not disjoint. In fact, u depends on the combination ϵ t as well as on the individual t and ϵ. Carrying out the expansion to higher order, we find that u depends on the combinations ϵ t, (ϵ^2) t,  (ϵ^3) t,⋯ as well as on the individual t and ϵ. Hence, u = u(ϵ t, (ϵ^2) t,  (ϵ^3) t,⋯; ϵ). For small ϵ, (Subscript[T, n]) = (ϵ^n) t, for n = 1, 2, 3, ⋯, represent different time scales. For example, if ϵ = (1/60), variations on the time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]) can be observed, respectively, on the second, minute, and hour arms of a watch. Thus, (Subscript[T, 0]) represents a fast time scale, (Subscript[T, 1]) represents a slower time scale, (Subscript[T, 2]) represents an even slower time scale, and so on. Since the dependence of u on t and ϵ occurs on different time scales, we imagine that we have a watch and attempt to observe the behavior of u using the different scales of the watch. Therefore, instead of determining u as a function of t, we determine u as a function of (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2]), ⋯. To this end, we change the independent variable t in the governing equation from t to (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2]), ⋯ , which are symbolized by

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

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

We define the derivative operators as

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;

and treat the time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]) as independent variables.

Next, we use the method of multiple scales to determine an approximate solution for the DuffingEq. In Section 2.4.1, we treat the DuffingEq in its second-order form; in Section 2.4.2, we transform it first into a system of two real-valued equations; and in Section 2.4.3, we transform it first into a single complex-valued equation.

2.4.1 Second-Order Real-Valued System

DuffingEq=u''[t]+u[t]+ϵ u[t]^3==0;

To determine a first-order uniform expansion of the solution of the Duffing equation, we use the two time scales (Subscript[T, 0]) and (Subscript[T, 1]) and hence let

maxOrder = 1;

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

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

Consequently, we have replaced the original ordinary-differential equation by a partial-differential equation, and it appears that, as a result, we have complicated the problem. This is partially true, but experience with this method has shown that the disadvantages of introducing this complication are far outweighed by the advantages. Not only does this method provide a uniform expansion, but it also provides all the various nonlinear resonance phenomena, as we shall see in subsequent chapters.

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

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

Substituting this expansion into eq241a, expanding the result, and discarding terms of order higher than ϵ, we obtain

eq241b = 
      (eq241a /. solRule // ExpandAll) /. 
        (ϵ^(n_ /; n > maxOrder)) -> 0

Equating coefficients of like powers of ϵ in eq241b yields

eqEps = CoefficientList[eq241b[[1]], ϵ] == 0 // Thread

The general solution of eqEps[[1]] can be written in a pure function form as

sol0 = (Subscript[u, 0]) -> (a[#2] Cos[#1 + β[#2]] &);

where #1 stands for the scale (Subscript[T, 0]) and #2 stands for the scale (Subscript[T, 1]). The functional dependence of a and β on (Subscript[T, 1]) is not known at this level of approximation; it is determined at subsequent levels of approximation by eliminating the secular terms.

Substituting sol0 into eqEps[[2]] and moving the nonhomogeneous parts to the right-hand side of the resulting equation yields

eq241c = 
      (eqEps[[1, 1]] /. (Subscript[u, 0]) -> (Subscript[u, 1])) - # & /@ 
        Reverse[eqEps[[2]] /. sol0 // Simplify]

Next, we expand the right-hand side of eq241c in a Fourier series using trigonometric identities and obtain

eq241c[[2]] = eq241c[[2]] // TrigFactor // Expand

To produce a uniform expansion, we eliminate the terms that produce secular terms from eq241c[[2]]; that is, we choose ((Subscript[a, 1])^′)[(Subscript[T, 1])] and (β^′)[(Subscript[T, 1])] to annihilate each of the coefficients of Sin[(Subscript[T, 0]) + β[(Subscript[T, 1])]] and Cos[(Subscript[T, 0]) + β[(Subscript[T, 1])]]. The result is

eq241d = 
      Coefficient[eq241c[[2]], #] == 0 & /@
         {Sin[(Subscript[T, 0]) + β[(Subscript[T, 1])]], Cos[(Subscript[T, 0]) + β[(Subscript[T, 1])]]}

Solving eq241d for ((Subscript[a, 1])^′)[(Subscript[T, 1])] and (β^′)[(Subscript[T, 1])] when a[(Subscript[T, 1])] ≠ 0 yields

SCond=Solve[eq241d,{a'[(Subscript[T, 1])],β'[(Subscript[T, 1])]}][[1]]

The solution of SCond[[1]] is a = (Subscript[a, 0]) = constant. Then, the solution of SCond[[2]] is

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

With the above solvability conditions, eq241c becomes

eq241e = eq241c /. SCond

The particular solution of eq241e can be obtained as

sol1 = 
      DSolve[eq241e, (Subscript[u, 1])[(Subscript[T, 0]), (Subscript[T, 1])], {(Subscript[T, 0]), (Subscript[T, 1])}][[1]] /. 
          C[_] -> (0&) // Simplify

eq241f = 
      u[(Subscript[T, 0]), (Subscript[T, 1])] == 
        (u[(Subscript[T, 0]), (Subscript[T, 1])] /. solRule /. sol0 /. sol1 /. 
              a[(Subscript[T, 1])] -> (Subscript[a, 0]) /. betaRule)

where (Subscript[a, 0]) and (Subscript[β, 0]) are constants to within the order of the error indicated. In terms of the original variable t, eq241f can be expressed as

u[t] == eq241f[[2]] /. {(Subscript[T, 0]) -> t, (Subscript[T, 1]) -> ϵ t}

in agreement with the expansion obtained by using the Lindstedt-Poincaré technique.

In the higher-order approximations, we include the scales (Subscript[T, 0]), (Subscript[T, 1]),  ⋯, (Subscript[T, N]) but we do not include the term O[(ϵ^N)] in an N th-order expansion.

Before closing this section, we present an alternate representation of the solutions of the perturbation equations. Instead of the real-valued form sol0, we represent the solution of eqEps[[1]] in a complex-valued form; that is,

\(\(sol0 = 
      u\_0 -> \((A[#2]\ Exp[I\ #1] + A\&_[#2]\ Exp[\(-I\)\ #1]\ &)\); \)\)

where \(TraditionalForm\`A\&_\) is the complex conjugate of A and

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

Substituting sol0 into eqEps[[2]] and moving the nonhomogeneous terms to the right-hand side of the resulting equation, we have

eq241g=(Subscript[u, 1])[(Subscript[T, 0]), (Subscript[T, 1])]+(Subscript[u, 1]^((2, 0)))[(Subscript[T, 0]), (Subscript[T, 1])]-# & /@ (eqEps[[2]] /. sol0)//ExpandAll//Reverse

We note that the terms proportional to Exp[I (Subscript[T, 0])] and Exp[-I (Subscript[T, 0])] produce secular terms in the particular solution of (Subscript[u, 1]). Thus, to produce a uniform expansion, we set the coefficient of each of these functions equal to zero; that is,

eq241h = 
      Coefficient[eq241g[[2]], #] == 0 & /@
         {Exp[I (Subscript[T, 0])], Exp[-I (Subscript[T, 0])]}

These two equations are not independent because they are complex conjugates of each other. Hence, if one of them is satisfied, the other is automatically satisfied.

To analyze the solutions of eq241h[[1]], we multiply it by Exp[-I β[(Subscript[T, 1])]] and replace A with its polar form ruleA. The result is

eq241i = 
      Expand[eq241h[[1, 1]] Exp[-I β[(Subscript[T, 1])]] /. ruleA] == 0

Next, we separate the real and imaginary parts of eq241i. To accomplish this, we define the following rule:

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

Using the realRule, we find that the imaginary part of eq241i is

ampEq = Im[eq241i[[1]]] == 0 /. realRule

which governs the modulation of the amplitude. Similarly, the real part of eq241i is

phaseEq = Re[eq241i[[1]]] == 0 /. realRule

which governs the modulation of the phase. These modulation equations are in agreement with those obtained above by expressing the solution in real form.

Comparing the complex-valued and the real-valued representations, we find it more convenient to use the complex-valued form. Therefore, the complex-valued form is used in the remainder of this book.

2.4.2 First-Order Real-Valued System

We let v[t] = (u^′)[t] and transform the DuffingEq into the following set of two first-order equations:

eq242a={(u^′)[t] == v[t],v'[t]+u[t]==-ϵ u[t]^3};

We seek a first-order uniform expansion of the solution of eq242a in terms of the two time scales (Subscript[T, 0]) and (Subscript[T, 1]) in the form

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

Substituting solRule into eq242a, transforming the derivative with respect to t in terms of the derivatives with respect to (Subscript[T, 0]) and (Subscript[T, 1]), and discarding terms of order higher than ϵ, we obtain

eq242b=(eq242a/.{(u_)'[t]->dt[1][u[(Subscript[T, 0]), (Subscript[T, 1])]],u_[t] -> u[(Subscript[T, 0]), (Subscript[T, 1])]}/.solRule//ExpandAll)/.(ϵ^(n_ /; n > 1)) -> 0

Equating coefficients of like powers of ϵ, we have

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

Zeroth-Order Problem: Linear System

The zeroth-order problem is given by the linear system

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

To determine the solution of this linear system, we seek a solution proportional to Exp[I (Subscript[T, 0])] as

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

Next, we determine the coefficient matrix as

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

To determine the eigenvalues and eigenvectors of coefMat, we define the following conjugate rule and Hermitian matrix:

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

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

Then the left and right eigenvectors of coefMat are given by

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

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

To express the solution of the zeroth-order problem, we introduce the following basic function:

basicH = A[(Subscript[T, 1])] (E^(I (Subscript[T, 0])));

In terms of this function and the right eigenvector, the solution of the zeroth-order problem can be expressed in terms of

sol0Form = rightVec basicH

and its complex conjugate as

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

We have expressed the solution of the zeroth-order problem in function form so that its partial derivatives can be readily evaluated.

First-Order Problem

Substituting the zeroth-order solution sol0 into the first-order problem eqEps[[2]] and moving the nonhomogeneous terms to the right-hand sides of the resulting equations, we obtain

\(order1Eq = 
      \((linearSys /. u_\_0 -> u\_1)\) == 
          \((linearSys /. u_\_0 -> u\_1)\) - 
            \((\(Subtract\ @@\ #\ &\)\ /@\ eqEps[\([2]\)] /. sol0 // Expand)
              \) // Thread\)

Eliminating the terms that produce secular terms in order1Eq demands that their right-hand sides are orthogonal to every solution of the adjoint homogeneous problem. To determine these conditions, we first calculate the coefficient vector of Exp[I (Subscript[T, 0])] as

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

Then, we demand that this vector is orthogonal to the left eigenvector and obtain the solvability condition

SCond = Expand[Conjugate[leftVec] . STerms] == 0

which is in agreement with that obtained by treating the second-order form of the DuffingEq. Repeating the procedure for the coefficient vector of Exp[-I (Subscript[T, 0])], we obtain the complex conjugate of SCond.

2.4.3 First-Order Complex-Valued System

In this section, we determine a first-order uniform expansion of the solution of the DuffingEq by transforming it first into a single first-order complex-valued equation. To this end, we introduce the transformation

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

whose inverse is

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

It follows from the DuffingEq that the acceleration is related to the displacement according to

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

Differentiating ζ[t], zetaRule[[1]], with respect to t and substituting for u[t], (u^′)[t],  and (u^″)[t] using the acceleration and transfRule, we obtain

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

Next, we use the method of multiple scales to determine a first-order uniform expansion of the solution of eq243a in the form

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

Substituting this expansion into eq243a, expanding the result, and discarding terms of order higher than ϵ, we have

eq243b=(eq243a/.{ζ[t] -> ζ[(Subscript[T, 0]), (Subscript[T, 1])],\(ζ\&_[t] -> ζ\&_[T\_0, T\_1]\),ζ'[t]->dt[1][ζ[(Subscript[T, 0]), (Subscript[T, 1])]]}/.solRule//ExpandAll)/.(ϵ^(n_ /; n > 1)) -> 0

Equating coefficients of like power of ϵ in eq243b yields

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

The solution of the zeroth-order problem, eqEps[[1]], can be expressed as

sol0Form = 
      DSolve[eqEps[[1]], (Subscript[ζ, 0])[(Subscript[T, 0]), (Subscript[T, 1])], {(Subscript[T, 0]), (Subscript[T, 1])}][
          [1]] /. C[1] -> A

which we express in function form as

sol0 = {
        (Subscript[ζ, 0]) -> 
          Function[{(Subscript[T, 0]), (Subscript[T, 1])}, sol0Form[[1, 2]] // Evaluate]}

To evaluate the complex conjugates of A and ζ, we define the conjugate rule

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

Substituting sol0 into the first-order problem, eqEps[[2]], using conjugateRule, and moving the nonhomogeneous terms to the right-hand side of the resulting equation, we obtain

order1Eq = 
      (eqEps[[1, 1]] /. (Subscript[ζ, 0]) -> (Subscript[ζ, 1])) == 
        (eqEps[[1, 1]] /. (Subscript[ζ, 0]) -> (Subscript[ζ, 1])) - 
          (eqEps[[2, 1]] /. sol0 /. (sol0 /. conjugateRule))

Elimintaing the terms that lead to secular terms from order1Eq yields the solvability condition

SCond = Coefficient[order1Eq[[2]], (E^(I (Subscript[T, 0])))] == 0

in agreement with that obtained by treating the second-order form of the DuffingEq as well as that obtained by transforming it first into two first-order real-valued equations.

2.5 Variation of Parameters

In the next section, we use the method of averaging to determine a first-order uniform expansion of the solution of the Duffing equation. To this end, we use the method of variation of parameters to transform it into a system of two first-order equations. The Duffing equation is repeated here; that is,

DuffingEq=u''[t]+u[t]+ϵ u[t]^3==0;

When ϵ = 0, the solution of this equation can be written as

usolEq = u[t] == a Cos[t + β];

where a and β are constants, which are sometimes referred to as parameters. It follows from usolEq that

cond01 = D[usolEq, t]

When ϵ ≠ 0, we assume that the solution of DuffingEq is still given by usolEq but with time-varying a and β; that is,

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

and hence

usolEqt = usolEq /. tdepRule

In other words, we consider usolEq as a transformation from u[t] to a[t] and β[t]. This is why this approach is called the method of variation of parameters. Using this view, we note that we have two equations, namely DuffingEq and usolEq, for the three unknowns u[t], a[t], and β[t]. Hence, we have the freedom of imposing a third condition (third equation). This condition is arbitrary except that it must be independent of DuffingEq and usolEq. This arbitrariness can be used to advantage, namely to produce a simple and convenient transformation. Out of all possible conditions, we choose to impose the condition cond01, thereby assuming that u[t] as well as (u^′)[t] have the same form as in the linear case. This condition leads to a convenient transformation; it leads to a set of first-order rather than second-order equations for a[t] and β[t].

Differentiating usolEq with respect to t and recalling that a and β are functions of t, we have

cond02 = D[usolEqt, t]

Comparing cond02 with cond01, we conclude that

cond1 = 
      Expand[cond02[[2]] - (cond01[[2]] /. tdepRule)] == 0

Differentiating cond01 with respect to t, we obtain

cond03 = D[cond01 /. tdepRule, t]

Substituting for u[t] and (u^″)[t] from usolEqt and cond03 into DuffingEq, we have

cond2 = 
      DuffingEq /. ({usolEqt, cond03} /. Equal -> Rule) // ExpandAll

Solving cond1 and cond2 for (a^′)[t] and (β^′)[t], we obtain the desired two first-order equations

transformedEq=(Solve[{cond1, cond2},{a'[t],β'[t]}][[1]]//Simplify)/.Rule -> Equal

if a ≢ 0. Thus, the original second-order DuffingEq for u[t] has been replaced by two first-order equations for a[t] and β[t]. We emphasize that no approximations have been made in deriving transformedEq. Comparing transformedEq with DuffingEq, we find that the transformed equations are more nonlinear than the original equation. Then, the question arises what is the value of this transformation? The answer depends on the value of ϵ. If ϵ is small, the major parts of a and β vary more slowly with t than u. This fact can be used to advantage analytically and numerically. The analytical advantage is utilized in the method of averaging, as discussed in the next section. Numerically, it is advantageous to solve the transformed equations instead of the original equation because a large step size can be used in the integration. This is the reason why astronomers numerically solve the variational equations rather than the original equations. Usually, astronomers and celestial mechanicians refer to this approach as the special method of perturbations.

2.6 The Method of Averaging

In this section, we determine a first-order approximation to the transformed equations obtained in the preceding section. To this end, we rewrite transformedEq as

ampEq = TrigReduce /@ transformedEq[[1]]

phaseEq = TrigReduce /@ transformedEq[[2]]

The major parts of a and β are slowly varying functions of time if ϵ is small. Hence, they change very little during the time interval π (the period of the circular functions) and, to the first approximation, they can be considered constant in the interval [0, π]. Hence, we replace a[t] and β[t] on the right-hand sides of ampEq and phaseEq with time-independent a and β and obtain

ampEq[[2]] = ampEq[[2]] /. (Reverse /@ tdepRule)

phaseEq[[2]] = phaseEq[[2]] /. (Reverse /@ tdepRule)

We average the right-hand sides of ampEq and phaseEq over the interval [0, π] and obtain

list26a = 
      (1/π) Integrate[#, {t, 0, π}]& /@
           {ampEq[[2]], phaseEq[[2]]} // Expand

It follows from the above result that

averagingEq = {(a^′)[t], (β^′)[t]} == 
          list26a // Thread

This averaging method is usually referred to as the Krylov-Bogoliubov or van der Pol technique.

Solving the averaged equations yields

rule26a = 
      DSolve[averagingEq, {a[t], β[t]}, t][[1]] /. {
          C[1] -> (Subscript[a, 0]), C[2] -> (Subscript[β, 0])}

Substituting for a[t] and β[t] into usolEq, we obtain, to the first approximation, that

usolEq /. tdepRule /. rule26a /. 
      Cos[arg_] :> Cos[Expand[arg]]

in agreement with the solutions obtained by using the Lindstedt-Poincaré technique and the method of multiple scales.

Before closing this section, we note that one can arrive at the final results in ampEq and phaseEq without going through the averaging process. The right-hand sides of ampEq and phaseEq are the sum of two groups of terms — a group that is a linear combination of fast varying terms and a group that is a linear combination of slowly varying terms. Then, to the first approximation, (a^′)[t] in ampEq is equal to the slowly varying group on its right-hand side, which is zero. And, to the first approximation, (β^′)[t] in phaseEq is equal to the slowly varying group on its right-hand side, which is (3/8) ϵ (a^2).