Applying the perturbation method of Multiple-Scales. 

	The authors present two examples of the application of the Multiple Scales solution procedure to a single-degree-of-freedom systems  (the pendulum equation, and Duffing equation) using Mathematica .

     by Raya Khanin and Matthew P.Cartmell

Very few nonlinear problems can be solved directly so that an exact solution is obtained. Usually, one has to make approximations and apply perturbation analysis to find the behavior of the nonlinear system. Perturbation methods work by applying small nonlinear perturbations to known linearised solutions. Their correct application is resticted to weakly nonlinear systems, so the nonlinear terms are small compared to the linear terms. This is usually the case, when the motions are finite but not very large. The correctness of perturbation analysis decreases for growing amplitudes of motion.
	 One might question the use of  perturbation methods to find solutions of differential equations, since powerful computers one can now simulate (almost) any system numerically.  The reason is that even the most powerful and fast numerical integration will never provide a researcher or engineer with an insight into what is going on in the system.  The intrinsic problem with all numerical methods is that the complex interactions of the system and the underlying physics of the problem are usually hidden.  Computers provide highly specific answers to highly specific questions. One obtains the response of a system only for particular sets of parameters, and for particular intervals of time. For nonlinear systems, one cannot just interpolate or exrapolate from discrete data. So — in particular for nonlinear systems — a theoretical insight into the physical nature of the probems is required for turning the huge amount of computer-generated information into useful understanding. 
	There are several perturbation methods for solving weakly non-linear equations  (for example, the straightforward expansion, the method of harmonic balance, several methods of averaging,  the method of strained coordinates,  the method of matched and composite solutions, and the method of multiple scales).  For a general overview on perturbation methods see, for example, Nayfeh [Nayfeh 1973].  To some extent, the choice of particular perturbation method is a matter of personal taste and habit. None of the perturbation methods is superior in general to any other. In difficult cases it is possible to modify a perturbation method, or to take a combination of different methods. If one method fails, another method might work, but generally speaking this approach is tentative, and considerable skill and intuition is necessary to get the desired result. 
	In this paper the essentials of perturbation methods (and Mathematica code to implement them) are presented and illustrated on a simple example of the pendulum equation, by first considering the method of straghtforward expansion. After demonstrating that this elementary method does not work for such a simple system, we turn to the more advanced perturbation method of Multiple Scales.  The implementation of the Multiple Scales method using Mathematica is given in detail for the pendulum and Duffing equations.

The pendulum equation

The pendulum equation is a standard model for demonstration of perturbation techniques. It has already been used for demonstrating how Mathematica  can be used to find exact power series solution (to arbitrary order) of differential equations by R.Pickett,  G.Lindgren and R.Anderson  in their notebook   "Time dependent solutions of the pendulum equations" submitted to MathSource.  
	We first consider the equation of an undamped and unforced pedulum  sin(x(t)) (ω^2)+((x(t)))''=0 subjected to initial conditions x (0) = (Subscript[x, 0]), x' (0) = (Subscript[((x^′)), 0]). Here x is the inclination,and  ω is the natural circular frrequency for small-amplitude swings. An attempt to find an exact solution to this equation by using DSolve[] fails. Assuming the pendulum rotations to be finite but not very large  sin(x) = x - ((x^3)/6) , yields

eqn = x''[t] +  ω^2x[t] - 1/6 ϵ x[t]^3 ==0;

A small parameter ϵ is introduced here simply to indicate that the nonlinear term is assumed to be small compared to the linear terms (see Thomsen, p.84). (In this case, DSolve[] will return the solution to this nonlinear equation. This lengthy solution contains special functions and is difficult to interpret physically. So, it is better to try to find an approximate and more informative solution.)

Straighforward expansion

