Methods For Solving Singular Perturbation Problems Arising In Science.pdf

Mathematical and Computer Modelling 54 (2011) 556–575

Contents lists available at ScienceDirect

Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm

Methods for solving singular perturbation problems arising in science and engineering Manoj Kumar, Parul ∗ Department of Mathematics, Motilal Nehru National Institute of Technology, Allahabad-211004 (U.P.), India

article

info

Article history: Received 26 August 2010 Received in revised form 25 February 2011 Accepted 25 February 2011 Keywords: Singular perturbations Asymptotic expansion method Matched asymptotic expansion method Multiple-scale method WKB method Poincare–Lindstedt method

abstract Singular perturbation problems are of common occurrence in all branches of applied mathematics and engineering. These problems are encountered in various fields such as solid mechanics, fluid dynamics, quantum mechanics, optimal control, chemical reactor theory, aerodynamics, reaction–diffusion processes, geophysics etc. In this paper, the basic methods and literature for solving the singular perturbation problems have been presented with their comparative study. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Perturbation theory is a vast collection of mathematical methods used to obtain approximate solution to problems that have no closed-form analytical solution. The methods work by reducing a hard problem to an infinite sequence of relatively easy problems that can be solved analytically. Perturbation problems depend on a small positive parameter. This parameter affects the problem in such a way that the solution varies rapidly in some region of the problem domain and slowly in other parts. The region where the solution varies rapidly is called the inner region. It is a well-known fact that the singularly perturbed boundary value problem possesses boundary or interior layer i.e. regions of rapid change in the solution near the end points or some interior points with width O(1) as ε → 0. Hence the numerical treatment of singularly perturbed differential equations gives rise to computational difficulties. The problem just considered is a regular perturbation problem, that is, its solution varies smoothly as the perturbation parameter ε approaches zero. The problems of interest here are singular perturbation problems, their solutions change abruptly in some way as ε reaches zero. Some or all solutions either might cease to exist or might become infinite or degenerate. A singular perturbation generally occurs when a problem’s small parameter multiplies its highest operator. Thus naively taking the parameter to be zero changes the very nature of the problem. In the case of differential equations, boundary conditions cannot be satisfied whereas in algebraic equations, the possible numbers of solutions are decreased. In recent years, a large number of special purpose methods have been developed to provide accurate numerical solutions. For details, one may refer to the books of Farrell et al. [1], Roos et al. [2], Miller et al. [3], Smith [4] and Morton [5] and the references therein. Zheng [6] proposed the use of control points of the Bernstein Bezier form to solve differential equations numerically, Kumar et al. [7–9], Reddy et al. [10] discussed initial and boundary value techniques for singularly perturbed two-point boundary value problems and many others. Singular perturbation theory is a vast and rich ongoing area of exploration for mathematicians, physicists, and other researchers. There are various methods which are used to tackle problems in this field. The more basic of these include the asymptotic

∗

Corresponding author. E-mail addresses: [email protected] (M. Kumar), [email protected] (Parul).

0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.02.045

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

557

expansion, method of matched asymptotic expansions and WKB approximation for spatial problems and in time, the method of multiple scales, and the Poincare–Lindstedt method. For a detailed description of singular perturbation in ordinary and partial differential equations we refer to [11–18]. A very readable introduction can also be found in [19]. The main purpose of this paper is to present a comprehensive overview of perturbation theory and explain mathematical methods for obtaining approximate analytical solutions to singular perturbation problems that cannot be solved exactly. Our objective is to help young and also established scientists and engineers to build the skills necessary to analyze equations that they encounter in their work. Our presentation is aimed at developing the insights and techniques that are most useful for attacking new problems. A variety of singularly perturbed problems arising in physical problems are also mentioned. However the matter presented in this paper is available in different books and research articles but we summarized the important useful material in an effective manner which can serve as an introduction to new researchers. The organization of this paper is as follows Section 2: Singular perturbation problems are explained by considering examples of different variety. Section 3: Methods of solving singular perturbation problems are described in detail and the recent development in the existing literature has been mentioned. Section 4: Some worked out examples of singular perturbation problems arising in physical problems are given and elaborated properly. Section 5: Comparative study of the methods with further direction of research and concluding remark. 2. Regular and singular perturbation theory Perturbation theory also occurs in two varieties. We define a regular perturbation problem as one whose perturbation series is a power series in ε having a nonvanishing radius of convergence. A basic feature of all regular perturbation problems is that the exact solution for small but nonzero |ε| smoothly approaches the unperturbed or zeroth-order solution as ε → 0. We define a singular perturbation problem as one whose perturbation series either does not take the form of a power series or, if it does, the power series has a vanishing radius of convergence. In singular perturbation theory there is sometimes no solution to the unperturbed problem i.e. the exact solution as a function of ε may cease to exist when ε = 0, when a solution to the unperturbed problem does exist, its qualitative features are distinctly different from those of the exact solution for arbitrarily small but nonzero ε . In either case, the exact solution for ε = 0 is fundamentally different in character from the ‘‘neighboring’’ solutions obtained in the limit ε → 0. If there is no such abrupt change in character, then we would have to classify the problem as a regular perturbation problem. When dealing with a singular perturbation problem, one must take care to distinguish between the zeroth-order solution (the leading term in the perturbation series) and the solution of the unperturbed problem, since the latter may not even exist. There is no difference between these two in a regular perturbation theory, but in a singular perturbation theory the zeroth-order solution may depend on ε and may exist only for nonzero ε . Here are some examples of regular and singular perturbation problems. Example 2.1. The initial value problem y′′ + (1 − ε x)y = 0,

y(0) = 1,

y′ (0) = 0,

(1)

is a regular perturbation problem in ε over the finite interval 0 ≤ x ≤ L. In fact, the perturbation solution converges for all x and ε , with increasing rapidity as ε → 0+ for fixed x. It is a singular perturbation problem on the infinite interval 0 ≤ x ≤ ∞ because the perturbed solution (ε > 0) is not close to the unperturbed solution (ε = 0). Example 2.2. The nonlinear oscillator problem y′′ + y + ε y3 = 0,

y(0) = 1,

y′ (0) = 0

(2)

is a regular perturbation problem because when ε → 0, no change in the order of differential equation and the exact solution for small but nonzero |ε| smoothly approaches the unperturbed or zeroth-order solution as ε → 0. Example 2.3. The boundary value problem

ε y′′ − y′ = 0,

y(0) = 0,

y(1) = 1,

(3)

is a singular perturbation problem because the associated unperturbed problem

− y′ = 0,

y(0) = 0,

y(1) = 1,

(4)

has no solution. The solution to this first order differential equation, y = constant, cannot satisfy both boundary conditions.

558

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

Example 2.4. The boundary value problem

ε y′′ + y = 0,

y(0) = 0,

y(1) = 1,

(5)

is a singular perturbation problem because when ε = 0, the solution to the unperturbed problem, y = 0, does not satisfy the boundary condition y(1) = 1. 3. Methods for solving singular perturbative problems In this section, we describe the general methods which are used to solve the singular perturbation problems in the existing literature. 3.1. Method of asymptotic expansion In mathematics, an asymptotic expansion is a series expansion of functions which may converge or diverge but whose partial sum can be made an arbitrary good approximation to a given function. Typically, the best approximation is obtained when the series is truncated at the smallest term. The convergent series expansions are asymptotic, while asymptotic series expansions are, in general, divergent. Asymptotic expansions are used in analysis to describe the behavior of a function in a limiting situation. To explain the asymptotic expansion, we consider the following. Let a set of functions {φn (x)}, n = 0, 1, 2, . . . is an asymptotic sequence as x → a if φn+1 (x) = o(φn (x)) as x → a for every n ≥ 0, where a is the limit point. For example, in sequence {xn }, as x → 0 and {ln[1 + (1 − x)n ]} as x → 1both sequences conform that limx→0 [φn+1 (x)/(φn (x))] = 0 and limx→1 [φn+1 (x)/(φn (x))] = 0 for every n ≥ 0. The series of terms written as N −

An φn (x) + O(φn+1 )

(6)

n =0

where the An is constant, and is an asymptotic expansion of the function f (x), with respect to the asymptotic sequence {φn (x)} if, for every N ≥ 0, f (x) −

N −

An φn (x) = o[φN (x)]

as x → a.

(7)

n =0

If this expansion exists, it is unique in that the coefficients, An , are completely determined. The most common type of asymptotic expansion is a power series in either positive or negative powers. Asymptotic expansions arise in the approximation of integrals—Laplace method, Watson’s lemma and the method of steepest descents. Integration by parts is a particularly easy procedure for developing asymptotic approximations to many kinds of integrals. The asymptotic expansion of integral representations is an extremely important technique because the entire special functions (Bessel, Airy, Gamma, Parabolic cylinder, Hyper geometric, Error function) and Exponential integral commonly used in mathematical physics and applied mathematics has integral representations. Some examples of asymptotic expansions are given as Gamma function :

ex

xx

√

2π x

0 (x + 1) ∼ 1 +

Exponential integral : xex E1 (x) ∼



12x

∞ − (−1)n n!

xn

n =0

Error function :

1

π xex2 erf c (x) ∼ 1 +

+

1 288x2

−

139 51 840x3

− ···

as x → ∞

as x → ∞

∞ − (2n)! (−1)n n!(2x)2n n=1

as x → ∞.

To illustrate the method of asymptotic expansion in detail we consider the following example

ε2 + 1+ 2 dx x + ε2

dy





y + ε y2 = 0,

0 ≤ x ≤ 1, with y(1; ε) = 1.

(8)

As ε → 0, we are interested in finding the asymptotic solution of (8) in the form y(x; ε) ∼

N − n=0

ε n y n ( x)

(9)

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

559

and we will need to find at least three terms y0 , y1 and y2 . From (9) and (8) equations for the yn (x) are y0

+ 2y0 y1 = 0, (10) x2 and so on, the boundary condition requires that y0 (1) = 1, yn (1) = 0 for n ≥ 1. Here, we should expect that evaluation of the expansion on x = 1 is allowed i.e. all terms are defined for x = O(1), but we must anticipate difficulties as x → 0. The solutions for the functions y0 (x) and y1 (x) are y′0 + y0 = 0,

y′1 + y1 + y20 = 0,

y0 (x) = e1−x ,

y′2 + y2 +

y1 (x) = e2−2x − e1−x .

(11)

The 2-term asymptotic expansion, y ∼ y0 + ε y1 is uniformly valid for ∀x ∈ [0, 1]. To find the next term in the expansion, we can write a third equation of (10) from (11) e1−x

+ 2e1−x (e2−2x − e1−x ) = 0 with y2 (1) = 0. x2 On introducing the integrating factor ex , we have y′2 + y2 +

(ex y2 )′ +

e

e1−x

+ e1−x (e2−2x − 2e1−x ) + Ae−x . (13) x To satisfy y2 (1) = 0 the arbitrary constant A must be zero. This third term in the asymptotic expansion is very different from the first two, as it is not defined on x = 0, so we must expect a breakdown. The expansion, to this order, is now x2

y(x, ε) = e

+ 2e(e2−2x − e1−x ) = 0 and so y2 (x) =

(12)

+ ε(e

1−x

2−2x

−e

1−x

)+ε

2



e1−x x

+e

3−3x

2−2x



− 2e

(14)

as ε → 0 for x = O(1); as x → 0, we clearly have a breakdown where the second and third terms in the expansion become the same size i.e. ε = O(ε 2 /x) or x = O(ε). Note that this breakdown occurs for a larger size of x (as x is decreased from O(1)) than the breakdown associated with the first and third terms, so we must consider x = O(ε). This example for x = O(ε) is formulated by writing x = εX

y(ε X ; ε) ≡ Y (X ; ε),

and

(y = O(1) for x = O(ε)).

The original Eq. (8), expressed in terms of X and Y , requires the identity dy

=

dY

=

d

dx dx dx and then we obtain

ε −1

dY dX

 + 1+

dY Y (x/ε; ε) = ε −1 dX

ε2 ε2 X 2 + ε2



Y + εY 2 = 0

or

dY dX

 +ε 1+

1



1 + X2

Y + ε 2 Y 2 = 0,

(15)

but the boundary condition is not available, because this is specified where x = O(1). Eq. (15) suggests that we seek a solution in the form Y (X ; ε) ∼

N −

ε n Yn (X )

(16)

n =0

which gives Y0′ = 0;

Y1′ +

 1+

1 1 + X2



Y0 = 0,

and so on.

(17)



Immediately we obtain Y0 = A0 which is an arbitrary constant, and then second Eq. (17) becomes Y1′ = − 1 + 1+1X 2 which integrates to give Y1 (X ) = −A0 (X + arctan X ) + A1



A0

(18)

where A1 is a second arbitrary constant. The resulting 2-term expansion is therefore Y (X ; ε) ∼ A0 + ε{A1 − A0 (X + arctan X )},

x = O(1)

(19)

the two arbitrary constants are determined by invoking the matching principle, (19) and (14) are to match. Thus we write the terms in (14) as functions of X , let ε → 0 (for x = O(1)) and retain terms O(1) and O(ε) (which are used in (19)); conversely, write (19) as a function of x, expand and retain terms O(1), O(ε) and O(ε 2 ). From (14) we construct e1−εX + ε(e2−2εX − e1−εX ) + ε 2



e1−εX

εX

+ e3−3εX − 2e2−2εX



∼e+ε

e X

+ e2 − e − eX



(20)

560

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

for x = O(1); and from (19) we write



A0 + ε A1 − A0

x ε

+ arctan

x 

ε

 π  ε 2 A0 ∼ A0 − A0 x + ε A1 − A0 + 2

x

for x = O(1).

(21)

The Eq. (21) requires the following formula arctan X ∼

π 2

−

1 X

+

1 3X 3

as X → +∞.

(22)

Eqs. (20) and (21), match if we consider A0 = e and A1 = e2 − e + e2π ; then Eq. (19) for x = O(1) can be written as Y (X , ε) ∼ e + ε

 eπ 2

 + e2 − e − e (X + arctan X ) .

(23)

From the above we observe that, although the expansion (14) is not defined on x = 0, the expansion valid for x = O(ε) does allow evaluation on x = 0 i.e. X = 0; in fact, from (23), we obtain



y(0; ε) = Y (0; ε) ∼ e 1 + ε

π 2

 +e−1

as ε → 0.

(24)

As a conclusion, the procedure involves the construction of an asymptotic expansion valid for x = O(1) and applying the boundary condition(s) if the expansion remains valid here. The expansion is then examined for ∀x ∈ D, seeking any breakdowns, rescaling and hence rewriting the equation in terms of the new, scaled variable; this example is then solved as another asymptotic expansion, matching as necessary. In Refs. [20–34] the method of asymptotic expansion has been used to solve the singular perturbation problems of different kinds. 3.2. Method of matched asymptotic expansion Asymptotic matching is mostly used to determine a uniform and accurate approximation to the solution of a singularly perturbed differential equations and to find global properties of differential equations such as eigen values. In a large class of singularly perturbed boundary value problems, the interval may be divided into two overlapping subintervals namely an inner region and an outer region. Then, each subinterval is used to obtain an asymptotic approximation to the solution of differential equation valid on the interval. Finally, the matching is done by assuring that the asymptotic approximations have the same functional form on the overlap of every pair of intervals and that it satisfies all the boundary conditions given at various points on the interval. This method is a generalization of Prandtl’s boundary layer theory. The outer and inner solutions are then combined through a process called ‘‘matching’’ in such a way that an approximate solution for the whole domain is obtained. To explain this method we consider the following equation.

ε y′′ + (1 + ε)y′ + y = 0,

y(0) = 0,

y(1) = 1,

and

0 < ε ≪ 1.

(25)

Outer and inner solutions. Since ε is very small, our first approach is to find the solution to the problem y′ + y = 0, which is yo = Ae−t where A is any constant. Applying the boundary condition y(0) = 0, we would have A = 0 and applying the boundary condition y(1) = 1, we would have A = e, i.e. one of the boundary conditions cannot be satisfied. From this we infer that there must be a boundary layer at one of the endpoints of the domain. Suppose the boundary layer is at t = 0. If we make the scaling τ = t /ε, the problem becomes 1 ′′ 1 y (τ ) + (1 + ε) y′ (τ ) + y(τ ) = 0.

ε

(26)

yI = B(1 − e−τ ) = B(1 − e−t /ε ).

(27)

ε Multiplying Eq. (26) by ε and taking ε = 0 gives y′′ + y′ = 0 with solution y = B − C e−τ for some constants B and C . Since y(0) = 0, we have C = B, so the inner solution is

Matching. Remember that we have assumed the outer solution to be y0 = e1−t . The idea of matching is for the√ inner and outer solutions to agree at some value of t near the boundary layer as ε decreases. For example, if we fix t = ε , we have the matching condition as

√

√

lim yI ( ε) = lim y0 ( ε) so B = e.

ε→0

ε→0

(28)

To obtain our final matched solution, valid on the whole domain, we add the inner and outer approximations and subtract their common value y(t ) = yI + y0 − e = e(1 − e−t /ε ) + e1−t − e = e(e−t − e−t /ε ).

(29)

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

561

Substituting the matched solution in the problem (25) a differential equation yields

ε y′′ + (1 + ε)y′ + y = ε e1−t

(30) 1−1/ε

which converges to zero as ε decreases, uniformly in t. As to the boundary conditions, y(0) = 0 and y(1) = 1 − e obtained from (29), which quickly converges to the value given in the Eq. (25). Our solution approximately solves the problem (25) at hand and at the same time closely approximates the problem’s exact solution. It can be seen that this particular problem is easily found to have an exact solution given by e−t − e−t /ε y(t ) = −1 , e − e−1/ε which is very close to our approximation for reasonably small ε .

(31)

Location of boundary layer. Conveniently, we can see that the boundary layer, where y′ and y′′ are large, is near t = 0, as we supposed earlier. If (1−t ) we had supposed it to be at the other endpoint and proceeded by making the rescaling τ = ε , we would have found it impossible to satisfy the resulting matching condition. For different problems, this kind of trial and error is the only way to determine the true location of the boundary layer a priori. For a detailed study we refer to [35,16,36–41] where the method of matched asymptotic has been applied to solve a variety of singular perturbation problems. 3.3. Method of multiple-scale analysis The method of multiple scales comprises techniques used to construct uniformly valid approximations to the solutions of perturbation problems, both for small as well as large values of the independent variables. This is done by introducing fastscale and slow-scale variables for an independent variable. Thereafter in the solution process of the perturbation problem, the resulting additional freedoms, introduced by the new independent variables are used to remove unwanted secular terms. In general, secular terms always appear whenever the inhomogeneous term is itself a solution of the associated homogeneous constant coefficient differential equation. In 1882–1883 Lindstedt introduced this method and used such scaled variables to eliminate spurious secular (resonant) terms in perturbation expansions in celestial mechanics. As an example to illustrate the method of multiple-scale analysis, we consider the undamped and unforced Duffing equation which is the regular perturbation problem given in [41]. d2 y