The most elementary version of the perturbation method is to represent the solution in the power series of ϵ:  \(TraditionalForm
      \`x(t) = ∑\+\(i = 0\)\%Order ϵ\^i\ \(\((x(i))\) (t)\)\). This expansion should be uniformly valid for all times t>0. (x(i))(t) are unknown functions to be determined. To form perturbation equations, this power series is substituted into the original equation. Since the functions (x(i))(t)  are independent of ϵ, the coefficients of powers of ϵ must balance. The  straghtforward perturbation method has been implemented by Stephan Kaufmann. In addition to this method, his Mathematica  package covers the methods of strained coordinates and matched and composite solutions (MathSource: ¬http://www.mathsource.com/cgi-bin/MathSource/Applications/Mathematics).
	As will be seen from below this simple approach does not yield the correct results. It is included here to illustrate the fundamental principle of applying perturbations, and to show why other, more advanced, perturbation methods are needed.

expandedeqn=Expand[(eqn /.  x -> (x[0][#1]  +  ϵ x[1][#1]  +  
                          ϵ^2 x[2][#1]&))][[1]];
perteqns = 
        Thread[Equal[
            Table[Coefficient[expandedeqn,  ϵ,  i], 
               {i,  0,  2}],  0]]

The solution of the zero-th order of ϵ perturbation equation is

zerosol  =  
      DSolve[{perteqns[[1]],  x[0][0] == x0,  x[0]'[0] == x0'}, 
        x[0][t],  t]

Subtituting this solution into the first order of ϵ pertubation equation and solving it, yields

firstsol = 
      Simplify[x[1][t]  /.  
          Flatten[
              DSolve[Expand[{perteqns[[2]],  x[1][0] == 0,  
                      x[1]'[0] == 0}  /.  Flatten[zerosol]] ,  x[1], 
                 t]][[1]]]

The solution of the perturbation equations to first order ϵ contains terms of the type tCos[t] and tSin[t]. This implies that this solution is non-periodic whireas one would expect a free pendulum to oscillate periodically. Secondly, the terms of this type grow unbounded with time violating the assumption that the expansion  is uniformly valid. Proceeding to calculate the terms to higher order ϵ will not cure this problem since then the terms t^2, t^3, and so forth, will  enter the solution. Solution terms of the type t Cos[t+ϕ0] and t Sin[t+ϕ0] are called secular.
	 Thus, the perturbation  method of straightforward expansion does not work here. Other pertubation methods should be used to cure these problems and to produce proper solutions. For example, the method of harmonic balance or one of the methods of averaging. In this paper, however,  we concentrate on the method of Multiple Scales.

Method of Multiple Scales: theory

As its name suggests, the method of Multiple Scales permits the solution to a be a function of multiple independent variables, or scales. This powerful perturbation method is commonly applied in Structural Dynamics (Nayfeh, 1973). It has been also applied to a wide variety of problems in physics and applied mathematics where weakly non-linear equations are encountered. 
	The main idea of the Multiple Scales (MS) method is that the single  independent variable t  is  split up into several new independent variables (Subscript[T, 1]),  (Subscript[T, 2]),  ..., (Subscript[T, M]), where (Subscript[T, i]) = (ϵ^i) t .   Although such variables are not precisely  independent  (each being related to the others due to a simple scaling relationship via  the perturbation parameter, ϵ) they are treated as such. The dependent variable , x , is assumed in the form of the uniformly valid expansion \(TraditionalForm
      \`x(t; \ ϵ)\  = \ 
        ∑\+\(i = 0\)\%M\( ϵ\^i\) 
              \(\(x\_i\)(T\_0, \(T\_1...\) T\_M)\)\  + \ O(ϵT\_M)\). The error in this expansion is stated as O(ϵ(Subscript[T, M])) as a reminder that this expansion is valid for times up to\[RawEscape]O((ϵ^(-M))) . Beyond these times, we must use other time scales to keep the expansion uniformly valid.
	By using the chain rule, the time derivative is transformed according to  ((\[RawEscape]d)/dt) = 
              (∂/∂)(Subscript[T, 0]) +  
                  ϵ (∂/∂)(Subscript[T, 1])  +  
                  (ϵ^2) (∂/∂)(Subscript[T, 2]) +  ...  and \(TraditionalForm
      \`d\^2\/dt\^2 = 
        ∂\^2\/∂T\_0\^2 + \ 
            2  ϵ \(∂\/∂T\_0\) 
              ∂\/∂T\_1 + 
            \(ϵ\^2\)(
              ∂\^2\/∂T\_1\^2 + 
                \(∂\/∂T\_0\) ∂\/∂T\_2)
              \  + \ ...\) .  Substituting these formulas into the original differential equations and equating coefficients of like powers of ϵ provides equations for determining  the functions (Subscript[x, 0]),  (Subscript[x, 1]),  ... (Subscript[x, M]). The solutions of these functions contain arbitrary functions of the time scales (Subscript[T, 1]),  (Subscript[T, 2]),  ..., (Subscript[T, M]). In order to determine these functions, additional conditions need to be imposed. If the expansion of x(t;  ϵ) is to be uniformly valid, i.e. valid  for times as large as (ϵ^(-M)) , then (ϵ^i) (Subscript[x, i]) should be a small correction to (ϵ^((i - 1)))(Subscript[x, i - 1]) for any i.  Thus, we require that\(x\_i\/x\_\(i - 1\) < ∞  \)for  all  (Subscript[T, 1]),  (Subscript[T, 2]),  ..., (Subscript[T, M]) .   Note that this condition does not mean that each  (Subscript[x, i]) is bounded. It is equivalent to the elimination of secular terms. This provides the necessary equations for the arbitrary functions to be determined.

Applying Multiple Scale method to the Pendulum equation

The analysis presented in this notebook closely follows analysis presented in the book ``Vibrations and Stability'' by J.J. Thomsen (1997). Consider the pendulum equation as above. In pure function form this equation can be written as

InputEq = 
  D[#1, #2, #2] + ω^2#1 - 1/6*ϵ*#1^3&;

where #1 and #2 are formal parameters representing x[t] and t respectively. Using upvalues for symbols we declare all system parameters to be real, and positive. For the work with complex functions we load the package

\(<< Algebra`ReIm`; 
		
Parameters\  = \ {ω, \ ϵ}; 

    Map[UpSet[Sign[#], \ 1]&, \ Parameters]; 

    Map[UpSet[Im[#], \ 0]&, \ Parameters]; 
Im[T[0]] ^= 0; \ Im[T[1]] ^= 0; 
    \)

We seek a first-order approximate solution:

ExpansionOrder = 1;

Function InPowers[InputEq_, u_List, t_, T_] substitutes the expansion of dependent variable x in the form \(TraditionalForm
      \`x[t]\  = \ 
        \(x[0]\)[T[0], \ T[1]]\  + \ ϵ\ \(x[1]\)[T[0], T[1]]\  + \ 
          O(ϵ\^2)\), (T[0] = t,  T[1] = ϵ T[0]) into the original equation and by equating like powers of ϵ returns a list of perturbation equations.

\(Options[InPowers]\  = 
      \ {ExpOrder -> Global`ExpansionOrder, \ 

        SmallPara -> Global`ϵ}; \ 

    InPowers[InputEq_, \ u_, \ t_, \ T_, \ opts___] := 

      Module[{utemp, \ i, \ ExpandedEqn, \ expansionorder, \ smallpara, \ 
          arguments, \ bigarguments, \ trules}, 
        
{expansionorder, \ smallpara}\  = 
          \ {ExpOrder, \ SmallPara} /. Flatten[{opts, \ Options[InPowers]}]; 
        
arguments\  = \ Table[t*smallpara^i, \ {i, \ 0, \ expansionorder}]; 
        
utemp[t]\  = \ 
          Sum[smallpara^i*\(u[i]\)[Sequence@@arguments], 
            \ {i, \ 0, \ expansionorder}]; 

        bigarguments\  = Table[T[i], \ {i, \ 0, \ expansionorder}]; 

        trules\  = \ MapThread[Rule, \ {arguments, \ bigarguments}]; 		
	
        ExpandedEqn\  = \ InputEq[utemp[t], \ t] /. trules; \ 

        Table[Coefficient[ExpandedEqn, \ smallpara, \ j], 
          \ {j, \ 0, \ 
expansionorder}]
]; \)

Applying InPowers[] to our input equation InputEq results in the following perturbation equations:  to order (ϵ^0)

\(\({zeroordereqn, \ firstordereqn} = InPowers[InputEq, \ x, \ t, \ T]; 

    zeroordereqn\  == \ 0\ \)\)

and to order (ϵ^1)

firstordereqn == 0

The next step is finding the general solution to the zero-th order equation. Note, that in order for the MS method to be  succesfully applied, the zero-th order equation should be a linear ode with constant coefficients, so that its solution can be readily generated.  In this case, the perturbation equation to zero-th order ϵ has the form of the undamped and unforced linear pendulum. If only  a single time-scale T(0) was present, the solution would have the familiar form x(0)=aCos(ϕ+T(0)). With two independent variables T(0) and T(1), the same solution applies with the difference that now both amplitude and phase are functions of T(1). Thus, the solution of this equation is x(0) = a (T(1) cos(T(0)+ϕ(T(1)))). With the method of Multiple Scales it is convenient to employ complex exponentials, as below

zerosol = DSolve[zeroordereqn == 0, x[0], {T[0], T[1]}]/.{C[_][#2]*Cos[omega_.*#1]+C[_][#2]*Sin[omega_.*#1] -> A[#2]*Exp[I*omega*#1]+Conjugate[A[#2]]*Exp[-I*omega*#1]}

It should be noted that old versions of DSolve[] will not solve this type of equation considering it to be a partial differential equation.  Here the formal parameter #1 stands  for  fast  time  T[0]  and  the  formal parameter #2   stands  for slow time  T[1].  For further analysis let us define the natural frequency of the zero-th order perturbation equation as

natfreq  =  ω;

The zeroth order ϵ solution is  then substituted into the first-order ϵ equation

firsteqn = Collect[Expand[TrigToExp[firstordereqn/.zerosol[[1]]]], E^n_];
firsteqn == 0

Terms containing the factor  \(TraditionalForm
      \`e\^\(i\ ω\ \(T(0)\)\)\) are resonant because they correspond to harmonic excitation at the natural  frequency of the linear system. These resonant excitation terms will produce secular terms of the form \(TraditionalForm
      \`\(T(0)\)\ e\^\(i\ ω\ \(T(0)\)\)\) in the solution for the perturbation equation to first order ϵ. However, we are free to choose the arbitrary functionA(T(1))  from the requirement that the solution for x(1) should be free of secular terms. By requiring all the coefficients at  (e^(ΙωT(0)))to vanish the solvability condition is obtained. (Note that all complex functions (like (E^Ι)ωT[0]) that appear in the equation have their conjugates that appear in the equation in exactly the same way so that we only need to consider one of them).

\(solveq1\  = \ 
      DeleteCases[Coefficient[firsteqn, \ E^\((I*T[0]*natfreq)\)], \ 
        a_*E^b_]; 
solveq1\  == 0\)

With the above condition fulfilled, the secular terms are removed from the equation, which now constitutes a linear oscillator with a single harmonic forcing term \(TraditionalForm
      \`1\/6\ \(A(T(1))\)\^3\ e\^\(3\ ω\ \(T(1)\)\)\). When the excitation is harmonic, one is primarily interested in the steady-state behavior. The standard procedure is to ignore the homogeneous part of the equation. We thus find the particular solution by first applying DSolve[] and then deleting terms containing arbitrary constants

resteqn1 = DeleteCases[firsteqn, a__*E^(I*T[0]*ω) | a__*E^(-I*T[0]*ω)];

statres = x[1][T[0], T[1]]/.DSolve[resteqn1 ==0 , x[1], {T[0], T[1]}][[1]];

statres = TrigToExp[DeleteCases[statres, a__*C[_][_]]]/. 
      E^n_ :> E^(Collect[Expand[n], T[0]])

For further work with derivatives of complex functions we should add at least one rule:

\(\(Unprotect[Im, \ Derivative]; 

    Im[\(\(Derivative[n_]\)[f_]\)[x_]]\  := \ 
      \(\(Derivative[n]\)[Im[f[x]]]\)[x]\  /; \ Im[x] === 0; 

    Protect[Derivative, Im]; 		\)\)

The arbitrary function \[RawEscape]A(T(1))   is to be determined from the solvability condition. The  usual trick is to express it in polar form \(TraditionalForm
      \`A(T(1)) = 
        1\/2\ \(a(T(1))\)\ 
          e\^\(Ιϕ(T(1))\)\), where a(T(1)) ,ϕ(T(1))  are real valued functions  for amplitude and phase of the  real argument\[RawEscape]T(1).   Substituting this  expression for A into the solvability conditions and separating the real and imaginary parts we obtain modulation equations for amplitude a and phase ϕ:

\(\(solveqtmp\  = \ 
      Expand[solveq1\  /. \ A -> \ \((1/2*a[#] Exp[I*ϕ[#]]&)\)]; 

    Im[T[1]] ^= 0; 
RealValued[ϕ]; RealValued[a]; 


    solveqtmp\  = \ 
      Expand[solveqtmp*E^\((\(-I\)*ϕ[T[1]])\)\  /. \ 
          E^n_\  :> \ E^\((Expand[n])\)]; 
s1\  = \ Re[solveqtmp]; 

    s2\  = \ Im[solveqtmp]; 
ColumnForm[{s1\  == 0, \ s2\  == 0}]\ \)\)

The above equations are the modulation equations of this perturbation analysis. They govern the modulations in slow time T(1) of the amplitudes and phases of those terms which are oscillating harmonically in time T(0). To determine the stationary response one might attempt to solve these differential equations to find out what is happening with the amplitudes and phases when time approaches infinity:

\(\(modsol = 
      DSolve[{s1 == 0, \ s2 == 0, \ a[0] == a0, \ ϕ[0] == ϕ0}\  /. 
          \ T[1] -> t, \ {a, \ ϕ}, \ t]
\)\)

Limit[{a[t],  ϕ[t]}  /.  modsol[[1]],  t -> ∞]

The approximate near-resonant solution can now be written as (x(t) = 
                O((ϵ^2)) + (x(0)) (T(0), T(1)) + 
                  ϵ (x(1)) (T(0), T(1)))if we recall all substitutions made so far. Function FinalAnswer[] makes all the required substitutions:

\(\(Clear[FinalAnswer]; 

    FinalAnswer[expr_, \ subst_] := 

      Collect[\(ExpToTrig[expr\  /. \ A -> \((1/2*a[#]*Exp[I*ϕ[#]]&)\)]\ 
             /. \ Flatten[{subst, \ T[0] -> t, \ T[1] -> ϵ*t}]\)\ 
           /. \ Cos[n_]\  :> \ Cos[Collect[Expand[n], \ t]], \ Cos[_]
		]; 
    
			\)\)

x[t;  ϵ]  =  
      (FinalAnswer[x[0][T[0],  T[1]] /.  zerosol[[1]],  modsol] 
             +  ϵ*FinalAnswer[statres,  modsol])  /.  modsol

As expected, the solution is periodic. The base frequency of the oscillations depends upon the oscillation amplitude a0, a typical non-linear phenomenon.

Duffing equation

Let us now consider a more complex example and elaborate more on finding resonance conditions and stationary responses. Consider the Duffing equation with parametric excitation (see Thomsen, 1997, p94):

\(eqn\  = \ 
      \(\(x'\)'\)[t]\  + \ 2\ ϵ\ β\ \(x'\)[t]\  + \ 
          \((1 - \ ϵ\ q\ Ω^2\ Cos[ω\ t])\) 
            x[t]\  + \ ϵ\ γ\ x[t]^3 == 0; 

    InputEq\  = \ 
\ \ 
      D[#1, \ #2, \ #2]\  + \ 2*β*ϵ*D[#1, \ #2] + \ 
          ω^2  #1\  - \ 
          ϵ*Ω^2*q*Cos[Ω*#2]*#1\  + 
          ϵ*γ*#1^3&; 

    Parameters\  = 
      \ {ϵ, \ β, \ q, \ γ, \ ω, \ 
        Ω, \ σ}; 

    Map[UpSet[Sign[#], \ 1]&, \ Parameters]; 

    Map[UpSet[Im[#], \ 0]&, \ Parameters]; \)

Equations of this type arise in numerous problem of nonlinear mechanics. With γ = -1/6 this equation becomes the pendulum equation with forcing and damping included. 
The initial steps of the perturbation analysis are straightforward:

{zeroordereqn,  firstordereqn} = InPowers[InputEq,  x,  t,  T]

zerosol = DSolve[zeroordereqn == 0, x[0], {T[0], T[1]}]/.{C[_][#2]*Cos[omega_.*#1]+C[_][#2]*Sin[omega_.*#1] -> A[#2]*Exp[I*omega*#1]+Conjugate[A[#2]]*Exp[-I*omega*#1]};

firsteqn = Expand[TrigToExp[firstordereqn/.zerosol[[1]]]];

As we have seen already, terms oscillating at a frequency equal to the natural frequency of the homogeneous equation are resonant, and their sum should be equated to zero to get rid of secular terms. Terms oscillating at a frequency near the natural frequency cause small divisor terms to appear in the solution for x[1], and small divisor terms are just as unacceptable as secular terms. Consequently, small divisor terms should be eliminated by equating their sum to zero. 
	The strength of the MS method is that once perturbation equations are set up it is easy to tell what special cases should be considered. The procedure for determining these special cases is performed in the function ResonanceConditions[eqn_, NatFreq_]. This function examines exponential frequencies of the first order perturbation equation, deletes those that cannot result in resonance conditions and returns a list of expressions which, when they are near zero, result in primary resonance conditions.

Unprotect[Exponent]; 
Exponent[expr_, patt_]:=Union[Join[Cases[expr, patt^_, -1], Cases[expr, patt]]/.Power[_, n_]:>n/.patt:>1];
Protect[Exponent];

Clear[ResonanceConditions];
Options[ResonanceConditions] = {Targ->Global`T[0]};
ResonanceConditions[expr_, NatFreq_, opts___] := 
Module[{targ, freqq},
{targ} = {Targ}/.Flatten[{opts, Options[ResonanceConditions]}];
freqq = Factor[Union[DeleteCases[Exponent[expr, E], 1]]]/(I*targ);
freqq = DeleteCases[freqq, (a_?Positive + NatFreq /; FreeQ[a, NatFreq]) | (a_.*NatFreq /; FreeQ[a, NatFreq]) | (a_ + b_*NatFreq /; !Positive[a/(1-b)]) | (a_ /; !Positive[a])];
freqq = ExprForResCond[#, NatFreq]&/@freqq;
freqq = Select[freqq, !AtomQ[#]&];
Union[freqq]
];

Clear[ExprForResCond];
ExprForResCond[a_?Positive+ b_.*NatFreq_, NatFreq_] := Return[a-(1-b)*NatFreq] /; FreeQ[a, NatFreq];
ExprForResCond[a_?Negative+b_.*NatFreq_, NatFreq_]:=
Return[-a+(1-b)*NatFreq] /; FreeQ[a, NatFreq];
ExprForResCond[a_, NatFreq_] := 
Return[a -NatFreq] /; FreeQ[a, NatFreq];

freq = ResonanceConditions[firsteqn, natfreq]

Thus, for the Duffing equation, there are two special cases to consider: the non-resonant  case ( Ω is away from 2  ω), and the primary resonant case (Ω ∼ 2  ω).
	Once the resonance conditions for the problem are determined the analysis goes in several independent directions. For each resonance condition the rules for the detuning parameter should be formulated and then substituted into the original equations. The rest of the analysis for this (and other systems) closely follows analysis presented for the pendulum equation with the main difference that a new phase variable (or variables in the case of the multi-degree-of-freedom systems) will have to be introduced to obtain an autonomous set of modulation equations. 
	To measure the closeness of Ω to 2ω a detuning parameter σ is introduced: \(TraditionalForm
      \`\[RawEscape] Ω\  = \ 
        2  ω\  + \ ϵσ\). We should substitute an expression for the detuning parameter into the oscillatory terms. It is unneccesary (albeit not harmful) to substitute the detuning parameter into the non-oscillatory terms. In other words, all occurrences of Ω-2ω or Ω in the exponents will be changed to ϵσ or 2ω+ϵσ respectively. The function σrule[] introduces rules to perform these substitutions.

σrule[a_+b_.*NatFreq_, NatFreq_] := Return[{E^(x_.*(a + b*NatFreq) + y_.)->E^(x*ϵ*σ + y), E^(z_.*(x_.*a + 
y_.))->E^(z*(x*(-b)*NatFreq + x*ϵ*σ + y))}] /; 
(FreeQ[a, NatFreq]);

firsteqntmp = firsteqn /. σrule[freq[[1]], natfreq];

This is just a part of the Mathematica code for this function. A fuller version of this function has been developed in the package MultipleScale.m to deal with different types of resonance conditions. 
	A solvability condition is obtained by requiring all the cofficients of (E^I)ωT[0] in the first  ordeϵ perturbation equation to vanish:

\(\(solveq2\  = \ 
      DeleteCases[Coefficient[firsteqntmp, \ E^\((I*T[0]*natfreq)\)], \ 
        a__*E^b__\  /; \ \(! MatchQ[b, \ ϵ*T[0]*c___]\)]; 

    solveq2\  == 0\ \)\)

With the removal of secular terms the stationary response is given by

\(\(resteqn2\  = \ 
      DeleteCases[firsteqntmp, \ 
        a__*E^\((I*T[0]*ω)\)\  | \ 
          a__*E^\((I*T[0]*ω\  + \ ϵ*T[0]*c_.)\)\  | \ \ 
          a__*E^\((\(-I\)*T[0]*ω)\)\  | \ \ 
          a__*E^\((\(-I\)*T[0]*ω\  + \ ϵ*T[0]*c_.)\)]; 

    statres\  = \ 
      \(x[1]\)[T[0], \ T[1]] /. 
        \(DSolve[resteqn2\  == 0\ , \ x[1], \ {T[0], \ T[1]}]\)[\([1]\)]; 

    statres = 
      \(TrigToExp[DeleteCases[statres, \ a__*\(C[_]\)[_]]] /. \ 
          σ*ϵ -> freq[\([1]\)]\)\  /. \ 
        E^n_\  :> \ E^\((Collect[Expand[n], \ T[0]])\)\ \)\)

To construct the modulation equations the polar form of the arbitrary function A[T[1]]  is substituted into the solvability condition:

solveqtmp = solveq2 /. A-> (1/2*a[#]*Exp[I*ϕ[#]]&);
Im[T[1]]^=0; RealValued[ϕ];
solveqtmp = Expand[solveqtmp*E^(-I*ϕ[T[1]]) /. ϵ*T[0]->T[1] /. E^n_ :> E^(Expand[n])];
solveqtmp == 0

A change of variables is needed to obtain autonomous modulation equations, i.e. to eliminate time in the equations. Let us define  a new real valued function ψ by finding exponential frequencies in the solvability condition, and then substituting it back into the condition to eliminate the time variable

\(subst[t_]\  = \ 
      \((Expand[Exponent[solveqtmp, \ E]/I]\  /. \ ϵ*T[0] -> T[1])
            \)\  - \ ψ[T[1]]\  /. \ T[1] -> t; 
Im[T[1]] ^= 0; \ 
    RealValued[ψ]; 


    ϕsubst = 
      Flatten[{Solve[subst[t]\  == 0, \ ϕ[t]], \ 
            Solve[\(subst'\)[t]\  == 0\ , \ \(ϕ'\)[t]]}] /. \ t -> T[1]

    \)\(solveqtmp = 
      Expand[\(solveqtmp /. \ ϕsubst\)\  /. \ 
          E^n_\  :> \ E^\((Expand[n])\)]; 
    
{s1, \ s2} = {Re[solveqtmp], \ Im[solveqtmp]}; 
    
{s1\  == \ 0, \ s2 == 0}\)

The case when the modulation equations are readily solved is very rare in practice. In general, modulation equations have quite complicated form. They  do not, however, have to be solved to obtain stationary responses. Instead one  seeks  the singular points of the modulation equations by setting the derivatives of the amplitude a and the phase ψ to zero:

\({as1, \ as2}\  = 
      \ {s1, \ s2} /. 
        \ {\(a'\)[_] -> 0, \ \(ψ'\)[_] -> 0, \ a[_] -> a, \ 
          ψ[_] -> ψ}; 
ColumnForm[{as1\  == 0, \ as2\  == 0}]\)

To find all possible solutions for a and ψ one can directly apply Solve[] to the modulation equations: modsol = Solve[{as1 == 0,  as2  == 0},  {a, ψ}]. There will be 9 pairs of solutions in this case including  a zero solution a = 0 resulting in  x(t) = 0.
	To find solutions for nonzero amplitude a one can usea further approach which has proved to be quite efficient for a great number of different problems (see Nayfeh and Mook, 1979).  Note, however, that there is no general recipe for studying the singular points of the modulation equations. First we solve the modulation equations for Cos[ψ] and Sin[ψ], and then by taking (Cos[ψ]^2  +  Sin[ψ]^2) we  find a frequency response curve and expression for Tan[ψ]

cossinψ = Solve[{as1 ==0, as2 ==0}, {Cos[ψ], Sin[ψ]}];

freqresponsecurve = Simplify[Simplify[Cos[ψ] /. cossinψ[[1]]]^2 + Simplify[Sin[ψ] /. cossinψ[[1]]]^2]==1;
tangentvalue=(Simplify[Sin[ψ]/. cossinψ[[1]]])/(
    Simplify[Cos[ψ]/. cossinψ[[1]]]);
{freqresponsecurve, Tan[ψ]==tangentvalue}

This approach quite often results in a biquadratic equation for a. Thus, solving it for (a^2) we obtain

singpointa = Solve[freqresponsecurve /. a^2 -> b, b]/.b->a^2

Solutions for a exist only if 2 σ ω - 
                (Sqrt[-16 (β^2) (ω^2) + 
                    (q^2) (Ω^4)]) >= 0 && -16 (β^2) (ω^2) + (q^2) (Ω^4) 
               >= 0 \[RawEscape]
	The last step in the Multiple Scale perturbation analysis is to check the stability of the stationary solutions. The determination of solution stability reduces down to the determination of the stability of the singular points for an autonomous system of modulation equations. Let us examine the stability of the two nonzero solutions defined by singpointa and cossinψ. For this, we first define the derivatives \(TraditionalForm\`\(a\^'\)\) and \(TraditionalForm\`\(ψ\^'\)\) as functions of a and ψ from the modulation equations, calculate the Jacobian of these equations and find their eigenvalues  at the singular point determined by cossin :

Clear[ader, ψder];
MapThread[SetDelayed[#1, #2]&, {{ader[a_, ψ_], ψder[a_, ψ_]}, Simplify[{a'[T[1]], ψ'[T[1]]} /. Solve[{s1 ==0, s2 == 0}, {a'[T[1]], ψ'[T[1]]}][[1]] /. {a[_]->a, ψ[_]->ψ}]}];
Clear[auxilary, MSJacobian];
auxilary[y_, x_List]:= Map[D[y, #]&, x];
MSJacobian[y_List, x_List] := Map[auxilary[#, x]&, y];
lambdas = Eigenvalues[MSJacobian[{ader[a, ψ], ψder[a, ψ]}, {a, ψ}] /. cossinψ[[1]]]

The solution is unstable if at least one eigenvalue contains a positive real part, i.e. if

ineq = Thread[Re[lambdas] >  0]

These inequalities are true if -9 (a^4) (γ^2) + 
              12 (a^2) γ σ ω + 
              16 (β^2) (ω^2)>0\[RawEscape]and 
\(TraditionalForm
      \`\(\[RawEscape]\(-4\)\ β\ ω + 
          \@\(\(-9\)\ a\^4\ γ\^2 + 
              12\ a\^2\ γ\ σ\ ω + 
              16\ β\^2\ ω\^2\)\)\/\(4\ ω\)\)> 0  which is fulfilled if  \(TraditionalForm
      \`\(-9\)\ a\^4\ γ\^2 + 12\ a\^2\ γ\ σ\ ω > 0
        \). If \(TraditionalForm\`\[RawEscape] γ < 0  \)this  condition becomes

cond=-3 (a^2) γ + 4  σ ω < 0;

Substituting values for \(TraditionalForm\`\(a\^2\) \)into  this condition we obtain

\(\(\(ReplaceAll[#, \ expr___\  < 0\  :> \ Expand[expr]\  < \ 0]&\)/@cond
        \  /. \ singpointa
\)\)

Recalling conditions for the existing solutions to a, we obtain that for negative γ  the first solution is unstable, while the second is stable. For positive  γ it is vice versa. Although it is easy to perform the stability analysis in symbolic form for a set of two equations in Mathematica, it becomes more complicated as the number of equations grows. Thus, it is desirable to have a symbolic-numerical procedure which will examine solutions and determine under what conditions solutions exist and are stable.  Adding a new feature to function FinalAnswer[] to account for the resonance case, we obtain the solution in terms of the original variables

\(Clear[FinalAnswer]; 

    FinalAnswer[expr_, \ subst_, \ freq_] := 
      Collect[ExpToTrig[
          \(\(expr\  /. \ \ A -> \((1/2*a[#]*Exp[I*ϕ[#]]&)\)\)\  //. \ 
              Flatten[{subst, \ 
                  Solve[freq - ϵ*σ == 0, \ σ], \ 
                  T[0] -> t, \ T[1] -> ϵ*t, \ a[t_] -> a, \ 
                  ψ[t_] -> ψ}]\)\  /. \ 
            E^n_\  :> \ E^\((Expand[n])\)], \ Cos[n_]]; \)

\(\(	x[t; \ ϵ]\  = \ 
      FinalAnswer[\(x[0]\)[T[0], \ T[1]] /. \ zerosol[\([1]\)], \ 
          ϕsubst, \ freq]\  + 
        ϵ*FinalAnswer[statres, \ ϕsubst, \ freq]\)\)

Conclusions

This article summarises the intial part of a research project that is aimed at developing a parallelised computational package using Mathematica to analyse large-scale problems encountered in Structural Dynamics and also in other areas. In this paper we illustrated the implementation of the Multiple Scales solution procedure for a single-degree-of freedom system typically represented by an undamped pendulum or by Duffing system. We have shown how the Multiple Scales method works and how Mathematica can be used to perform this perturbation analysis.  
	For multiple-degree-of-freedom systems the Multiple Scales method is implemented in the same step-by-step manner. For this, the functions and pieces of code presented in this article are generalised and extended to deal with the many different cases that can be encountered.  
	Based on many existing examples we have developed a  Mathematica package MultipleScale.m  for the application of the Multiple Scales method to sets of weakly non-linear second-order ordinary differential equations  (Khanin and Cartmell, 1998). The package has been tested on numerous problems in Vibration Theory.  Initial experiments with parallelisation of the Multiple Scales method have been also started on a multi-processor computer.

Acknowledgements

The authors gratefully acknowledge the help and guidance provided by Dr. Anthony Gilbert of the Department of Mathematics and Statistics, University of Edinburgh. We would also like to thank the EPSRC for financial support, Dr. Tom Wickam-Jones of Wolfram Research Inc. for his comments and Mr. Alistair MacLean of Consafe Engineering UK Ltd for his interest in this work and for the provision of case study materials. We are also grateful to the Editor of Mathematica in Education and Research for his suggestions for improving this paper.

References

Gilbert A.D., A geometric approach to two-timing.
Math. Proc. Camb. Phil Soc. 1980. Vol.87, 167-178.

Khanin  R. and M.P.Cartmell., 
Computational Approach to applying the Multiple Scales Method to Weakly Nonlinear Problems. Internal Report D/1998/MPC/RK2. Department of Mechanical Engineering, University of Glasgow, Glasgow, G12 8QQ, Scotland, UK, 1998.

Nayfeh A.H., Perturbation methods.
John Wiley&Sons. New York, 1973

Nayfeh A.H. and Mook D.T., Nonlinear oscillations.
John Wiley&Sons. New York, 1979 

Thomsen J.J.,  Vibrations and stability. Order and chaos.
McGraw-Hill, 1997

About the authors

Raya Khanin received her Ph.D. in Applied Mathematics from the Weizmann Institute of Science, Israel in 1996. She now works as a Research Fellow in the Department of Mechanical Engineering, University of Glasgow, Scotland, UK.  

Matthew Cartmell graduated with a Ph.D. in Nonlinear Structural Dynamics from the University of Edinburgh in 1984 and is now Professor of  Applied Dynamics in the Department of Mechanical Engineering, University of Glasgow, Scotland, UK.

Department of Mechanical Engineering
James Watt Building
University of Glasgow
Glasgow G12 8QQ,
Scotland, UK

emails: raya@mech.gla.ac.uk and m.cartmell@mech.gla.ac.uk

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file MultScales.nb containing Mathematica code for the material described in this column.