+ y + ε y3 = 0, y(0) = 1, y′ (0) = 0, (32) dt 2 which is a second order ordinary differential equation describing a nonlinear oscillator. A solution y(t ) is sought for small values of the positive nonlinearity parameter 0 < ε ≪ 1. First eliminating the most secular terms to all orders begins by introducing a new variable τ = ε t. τ defines a long time scale because τ is not negligible when t is of order 1/ε or larger. Even though the exact solution y(t ) is a function of t alone, multiple-scale analysis seeks solutions which are functions of both variables t and τ treated as independent variables. We emphasize that expressing y as a function of two variables is an artifice to remove secular effects, the actual solution has t and τ related by τ = ε t so that t and τ are ultimately not independent. The formal procedure consists of assuming a perturbation expansion of the form y(t ) = Y0 (t , τ ) + ε Y1 (t , τ ) + · · · .

(33)

We use the chain rule for partial differentiation to compute a derivative of y(t ) dy dt

 =

∂ Y0 ∂ Y0 dτ + ∂t ∂τ dt



+ε



∂ Y1 ∂ Y1 dτ + ∂t ∂τ dt



+ ···.

(34)

However, since τ = ε t , ddtτ = ε . Thus,

  ∂ Y0 ∂ Y0 ∂ Y1 = +ε + + O(ε 2 ). dt ∂t ∂τ ∂t

dy

(35)

Also, differentiating with respect to t again gives

 2  ∂ 2 Y0 ∂ Y0 ∂ 2 Y1 = +ε 2 + + O(ε2 ). dt 2 ∂t2 ∂τ ∂ t ∂t2

d2 y

(36)

Substituting (36) into (32) and collecting powers of ε gives

∂ 2 Y0 + Y0 = 0, ∂t2 ∂ 2 Y1 ∂ 2 Y0 3 + Y = − Y − 2 . 1 0 ∂t2 ∂τ ∂ t

(37) (38)

562

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

The general solution of (37) is Y0 (t , τ ) = A(τ )e+it + A∗ (τ )e−it ,

(39)

where A(τ ) is an arbitrary complex function of τ . A(τ ) will be determined by the condition that secular terms do not appear in the solution to (38). From (39), the right-hand side of (38) is

[

]

dA

e+it −3A2 A∗ − 2i

dτ

[ ] dA∗ + e−it −3A(A∗ )2 + 2i − e+3it A3 − e−3it (A∗ )3 . dτ

(40)

To preclude the appearance of secularity, we require that the as yet arbitrary function A(τ ) satisfy

−3A2 A∗ − 2i

dA dτ

−3A(A∗ )2 + 2i

= 0,

dA∗

(41)

= 0.

dτ

(42)

To solve the Eq. (41) for A(τ ), we represent A(τ ) in polar coordinate form as A(τ ) = R(τ )eiθ(τ )

(43)

where R(τ ) and θ (τ ) are real. Also, differentiating with respect to τ gives dA dτ

=

dR iθ dθ i θ e + Ri e . dτ

(44)

dτ

Substituting into Eq. (41) and equating real and imaginary parts gives dR dτ

dθ

= 0 and

dτ

=

3 2

R2 .

2 Therefore, A(τ ) = R(0)eiθ(0)+3iR (0)τ /2

(45)

and the zeroth order solution (equating real part) (39) is Y0 (t , τ ) = 2R(0) cos[θ (0) +

3 2

R2 (0)τ + t ]

(46)

′ R(τ ) and θ (τ ) are determined by the initial condition   y(0) = 1, y (0) = 0.

These conditions imply that Y0 (0, 0) = 1,

θ(0) = 0. Therefore the zeroth order solution is ] [ 3 Y0 (t , τ ) = cos t + τ .

∂ Y0 ∂t

(0, 0) = 0, putting these values in Eq. (46) we get R(0) =

8

1 2

and

(47)

Finally, since τ = ε t, By Eq. (33) we get y(t ) = cos

[ 1+

3 8

 ]

ε t + O(ε),

ε → 0+, ε t = O(1).

(48)

Higher-order solutions, using the method of multiple scales, require the introduction of additional slow scales, i.e. t2 = ε2 t , t3 = ε 3 t etc. However, this introduces possible ambiguities in the perturbation series solution, which require a careful treatment as given in [14]. In Refs. [42,43,41,44–50] multiple scale method in different forms for solving singular perturbation problems is given. 3.4. Method of WKB approximation The WKB theory is a powerful tool for obtaining a global approximation to the solution of a linear singularly perturbed differential equation of any order, eigen value problems, initial and boundary value problems and is also used to evaluate integrals of the solution of a differential equation. The most familiar example of a semi classical calculation in quantum mechanics in which the wave function is recast as an exponential function, semi classically expanded, and then either the amplitude or the phase is taken to be slowly changing. The limitation of this technique is that they are only useful for linear equations. The name of this method is an acronym for the Wentzel–Kramers–Brillouin approximation. Other often used acronyms for the method include the JWKB approximation and the WKBJ approximation, where the ‘‘J’’ stands for Jeffreys.

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

563

Generally, WKB theory is a method for approximating the solution of a differential equation whose highest derivative is multiplied by a small parameter ε . In this method we assume a solution of the form of an asymptotic series expansion as

 y(x) ∼ exp

∞ 1−

δ

 δ S n ( x)

as δ → 0

n

(49)

n =0

where S (x) is a phase function and depends implicitly on δ . This expression is the starting formula for which all WKB approximations are derived. WKB theory is a special case of multiple scale analysis [14,41,51]. As an example, consider that the second-order homogeneous linear differential equation is in Schrodinger form if the y′ term is absent.

ε2

d2 y dx2

= Q (x)y where Q (x) ̸= 0.

(50)

Differentiating (49) twice we have

 y ( x) ∼



∞ 1−



∞ 1−



δ (x) exp δ Sn (x) as δ → 0 δ n =0 δ n=0     2  ∞ ∞ ∞ − − − 1 1 1 δ n Sn′ + δ n Sn′′  exp δ n Sn (x) y′′ (x) ∼  2 δ δ n=0 δ n =0 n=0 ′

n ′ Sn

n

(51)

as δ → 0.

(52)

Now, we substitute the values of y(x) and y′′ (x) into (50) and dividing by the exponential factors, we obtain

ε 2 ′2 2ε 2 ′ ′ ε 2 ′′ S S S + · · · = Q (x) as δ → 0. + S + 0 0 1 δ2 δ δ 0

(53)

The largest term of the left side of (53) is εδ 2 S0′2 . By dominant balance this term must have the same order of magnitude as Q (x) on the right side. Thus δ is proportional to ε and for simplicity we choose δ = ε in (53) then comparing powers of ε , 2

ε 0 : S0′2 = Q (x).

(54)

The equation for S0 is called the eikonal equation, its solution is S 0 ( x) = ±

∫ x

Q (t ) dt .

(55)

a

Looking at first-order powers of ε gives ε 1 : 2S0′ S1′ + S0′′ = 0.

(56)

The equation for S1 is called the transport equation, its solution is S 1 ( x) = −

1 4

ln(Q (x)).

(57)

Combining (55) and (57) gives a pair of approximate solutions to the Schrodinger Eq. (50) one for each sign of S0 . The general solution is a linear combination of the two equations 1

y(x) ∼ c1 Q − 4 (x) exp

x

[ ∫ 1

ε

dt



]

1

[

Q (t ) + c2 Q − 4 (x) exp −

a

1

ε

x

∫

dt



Q (t )

]

as ε → 0

(58)

a

where c1 and c2 are constants to be determined from initial or boundary conditions and a is an arbitrary but fixed integration point. This expression is the leading order WKB approximations to the solution of (50), it differs from the exact solution by terms of order ε in regions where Q (x) ̸= 0. A more accurate approximation to y(x) may be constructed from the higher-order terms in the WKB series. Explicitly 2S0′ Sn′ + Sn′′−1 +

n −1 −

Sj′ Sn′ −j = 0

for n ≥ 2.

j =1

In Refs. [52–56,51,57] WKB methods have been described for different types of singular perturbation problems.

(59)

564

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

3.5. Poincare–Lindstedt method In perturbation theory, the Poincare–Lindstedt method is a technique for uniformly approximating periodic solutions to ordinary differential equations, when regular perturbation approaches fail. The method removes secular terms (terms growing without bound) arising in the straightforward application of perturbation theory to weakly nonlinear problems with finite oscillatory solutions [58]. The method is named after Henri Poincare [59] and Anders Lindstedt [60]. The undamped, unforced Duffing equation [61] is given by x¨ + x + ε x3 = 0

for t > 0, 0 < ε ≪ 1

and

x(0) = 1,

x˙ (0) = 0.

(60)

A perturbation-series solution of the form x(t ) = x0 (t ) + ε x1 (t ) + · · · is sought. The first two terms of the series are x(t ) = cos(t ) + ε

[

1 32

] 3 (cos(3t ) − cos(t )) − t sin(t ) + · · · . 8

(61)

This approximation grows without bound in time, which is inconsistent with the physical system that the equation models. The term responsible for this unbounded growth, called the secular term, is t sin t. The Poincare–Lindstedt method allows for the creation of an approximation that is accurate for all time, as follows. In addition to expressing the solution itself as an asymptotic series, form another series with which to scale time t : τ = ωt where ω = ω0 + εω1 + · · ·. For convenience, take ω0 = 1 because the leading order of the solution’s angular frequency is 1. Then the original problem becomes

ω2 x′′ (τ ) + x(τ ) + ε x3 (τ ) = 0

(62)

with the same initial conditions. Now search for a solution of the form x(τ ) = x0 (τ ) + ε x1 (τ ) + · · · .

(63)

The following solutions for the zeroth and first order problem in ε are obtained x0 = cos(τ ),

x1 =

1 32

cos(3τ ) +

  3 ω1 − τ sin(τ ). 8

(64)

So the secular term can be removed through the choice ω1 = 3/8. Higher orders of accuracy can be obtained by continuing the perturbation analysis along this way. As of now, the approximation, correct up to first order in ε is x(t ) ≈ cos

 1+

3 8

      1 3 ε t + ε cos 3 1 + ε t . 32

8

(65)

In Refs. [62–67,58–61,68–77] the Poincare–Lindstedt method has been described for different types of singular perturbation problems. 3.6. Periodic averaging method In the study of dynamical system, the method of periodic averaging is used to study certain time-varying systems by analyzing easier, time-invariant systems obtained by averaging the original system. Consider a general, nonlinear dynamical system x˙ = ε f (t , x, ε)

(66)

where f (t , x) is periodic in t with period T . The evolution of this system is said to occur in two timescales: a fast oscillatory one associated with the presence of t in f and a slow one associated with the presence of ε in front of f . Assume, for simplicity, that f is periodic in t with period T . The corresponding (leading order in ε ) averaged system is x˙ a = ε

1 T

∫

T

f (τ , x, 0) dτ = f˜ (xa ).

(67)

0

Averaging mods out the fast oscillatory dynamics by averaging their effect through time integration in Eq. (67). In this way, the mean (or long-term) behavior of the system is retained in the form of the dynamical equation for the evolution for xa . Standard methods for time-invariant systems may then be employed to analyze the equilibria (and their stability) as well as other dynamical objects of interest present in the phase space of the averaged system. In Refs. [78–90] the Periodic averaging method has been described for different types of singular perturbation problems.

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

565

4. Some worked examples arising from physical problems The aim of this section is to present a number of worked examples taken from [18] that introduce, describe, develop, explain and solve practical problems in various fields arising in physical problems. Here, we will not dwell upon the purely technical aspects of finding the solution of a particular differential equation. Each problem is described with sufficient detail to enable it to be put into context, although it would be quite impossible to include all the background ideas for those altogether unfamiliar with the particular field. It is the effort to motivate researcher to know about the power of the technique which has versatility in various fields. 4.1. Mechanical and electrical systems Few examples of singular perturbation problems arising in mechanical and electrical systems are collected which are based on fairly simple mechanical or physical principles. 4.1.1. Projectile motion with small drag We consider a projectile which is moving in a two dimensional (x, z ) plane under the action of gravity (which is constant in the negative z-direction) and of a drag force proportional to the square of the speed (and acting back along the local direction of motion). The non-dimensional equations are most conveniently written as du

 = −ε u u2 + w 2 ,

dt

dw dt

 = −1 − εw u2 + w 2

(68)

d where (u, w) = dt (x, z ) and the initial conditions are given as u = cos α, w = sin α, x = z = 0 ∀t = 0 and α is the angle of projection (0 < α < π /2). The small parameter is ε (>0), and for motion such as that for a shot-put its value is typically about 0.01. In these projectile problems, the main interest is in estimating the range and maximum range of x for a given vertical displacement (z), which may be zero, here we find x where z = h.

4.1.2. Child’s swing In this example of real life, we can again feel our childhood memories with swings in parks or in our houses in a rainy time. We never put a thought that it could also be interesting to understand that particular fun as mathematical equation. It is a very simple model to see the effect on amplitude and frequency with a change in a small amount of length. We are all familiar with the child’s swing and the technique for increasing the arc (amplitude) of the swing. The process of swinging the legs (coupled with a small movement of the torso) causes the center of gravity of the body to be raised and lowered periodically. This can be modeled by treating the swing as a pendulum which changes its length l(t ) by a small amount. The model equation for this (in the absence of damping) is d2 θ dt 2

+

2 dl dθ l dt dt

+ sin θ = 0,

t ≥0

(69)

for θ(t ) the angle of the swing, given l(t ). We choose to represent the child’s movement on the swing by l = l0 (1 + ε sin ωt ) where l0 is a positive constant and ω is a constant frequency to be selected. Thus we approximate the Eq. (69) as d2 θ dt 2

+

(2εω cos wt ) dθ + θ = 0. (1 + ε sin ωt ) dt

(70)

To solve the above equation using the method of multiple scales, we take the fast scale as T = t (or more generally T = Ω t but there is no advantage in this, for we are led to the choice Ω = 1) and by virtue of the term in ddtθ , a slow scale τ = ε t, thus we have the identity

∂ ∂ +ε . (71) dt ∂T ∂τ Now Eq. (70), with θ (t ; ε) ≡ Θ (T , τ ; ε) becomes   2ω cos ωt ΘTT + 2ε ΘT τ + ε 2 Θτ τ + ε (72) (ΘT + εΘτ ) + Θ = 0 1 + ε sin ωt ∑∞ n and we seek a solution in the form Θ (T , τ ; ε) ∼ n=0 ε Θn (T , τ ) which is periodic in T . In this problem, we have yet to choose the frequency ω which the child can control in order to increase the amplitude of the swing. Without any damping in the model, we must anticipate that we can find an ε which allows the amplitude to grow without bound. We conclude d

≡

therefore that the adjustments provided by the child on the swing must be at twice the frequency of the oscillation of the swing, which is what we learnt as children.

566

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

4.1.3. Meniscus on a circular tube At school level it was a new phenomenon to observe a liquid rising in a small-diameter tube that penetrates (vertically) the surface of the liquid and is very familiar, as are the menisci that form inside and outside the tube. The basic model assumes that the mean curvature at the surface is proportional to the pressure difference across the surface. With the two principal curvatures of radii written as k1 and k2 then this assumption can be expressed as k1 + k2 ∝ z where z is the vertical coordinate and the pressure difference is proportional to the (local) height of the liquid above the undisturbed level far away from the tube. This relation is usually referred to as Laplace’s formula. In detail, written in non-dimensional form, this equation (for cylindrical symmetry) becomes 1 d

 

r dr



r (dz /dr ) 1 + (dz /dr )2

= εz

(73)

where the surface is z = z (r , ε) and r is the radial coordinate with r = 0 at the center of the tube and the tube wall (of infinitesimal thickness) is at r = 1, the liquid surface satisfies z → 0 as r → ∞. The non-dimensional parameter ε is inversely proportional to the surface tension in the liquid and proportional to the square of the tube’s radius. We will examine the problem of solving Eq. (73) for ε → 0+ with the boundary conditions dz = 0 at r = 0, z → 0 as r → ∞, dz = ∓ cot θ dr dr ± as r → 1 (i.e. both inside and outside) where θ is the given contact angle between the meniscus and the tube (measured relative to the upward vertical side of the tube). As we will see, this problem results in the construction of a uniformly valid expansion in the interior and a scaling and matching problem in the exterior. Interior problem. It is a familiar observation, at least for 0 < θ < π /2 that, the narrower the tube the higher the liquid rises in the tube. This suggests that the height of the liquid will increase as ε → 0. When this is coupled with the property (the confirmation of which is left as an exercise) that no relevant solution of Eq. (73) exists if we ignore the term ε z we are led to write the solution in the form z (r , ε) = ε −1 h(ε) + ζ (r , ε),

where h(0) = O(1), h > 0.

(74)

Now we seek an asymptotic solution h(ε) ∼

∞ −

ε n hn ,

ζ (r , ε) ∼

∞ −

ε n ζn ( r )

(75)

n =0

n =0

and so the leading-order problem becomes 1 d

 

r dr



r (dζ0 /dr ) 1 + (dζ0 /dr )2

= h0 with ζ0 (0) = 0,

ζ0′ (0) = 0,

ζ0′ → cot θ as r → 1− .

(76)

The prime denotes the derivative with respect to r. Exterior problem. Finding the solution for the shape of the surface outside the tube is technically a more demanding exercise. First, the deviation of the surface from its level (z = 0) at infinity is observed to be not particularly large, as compared with what happens inside. This suggests that we attempt to solve Eq. (73) directly, subject to z ′ (r , ε) → − cot θ as r → 1+ z → 0 and z ′ → 0 as r → ∞. Let us write z (r , ε) = z0 (r ) + O(1) as ε → 0+ . So that we are not committing ourselves, at this stage, to the size of the next term in the asymptotic expansion. In fact, we shall observe that logarithmic terms arise, although we will not pursue the details here. The equation for z0 (r ) is simply 1 d r dr



r (z0′ )



1 + (z0′ )2



 = 0 or

r (z0′ )





1 + (z0′ )2

=A

(77)

where the arbitrary constant A is determined as A = − cos θ . On further integration it yields z0 = − cos θ [ln(r sec θ +

√

r 2 sec2 θ − 1)] + B

(78)

where B is a second arbitrary constant. The complications alluded to earlier are now evident z0 (r ) ∼ − cos θ ln r as r → ∞, which can never accommodate the conditions at infinity, for any choice of B. 4.1.4. Drilling by laser Now we are dealing with an engineering application problem which has wide use in electronics and mechanical fields. This problem is a one-dimensional model for the process of drilling through a thick block of material using a laser. The laser heats the material until it vaporizes and we assume that the vapor is continuously removed. The essential character of this

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

567

problem is therefore one of heat transfer at the boundary and the bottom of the drill hole, which is moving. In suitable non-dimensional variables, the temperature relative to ambient conditions satisfies the classical heat conduction equation

∂T ∂ 2T = 2, ∂t ∂x

t > 0, 0 < x < ∞,

T (x, 0, ε) = 0,

T (x, t , ε) → 0

as x → ∞

(79)

which describes the initial state and the condition at infinity, respectively. T (X (t , ε), t , ε) = 1.

(80)

Eq. (80) is the vaporization condition at the bottom of the drill hole at x = X (t , ε). The speed of the drilling process is controlled by dX

=1+ε

∂T (X (t , ε), t , ε) with X (0, ε) = 0, ∂x

(81) dt which is based on Fourier’s law applied at the bottom of the hole. This set (79)–(81) is an example of a Stefan problem. This problem perhaps more than the previous three, shows how powerful these techniques can be in illuminating the details. 4.1.5. The van der Pol (Rayleigh) oscillator This example is taken from the electronics field and it requires a fairly routine application of the method of multiple scales to a nearly linear oscillator, although the solution that we obtain takes quite a dramatic form. The equation first came to prominence following the work of van der Pol in 1922 on the self-sustaining oscillations of a triode circuit. However, essentially the same equation had already been discussed by Rayleigh in 1883, as a model for ‘maintained’ vibrations in, for example, organ pipes. We will write the equation in the Rayleigh form x¨ + x − ε

 x˙ −

1 3

3



x˙

= 0,

ε > 0, t ≥ 0

(82)

for x(t , ε) in the context of the van der Pol problem, x˙ is proportional to the grid voltage V. The governing equations for this triode circuit are di dia dV + iR + V − M = 0, C = i, dt dt dt where µ and λ are positive constants. L

ia = λV (1 − µV 2 ),

(83)

4.1.6. A diode oscillator with a current pump This problem has an electrical background but with a different asymptotic structure containing two small parameters, we seek the initial condition which leads to a periodic solution. The circuit is represented by the equations i=C

dV dt

,

i1 = (eα V − 1)Is ,

and then Kirchhoff’s law gives

i2 = (1 − e−α V )Is C

dV dt

(84)

+ (eαV − e−αV )Is = I sin ωt

(85)

which leads to the non-dimensional equation (x ∝ eα V )

δ˙x= x sin t − x2 + ε,

x(0, ε) = a (1 > a > 0).

(86)

Typical values of the parameters are δ ≈ 0.03 and ε ≈ 10−10 . 4.1.7. Klein–Gordon equation The general form of the Klein–Gordon equation, written in one spatial dimension, is utt − uxx + V ′ (u) = 0,

(87)

where V (u) (which can be taken as a potential, in quantum-mechanical terms) is typically a function with nonlinearity more severe than quadratic. We will consider the problem for which V ′ (u) = u + ε u3 ,

ε > 0,

(88)

and so introduce a parameter which we will allow one to satisfy ε → 0 . The equation with this choice arises, for example, in the study of wave propagation in cold plasma. (The choice V (u) = − cos u gives rise to the so called sine-Gordon equation, a pun on the original name and this equation has exact ‘soliton’ solutions.) +

4.2. Celestial mechanics Here, we described the celestial mechanics problems which are modeled in the form of singular perturbation problems. We present three typical problems that arise from planetary or related motions.

568

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

4.2.1. The Einstein equation for Mercury In childhood it was very fascinating to see stars in the sky. As we went to higher classes, we learnt more about the stars and universe. According to our forefathers every planet is considered as a God of some skills. For example Mercury (budh) is known as the God of intelligence but by the end of school education we knew that it is simply a physical body like our Earth with some other physical properties. Now here we study the motion of planet. Classical Newtonian (Keplerian) mechanics leads to an equation for a single planet around a sun of the form d2 u dθ 2

+ u = h (> 0, constant)

(89)

where θ is the polar angle of the orbit, u is inversely proportional to the radial coordinate of the orbit and h measures the angular momentum of the planet. However, when a correction based on Einstein’s theory of gravitation is added, the equation becomes d2 u dθ 2

+ u = h + ε u2

(90)

where ε (>0) is a small parameter (about 10−7 for Mercury, the planet for which the equation was first introduced). This example has proved to be a particularly straightforward application of the method of multiple scales. 4.2.2. Planetary rings In a study of a model for differentially rotating discs, the radial structure of the azimuthal velocity component for a large azimuthal mode number (essentially ε −1 here) satisfies an equation of the form

ε

d dr

[ r

d dr

]   2  r (r v) = 1 − v, r0

0 < r < ∞.

(91)

The above equation clearly possesses a turning point at r = r0 . Let us examine the solution near this point first. For the neighborhood of r = r0 , we set r = r0 (1 + δ R) where δ(ε) → 0+ as ε → 0+ and ignore the scaling of v because the equation is linear, so Eq. (91) becomes

  ε d d (1 + δ R) [(1 + δ R)V ] = −δ R(2 + δ R)V , δ 2 dR dR

where v(r , ε) ≡ V (R, ε).

(92)

4.2.3. Slow decay of a satellite orbit The equations for a satellite in orbit around a primary mass (in the absence of all other masses), with a drag proportional to the (speed)2 are first written down in terms of polar coordinates (r , θ ). These are then transformed to r −1 = u(t ) and θ(t ), where t is time and finally, we obtain the non-dimensional equations d2 u dθ 2

+u=u

4



dt dθ

2

,

d dθ

 u

2

dt dθ



=ε

dt dθ



 u2 +

du

2

dθ

(93)

where ε (>0) is a measure of the drag coefficient on the satellite. This is Laplace’s important observation to u(θ ) and t (θ ). 4.3. Physics of particles and of light In this section, we will examine some problems that arise from fairly elementary physics. These will touch on quantum mechanics, light propagation and the movement of particles. 4.3.1. Perturbation of the bound states of Schrodinger’s equation Perturbation of the bound states of Schrodinger’s equation is a classical problem in elementary quantum mechanics. It involves the time-independent, one-dimensional Schrodinger equation

− ψ ′′ + V (x, ε)ψ = E ψ,

−∞ < x < ∞

(94)

where V (x, ε) is the given potential and E is the energy (i.e. the eigen-values of the differential equation). 4.3.2. Light propagating through a slowly varying medium Fermat’s principle states that light travels between any two points on a path which minimizes the time of propagation. If the path, in two dimensions is written as y = y(x) and the speed of light at any point is c (x, y) then y(x) must satisfy

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

c

d2 y dx2

 dy + cy − cx





1+

dx

dy

2 

dx

= 0.

569

(95)

This equation can be obtained either from the eikonal equation for rays or as the relevant Euler–Lagrange equation in the calculus of variations. In the special case where the medium varies only in x, so that c = c (x), we obtain the following equation γ ′ /c (x)  = constant. (96) 1 + γ ′2 If, for a light ray, we set γ ′ (x) = tan[α(x)], then we obtain the Snell’s law c (x) = constant. sin[α(x)] This problem of finding the path of a light ray has used the idea of multiple scales in a less routine way.

(97)

4.3.3. Raman scattering: a damped Morse oscillator Under certain circumstances, a small fraction of the incident light propagating through a medium may be scattered so that the wavelength of this light differs from that of the incident light, usually of greater wavelength. This is called Raman scattering. An example of this, which incorporates the Morse (exponential) model for the potential energy of atoms as a function of their separation, is x¨ + ε˙x + (1 − e−x ) e−x = 0,

t ≥ 0.

(98)

This equation also includes a (weak) linear damping term, which we will characterize by ε → 0+ . We wish to find a solution, subject to the initial conditions x(0, ε) = a, x˙ (0, ε) = 0 as Eq. (98) is nonlinear (although the underlying solution, valid for ε = 0, can be expressed in terms of elementary functions). 4.3.4. Quantum jumps: the ion trap In the study of the discontinuous emission or absorption of energy (quantum jumps), a single ion is trapped (in an electromagnetic device called a Paul trap) and its motion is governed by an equation of the form

ε(iφt + φxx ) = φv(x) cos(t /ε) (99) where ε > 0 is a parameter, v(x) is a given function (sufficient for the existence of φ ) and we seek the complex-valued function φ(x, t , ε) for ε → 0+ . In this case, we see that the oscillatory term on the right oscillates rapidly and use the method of multiple scale in the form T = t , τ = t /ε with φ(x, t , ε) ≡ Φ (x, T , τ , ε) which gives the equation iΦτ + ε(iΦT + Φxx ) = Φ v(x) cos τ . (100) 4.3.5. Low-pressure gas flow through a long tube The flow of a gas through a long circular tube i.e. small radius or length, where molecular collisions are assumed to occur only with the wall of the tube, can be represented by the Clausing integral equation N (x) = N0 (x) +

∫

L

K (x − γ )N (γ ) dγ .

(101)

0

Here, N (x)1x is the rate of molecular collisions (with the wall) between x and x + 1x, and N0 (x)1x the rate contributed by those molecules that have their first collision between these same stations. The kernel K (x − y) measures the probability that a molecule which has collided with the wall at x = y will collide again between the stations. This type of process is called a free-molecular or Knudsen flow. When we introduce the appropriate models for N0 (x) and K (x − y) non-dimensionalize and use the symmetry n(x) (i.e. n(x) + n(−x) = 0 so that n(0) = 0), we obtain the following equation



1 n(x, ε) = ε −1  2 1

 x+

 

− ε 1/

x+

1 2

+ ε2 −

2

 

1 2

+ ε 2 − 1/

2 −x  

1 2

 + ε 2 − 2x 2

−x

 + ε2 

 ε 2 |x − γ | +ε 1−  − n(γ , ε) dγ 2{(x − γ )2 + ε 2 }3/2 ( x − γ )2 + ε 2 −1/2 with the normalized boundary condition n(1/2, ε) = −1/2. −1

∫



2

2



4

1/2

1

|x − γ |

(102)

At first sight, Eq. (102) looks quite daunting and very different from anything we have examined so far in this text. However, the first terms on the right do indicate the presence of boundary layers near x = ±1/2, so perhaps our familiar techniques can be employed. This concludes all that we will write about this very different type of boundary-layer problem.

570

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

4.4. Semi and superconductors The study of semiconductors and of superconductors, as it has unfolded over the last fifty years or so, has thrown up a large number of interesting and important equations that describe their properties and design characteristics. 4.4.1. Josephson junction The Josephson junction between two superconductors, which are separated by a thin insulator, can produce an AC current when a DC voltage is applied across the junction (by virtue of the tunneling effect). An equation that models an aspect of this phenomenon is u¨ + ε(1 + a cos u) u˙ + sin u = ε b,

t ≥0

(103)

where a and b are given constants and u(0, ε) = u˙ (0, ε) = 0. We will construct the asymptotic solution, using the method of multiple scales, for ε → 0+ . Note that, in the absence of the term ε b, u ≡ 0 is a solution of the complete problem. We anticipate that the presence of ε b will force a non-zero solution which, if remains bounded, should be O(ε) for all time t. 4.4.2. A p–n junction A p–n junction is where two semiconducting materials meet, such junctions may perform different functions. The one that we describe is a diode. We analyze the device for x ∈ [0, 1] where the junction sits at x = 0 (and, by symmetry, it extends into −1 ≤ x < 0) and an ohmic contact is placed at x = 1. In suitable non-dimensional, scaled variables we have

ε ε ε

de dx dp dx dn dx

=p−n+1

(104)

= ep − ε I (x)

(105)

= −en + ε I (x)

(106)

where e is the electrostatic field, p the hole density and n the electron density. The term ‘+1’ in (104) is a constant ‘doping’ density and we will assume that the current density I (x) in (105) and (106) is given. Indeed, in this simple model, we take I (x) = constant. The boundary conditions are p = n at x = 0, n = p + 1 and np = k(>0) both at x = 1 and ε is our parameter with a small value typically about 0.001. It is evident that the set (104)–(106) exhibits the characteristics of a boundary-layer problem because the small parameter multiplies the derivatives in each equation. However, a neat manoeuvre allows one equation to be independent of ε . 4.4.3. Impurities in a semiconductor A significant issue in the design and operation of semiconductors is the presence and movement, of impurities. In particular, the level of impurities that diffuse from the outer surface of the material and move to occupy vacant locations within the structure can be modeled by the equations

  ∂c ∂ ∂c ∂v = v −c , ∂t ∂x ∂x ∂x

∂v ∂c ∂ 2v +α =ε 2, ∂t ∂t ∂x

x > 0, t > 0.

(107)

Here, c (x, t , ε) is the concentration of the impurities, v(x, t , ε) the concentration of vacancies (holes) and α(<1) and ε are positive constants. The boundary and initial are c (x, 0, ε) = 0, v(x, 0, ε) = 1, c (0, t , ε) = 1, v(0, t , ε) = β, c → 0 and v → 1 both as x → ∞ (t ≥ 0). These entire boundary values are constants, β(̸= (1 − α)) is positive and further more, the appearance of the same values at t = 0 and as x → ∞ suggests that we could use a relevant similarity solution. 4.5. Fluid mechanics The study of fluid mechanics is broad and deep and it often has far-reaching consequences. Many of the classical techniques of singular perturbation theory were first developed in order to tackle particular difficulties that were encountered in this field. Examples that are available are numerous and a few of them could be selected for discussion here. 4.5.1. Viscous boundary layer on a flat plate In this example, we consider an incompressible, viscous fluid (in y > 0) flowing over a flat plate, x ≥ 0, the flow direction at infinity being parallel with the plate. The governing equations are the Navier–Stokes equation (in the absence of gravity) and the equation of mass conservation

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575 1 ut + uux + ν uy = −px + R− e (uxx + uyy ),

vt + uvx + νvy = −py + Re (vxx + vyy ), −1

571

(108) ux + v y = 0

(109)

where Re is the Reynolds number (and we have used subscripts throughout to denote Re partial derivatives). We will consider the problem of steady flow (∂/∂ t ≡ 0) with Re → ∞. The boundary conditions for uniform flow at infinity are u → 1 and v → 0 as x → −∞ and as y → ∞, u = 0 and v = 0 on y = 0 in x > 0, which implies that the pressure p → constant away from the plate (and we will not analyze the nature of the flow near x = 0) the plate will extend to infinity (x ≥ 0). 1 The presence of the small parameter R− e multiplying the highest derivatives is the hallmark of a boundary-layer problem. In particular, the (inviscid) problem can satisfy v = 0 on y = 0 but not u → 0 as y → 0 (in x > 0) so we expect a boundary-layer scaling in y. 4.5.2. Very viscous flow past a sphere In this example, we take Re → 0+ (and since Re = Ud/ν , where U is typical speed of the flow, d is typical dimension of the object in the flow and v is the kinematic viscosity, this limit can be interpreted as ‘highly viscous’ or ‘slow flow’ or flow past a ‘small object’). We consider the axisymmetric flow, produced by a uniform flow at infinity parallel with the chosen axis, past a solid sphere. This could be used as a simple model for flow past a rain drop. It is convenient to introduce a stream function Ψ (usually called a Stokes stream function, in this context), which eliminate pressure from the Navier–Stokes equation and hence work with the (non-dimensional) equation

  ∂ ∂ 2 1 4 ψ − ψ + 2 ( cot θ )ψ − ψ D2 ψ = D ψ, θ θ r θ 2 r sin θ ∂r ∂θ r Re 1

1 ≤ r < ∞, 0 ≤ θ ≤ π ,

∂2 1 ∂2 cot θ ∂ + 2 2− 2 2 ∂r r ∂θ r ∂θ ψ 1 with ψ = ψr = 0 on r = 1 and → sin2 θ as r → ∞, 2 where D2 ≡

r 2 this latter condition ensuring that there is an axisymmetric flow of speed one at infinity.

(110)

(111) (112)

4.5.3. A piston problem We consider the one-dimensional flow of a gas in a long, open-ended tube. The gas is brought into motion by the action of a piston at one end, which moves forward at a speed which is much less than the speed of sound in the gas. This is usually called the acoustic problem. The gas is modeled by the isentropic law for a perfect gas Pressure ∝ (Density)γ , 1 < γ < 2 and is described by the equations at + uax +

1 2

(γ − 1)aux = 0,

ut + uux +

2

γ −1

aax = 0,

(113)

where a is the (local) sound speed in the gas and u its speed along the tube. The initial  t and boundary conditions are u = 0 and a = 1 at t = 0, x > ε xp (t ) and u = ε V (t ) on x = ε xp (t ), t ≥ 0, with xp (t ) = 0 V (t ′ ) dt ′ and V (0) = 0. 4.5.4. A variable-depth Korteweg–de Vries equation for water waves We consider the one-dimensional propagation of waves over water (incompressible), which is modeled by an inviscid fluid without surface tension. The water is stationary in the absence of waves, but the local depth varies on the same scale (ε) that is used to measure the weak nonlinearity and dispersive effects in the governing equations. For right-running waves, the appropriate far-field coordinates are ξ = ε −1 χ (X )− t and X = ε x where χ (X ) is to be determined. The non-dimensional equations are

−uξ + ε[u(χ ′ uξ + ε uX ) + wuz ] = −(χ ′ pξ + ε pX ),

(114)

ε[−wξ + ε[u(χ wξ + εwX ) + wwz ]] = pz ,

(115)

′

with p = η and

and

χ uξ + ε uX + w z = 0 ′

w = −ηξ + u(χ nξ + ε nX ) on z = 1 + εη

w = ε uB (X ) on z = B(X ). ′

′

(116) (117)

Here (u, w) are the velocity components of the flow, p is the pressure in the fluid relative to the hydrostatic pressure (with pressure constant at the surface) and z = 1 + εη is the surface of the water. The bottom is represented by the function z = B(X ). 4.6. Extreme thermal processes This next group of problems concerns phenomena that involve explosions, combustion and the other phenomena of this kind.

572

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

4.6.1. A model for combustion A model that aims to describe ignition, followed by a rapid combustion, requires a slow development over a reasonable time scale that precedes a massive change on a very short time scale, initiated by the attainment of some critical condition. A simple (non-dimensional) model for such a process is the equation c˙ = c 2 (1 − c ) with c = ε at t = 0, where c (t , ε) is the concentration of an appropriate chemical that takes part in the combustive reaction. The whole process is initiated by the small disturbance ε(→ 0+ ) at time t = 0. It should be fairly apparent that this equation can be integrated completely to give the solution for c, but in implicit form and so the detailed structure as ε → 0+ is far from transparent. By virtue of the initial value (ε), we first seek a solution in the form c (t , ε) ∼

∞ −

ε n cn (t )

(118)

n =1

and then we obtain c˙1 = 0, c˙2 = c12 , c˙3 = 2c1 c2 − c12 , etc. with c1 (0) = 1 and cn (0) = 0 for n = 2, 3, 4, . . . . 4.6.2. Thermal runaway A phenomenon that can be encountered in certain chemical reactions involves the release of heat (exothermic) which increases the temperature and the temperature normally controls the reaction rate. It is possible, therefore, to initiate a reaction by releasing heat which raises the temperature and hence increases the rate of reaction which releases even more heat and so on, which is called thermal runaway. In the most extreme cases, there is no theoretical limit to the temperature, although physical reality intervenes e.g. the containing vessel might melt or the products explode. A standard model used to describe this is the equation Tt = ∆2 T + λ exp[T /(1 + ε T )]

(119)

where T is the temperature and λ is a parameter. In this example, we will examine the nature of the steady state temperature in one dimension i.e. the solution of Txx + λ exp[T /(1 + ε T )] = 0

as ε → 0+

(120)

for various λ. The boundary condition is that the external temperature is maintained. We will represent this by T = 0 on x = ±1 and hence we seek a solution for x ∈ [−1, 1] and assume symmetry of the temperature distribution about x = 0. 4.7. Chemical and biochemical reactions The examples that are presented here are intended to show that it is possible to model and describe using perturbation theory, some very complex processes that, perhaps fifty years ago, were thought to be mathematically unresolvable. Certainly, some extensive simplification is necessary in the development of the model, and this requires considerable skill and knowledge, but the resulting differential equations remain, generally, quite daunting. 4.7.1. Kinetics of a catalyzed reaction In a model for the kinetics of a particular type of catalyzed reaction, the concentration of the reactant varies according to the equation

ε

dc dx

=1−

2 xc

(1 − c )(γ + c ),

0<x≤1

(121)

with c (1, ε) = 0 and γ (where 0 < γ < 1) as a given constant. The parameter ε (>0) is a rate constant and for ε → 0+ , we observe that the equation reduces to a purely algebraic problem which is readily solved

c=−

1 2



1 2

x+γ −1

 +

1

 4γ +

2



1 2

2

x+γ −1

.

(122)

We have selected the positive sign so that c ≥ 0. For x = 1, this gives the value 1 2



1 2

 −γ +

γ 2 + 3γ +

1

 ̸= 0

4

(123)

and so we will require a boundary layer near x = 1. Note that from this solution, c = 1 at x = 0, so we can expect that c → 1 as x → 0 in the solution of (121). Let us introduce x = 1 − ε X and write c (1 − ε X , ε) ≡ C (X , ε) the equation for C then becomes

− C′ = 1 −

2(1 − C )(γ + C )

(1 − ε X )C

.

(124)

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

573

4.7.2. Enzyme kinetics A standard process in enzyme kinetics concerns the conversion of a substrate (x) into a product, by the action of an enzyme, via a substrate-enzyme complex (y). This is known as the Michaelis–Menton reaction. A model for this process is the pair of equations x˙ = −x + (x + α − β)y,

ε˙y = x − (x + α)y,

(125)

where α and β are positive constants and the initial conditions are x = 1 and y = 0 at t = 0. The parameter ε measures the rate of the production of y and we consider the problem posed above, with ε → 0+ . It is clear that we should expect a boundary layer structure in the solution for y, but not for x. A suitable boundary-layer variable is τ = ε −1 f (t ) rather than simply t /ε and we will take the opportunity to use the method of multiple scales here, using T = t and τ . We introduce x(t , ε) ≡ X (T , τ , ε) and y(t , ε) ≡ Y (T , τ , ε) and then equations (125) become

ε −1 f ′ Xτ + XT = −X + (X + α − β)Y ,

f ′ Yτ + ε YT = X − (X + α)Y

(126)

with X (0, 0, ε) = 1 and Y (0, 0, ε) = 0. 4.7.3. The Belousov–Zhabotinskii reaction This is a famous and much-studied phenomenon, first demonstrated by Belousov in 1951. He discovered that a steady oscillation of the concentration of a catalyst, between its oxidized and reduced states, was possible. In a suitable medium, this can be exhibited as a dramatic change in color, a color being associated with a state, with a period of a minute or so. A set of model equations for this chemical reaction is

ε˙x = x + y − xy − ελx2 ,

εµ˙y = v z − y − xy,

z˙ = x − z ,

(127)

where the concentration of the catalyst is represented by z. The constants λ, µ and v are positive and independent of ε , where measures the rate constant for the production of x, and gives the size of the corresponding constant for y and of the nonlinearity in the first equation of (127). We will consider the problem with ε → 0+ . The initial conditions are relatively unimportant here, but we will assume that they are sufficient to start the oscillatory process. It would be impossible in a text such as ours, with its emphasis on singular perturbation theory, to give a comprehensive description of the relevant solution of the set (127). This would involve, for example, a detailed (but local) stability analysis. We will content ourselves with a brief overview that emphasizes the various scales that are important. We have chosen to ignore many important elements in our presentation of the Belousov–Zhabotinskii- reaction, mainly because they require much that goes well beyond the methods of singular perturbation theory. An illuminating discussion of this process provides an excellent introduction to many chemical, biochemical and biological models and their solutions. 5. Conclusion and further directions for research Finally, it has been pointed out that the ever-increasing number of industrial applications of singular perturbation problems and the richness of behavior of the governing equations make this area a particularly rewarding one for mathematicians, engineers, and industrialists alike. In this paper, we have presented and described basic methods, a number of examples taken from the physical and chemical sciences and in each, the ideas of singular perturbation theory play a significant role. As we implied earlier, such a collection could not be exhaustive; indeed, we can hope to give only an indication of what is possible. Even if the particular applications offered here are of no specific interest to some readers, they do provide a set of additional worked examples that should help to reinforce the ideas that contribute to singular perturbation theory. Here, interested readers are encouraged to investigate the references to related material that have been provided throughout this paper. Acknowledgements The first author acknowledges the financial support provided by the University Grant Commission, New Delhi, Government of India, under research grant no. 37-515/2009(SR). The authors also extend their appreciation to anonymous reviewers for their valuable suggestions to revise this paper. References [1] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman and Hall, CRC Press, London, Boca Raton, FL, 2000. [2] H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, New York, 1996. [3] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [4] D.R. Smith, Singular-Perturbation Theory an Introduction with Applications, Cambridge University Press, Cambridge, 2009. [5] K.W. Morton, Numerical Solution of Convection–Diffusion Problems, in: Applied Mathematics and Mathematical Computation, vol. 12, Chapman and Hall, London, 1996.

574

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

[6] J. Zheng, T.W. Sederberg, R.W. Johnson, Least square methods for solving differential equations using Bezier control points, Applied Numerical Mathematics 48 (2004) 237–252. [7] M. Kumar, P. Singh, H.K. Mishra, An initial-value technique for singularly perturbed boundary value problems via cubic spline, International Journal of Computational Methods in Engineering Science and Mechanics 8 (2007) 1–9. [8] M. Kumar, H.K. Mishra, P. Singh, A boundary value approach for singularly perturbed boundary value problems, Advances in Engineering Software 40 (4) (2009) 298–304. [9] M. Kumar, P. Singh, H.K. Mishra, A recent survey on computational techniques for solving singularly perturbed boundary value problems, International Journal of Computer Mathematics 84 (2007) 1–25. [10] Y.N. Reddy, P.P. Chakravarthy, An initial-value approach for solving singularly perturbed two-point boundary value problems, Applied Mathematics and Computation 155 (2004) 95–110. [11] R.K. Mohanty, U. Arora, A family of non-uniform mesh tension spline methods for singularly perturbed two-point singular boundary value problems with significant first derivatives, Applied Mathematics and Computation 172 (2006) 531–544. [12] M.H. Holmes, Introduction to Perturbation Methods, Springer, 1995. [13] E.J. Hinch, Perturbation Methods, Cambridge University Press, 1991. [14] C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, Asymptotic Methods and Perturbation Theory, SpringerVerlag New York, Inc., 1999. [15] M.R. Owen, M.A. Lewis, How predation can slow, stop, or reverse a prey invasion, Bulletin of Mathematical Biology 63 (2001) 655–684. [16] A.H. Nayfeh, Perturbation Methods, in: Wiley Classics Library, Wiley-Interscience, 2000. [17] http://en.wikipedia.org/wiki. [18] S. Johnson, Singular Perturbation Theory, Mathematical and Analytical Techniques with Applications to Engineering, Springer, 2005. [19] M.J. Ward, Course notes on asymptotic methods, http://www.math.ubc.ca/~ward/teaching/math550.html. [20] R.K. Zeytounian, An asymptotic derivation of the initial condition for the incompressible and viscous external unsteady fluid flow problem, International Journal of Engineering Science 38 (2000) 1983–1992. [21] H.E. Williams, An asymptotic solution of the governing equation for the natural frequencies of a cantilevered, coupled beam model, Journal of Sound and Vibration 312 (2008) 354–359. [22] G. Iosilevskii, Asymptotic theory of an oscillating wing section in weak ground effect, European Journal of Mechanics, B/Fluids 27 (2008) 477–490. [23] A. Carpinteri, M. Paggi, Asymptotic analysis in linear elasticity: from the pioneering studies by Wieghardt and Irwin until today, Engineering Fracture Mechanics 76 (2009) 1771–1784. [24] F. Lebon, R. Rizzoni, Asymptotic analysis of a thin interface, the case involving similar rigidity, International Journal of Engineering Science 48 (2010) 473–486. [25] M.P. Fernando, A.A.M. Oliveira, F.F. Fernando, Asymptotic analysis of stationary adiabatic premixed flames in porous inert media, Combustion and Flame 156 (2009) 152–165. [26] L.C. Nitsche, W. Zhang, L.E. Wedgewood, Asymptotic basis of the L closure for finitely extensible dumbbells in suddenly started uniaxial extension, Journal of Non-Newtonian Fluid Mechanics 133 (2006) 14–27. [27] I.V. Andrianov, Vladyslav V. Danishevs’kyy, D. Weichert, Asymptotic determination of effective elastic properties of composite materials with fibrous square shaped inclusions, European Journal of Mechanics, A/Solids 21 (2002) 1019–1036. [28] D.J. Steigmann, Asymptotic finite-strain thin-plate theory for elastic solids, Computers and Mathematics with Applications 53 (2007) 287–295. [29] C. Pelissou, F. Lebon, Asymptotic modeling of quasi-brittle interfaces, Computers and Structures 87 (2009) 1216–1223. [30] S. Saintlos, J. Mauss, Asymptotic modelling for separating boundary layers in a channel, International Journal of Engineering Science 34 (2) (1996) 201–211. [31] A.P. Bunger, E. Detournay, Asymptotic solution for a penny-shaped near-surface hydraulic fracture, Engineering Fracture Mechanics 72 (2005) 2468–2486. [32] W. Kelley, Asymptotically singular boundary value problems, Mathematical and Computer Modelling 32 (2000) 541–548. [33] K.M. Yip, Model simplification by asymptotic order of magnitude reasoning, Artificial Intelligence 80 (1996) 309–348. [34] N. Popovic, P. Szmolyan, Rigorous asymptotic expansions for Lagerstrom’s model equation—a geometric approach, Nonlinear Analysis 59 (2004) 531–565. [35] F. Verhulst, Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics, Springer, 2005. [36] J. Priede, G. Gerbeth, Matched asymptotic solution for the solute boundary layer in a converging axisymmetric stagnation point flow, International Journal of Heat and Mass Transfer 50 (2007) 216–225. [37] O. Radulescu, P.D. Olmsted, Matched asymptotic solutions for the steady banded flow of the diffusive Johnson–Segalman model in various geometries, Journal of Non-Newtonian Fluid Mechanics 91 (2000) 143–164. [38] J. Margerit, O. Sero-Guillaume, Study of the evaporation of a droplet in its stagnant vapor by asymptotic matching, International Journal of Heat and Mass Transfer 39 (18) (1996) 3887–3898. [39] K. Gersten, The method of matched asymptotic expansions, Von Karman Institute for Fluid Mechanics, Mathematical Methods in Fluid Mechanics, 1980. [40] J. Mauss, J. Cousteix, Uniformly valid approximation for singular perturbation problems and matching principle, Comptes Rendus (CR) Mecanique 330 (2002) 697–702. [41] J. Kevorkian, J.D. Cole, Multiple Scale and Singular Perturbation Methods, Springer, 1996. [42] R. Khanin, M. Cartmella, A. Gilbert, A computerized implementation of the multiple scales perturbation method using Mathematica, Computers and Structures 76 (2000) 565–575. [43] D. Dessi, F. Mastroddi, L. Morino, A fifth-order multiple-scale solution for Hopf bifurcations, Computers and Structures 82 (2004) 2723–2731. [44] A.A.P. Zanoosi, A. Shooshtari, A multiple times scale solution for non-linear vibration of mass grounded system, Applied Mathematical Modelling 34 (2010) 1918–1929. [45] S. Zhu, P. Yu, Computation of the normal forms for general M-DOF systems using multiple time scales. Part I: autonomous systems, Communications in Nonlinear Science and Numerical Simulation 10 (2005) 869–905. [46] O. Coulaud, Multiple time scales and perturbation methods for high frequency electromagnetic–hydrodynamic coupling in the treatment of liquid metals, Nonlinear Analysis: Theory, Methods & Applications 30 (6) (1997) 3637–3643. [47] C.D. Munz, R. Fortenbach, M. Dumbser, Computational Aero Acoustics: from acoustic sources modeling to farfield radiated noise prediction multiple scale modelling of acoustic sources in low Mach-number flow, Comptes Rendus (CR) Mecanique 333 (2005) 706–712. [48] W.T. Van Horssen, On integrating vectors and multiple scales for singularly perturbed ordinary differential equations, Nonlinear Analysis 46 (2001) 19–43. [49] F. Lakrad, M. Belhaq, Periodic solutions of strongly non-linear oscillators by the multiple scales method, Journal of Sound and Vibration 258 (4) (2002) 677–700. [50] M. Belhaq, F. Lakrad, The elliptic multiple scales method for a class of autonomous strongly non-linear oscillators, Journal of Sound and Vibration 234 (3) (2000) 547–553. [51] Paul Filippi, Acoustics: Basic Physics, Theory and Methods, Academic Press, 1999, 171. [52] L. Gosse, N.J. Mauser, Multiphase semiclassical approximation of an electron in a one dimensional crystalline lattice-III. From ab initio models to WKB for Schrodinger–Poisson, Journal of Computational Physics 211 (2006) 326–346. [53] A.M. Wazwaz, A comparative study on a singular perturbation problem with two singular boundary points, Applied Mathematics and Computation 99 (1999) 179–193.

M. Kumar, Parul / Mathematical and Computer Modelling 54 (2011) 556–575

575

[54] N.B. Abdallah, O. Pinaud, Multiscale simulation of transport in an open quantum system-resonances and WKB interpolation, Journal of Computational Physics 213 (2006) 288–310. [55] A.E. Gill, Atmosphere Ocean Dynamics, in: International Geophysics Series, vol. 30, Academic Press, 1982. [56] R. Spigler, M. Vianello, A survey on the Liouville–Green (WKB) approximation for linear difference equations of the second order, in: Saber Elaydi, I. Gyori, G.E. Ladas (Eds.), Advances in Difference Equations: Proceedings of the Second International Conference on Difference Equations: Veszprém, Hungary, CRC Press, 1995, p. 567. [57] R. Krivec, V.B. Mandelzweig, Quasilinearization method and WKB, Computer Physics Communications 174 (2006) 119–126. [58] P.G. Drazin, Nonlinear Systems, Cambridge University Press, 1992, pp. 181–186. [59] H. Poincare, Les Méthodes Nouvelles de la Mécanique Célèste, vol. II, Dover Publication, New York, 1957, 123–128. [60] A. Lindstedt, Beitrag zur Integration der Differentialgleichungen der Storungstheorie, Abh. K. Akad. Wiss. St. Petersburg 31 (1882). [61] J.D. Logan, Applied Mathematics, second ed., John Wiley and Sons, 1997. [62] C.H. Yanga, S.M. Zhu, S.H. Chen, A modified elliptic Lindstedt–Poincare method for certain strongly non-linear oscillators, Journal of Sound and Vibration 273 (2004) 921–932. [63] H.M. Liu, Approximate period of nonlinear oscillators with discontinuities by modified Lindstedt–Poincare method, Chaos, Solitons & Fractals 23 (2005) 577–579. [64] H. Hu, Z.G. Xiong, Comparison of two Lindstedt–Poincare type perturbation methods, Journal of Sound and Vibration 278 (2004) 437–444. [65] O. Khrustalev, S. Vernov, Construction of doubly periodic solutions via the Poincare–Lindstedt method in the case of massless ϕ 4 theory, Mathematics and Computers in Simulation 57 (2001) 239–252. [66] T. Ozis, A. Yildirim, Determination of periodic solution for a u1/3 force by He’s modified Lindstedt–Poincare method, Journal of Sound and Vibration 301 (2007) 415–419. [67] R.R. Pusenjak, Extended Lindstedt–Poincare method for non-stationary resonances of dynamical systems with cubic nonlinearities, Journal of Sound and Vibration 314 (2008) 194–216. [68] P. Amore, A. Aranda, Improved Lindstedt–Poincare method for the solution of nonlinear problems, Journal of Sound and Vibration 283 (2005) 1115–1136. [69] S. Abbasbandy, Iterated He’s homotopy perturbation method for quadratic Riccati differential equation, Applied Mathematics and Computation 175 (2006) 581–589. [70] S. Benbachir, Lindstedt–Poincare method and periodic families of the Barbanis contopoulos Hamiltonian system, Mathematics and Computers in Simulation 51 (2000) 579–596. [71] A. Marasco, Lindstedt–Poincare method and mathematica applied to the motion of a solid with a fixed point, Computers and Mathematics with Applications 40 (2000) 333–343. [72] G.M. Abd, E.I. Latif, On a problem of modified Lindstedt–Poincare for certain strongly non-linear oscillators, Applied Mathematics and Computation 152 (2004) 821–836. [73] J.H. He, Modified Lindstedt Poincare methods for some strongly non-linear oscillations part II: a new transformation, International Journal of NonLinear Mechanics 37 (2002) 315–320. [74] S.H. Chen, J.L. Huang, K.Y. Szeb, Multidimensional Lindstedt–Poincare method for nonlinear vibration of axially moving beams, Journal of Sound and Vibration 306 (2007) 1–11. [75] J.F. Navarro, On the implementation of the Poincare–Lindstedt technique, Applied Mathematics and Computation 195 (2008) 183–189. [76] S.H. Chen, X.M. Yang, Y.K. Cheung, Periodic solutions of strongly quadratic nonlinear oscillators by the elliptic Lindstedt Poincare method, Journal of Sound and Vibration 227 (5) (1999) 1109–1118. [77] J.H. He, The homotopy perturbation method for nonlinear oscillators with discontinuities, Applied Mathematics and Computation 151 (2004) 287–292. [78] H. Yuce, C.Y. Wang, A boundary perturbation method for circularly periodic plates with a core, Journal of Sound and Vibration 328 (2009) 345–368. [79] J.F. Couchouron, M. Kamenski, R. Precup, A nonlinear periodic averaging principle, Nonlinear Analysis 54 (2003) 1439–1467. [80] M. Kamenskii, P. Nistri, An averaging method for singularly perturbed systems of semi linear differential inclusions with analytic semigroups, Nonlinear Analysis 53 (2003) 467–480. [81] L.D. Akulenko, Higher-order averaging schemes in systems with fast and slow phases, Journal of Applied Mathematics and Mechanics5 66 (2) (2002) 153–163. [82] J. Gine, J. Llibre, Limit cycles of cubic polynomial vector fields via the averaging theory, Nonlinear Analysis 66 (2007) 1707–1721. [83] G.M. Mahmoud, S.A. Aly, Periodic attractors of complex damped non-linear systems, International Journal of Non-Linear Mechanics 35 (2000) 309–323. [84] V. Marinca, N. Herisanu, Periodic solutions for some strongly nonlinear oscillations by He’s variational iteration method, Computers and Mathematics with Applications 54 (2007) 1188–1196. [85] G. Lin, Periodic solutions for van der Pol equation with time delay, Applied Mathematics and Computation 187 (2007) 1187–1198. [86] F. Lakrad, M. Belhaq, Periodic solutions of strongly non-linear oscillators by the multiple scales method, Journal of Sound and Vibration 258 (4) (2002) 677–700. [87] Z.H. Wang, H.Y. Hub, H.L. Wang, Robust stabilization to non-linear delayed systems via delayed state feedback—the averaging method, Journal of Sound and Vibration 279 (2005) 937–953. [88] F. Dohnal, Suppressing self-excited vibrations by synchronous and time-periodic stiffness and damping variation, Journal of Sound and Vibration 306 (2007) 136–152. [89] T. Ozis, A. Yildirim, Generating the periodic solutions for forcing van der Pol oscillators by the Iteration Perturbation method, Nonlinear Analysis. Real World Applications 10 (2009) 1984–1989. [90] J. Sanders, F. Verhulst, J. Murdock, Averaging Methods in Nonlinear Dynamical Systems, Springer, New York, 2007.