Perturbation Expansion to the Solution of Differential Equations

The main purpose of this chapter is to describe the application of perturbation expansion techniques to the solution of differential equations. Approximate expressions are generated in the form of asymptotic series. These may not and often do not converge but in a truncated form of only two or three terms, provide a useful approximation to the original problem. These analytical techniques provide an alternative to the direct computer solution. Before attempting to solve these problems numerically, one should have an awareness of the perturbation approach.


Introduction
The governing equations of physical, biological and economical models often involve features which make it impossible to obtain their exact solution. For instance, problems where we observe "a complicated algebraic equations", "the occurrence of a complicated integral", in case of differential equations (DE), "a varying coefficients or nonlinear term" sometimes problems with an awkwardly shaped boundary are tough to solve with the limited methods for finding analytical solutions. The main purpose of this chapter is to describe the application of perturbation expansion techniques to the solution of DE. Approximate expressions are generated in the form of asymptotic series. These may not and often do not converge but in a truncated form of only two or three terms, provide a useful approximation to the original problem. These analytical techniques provide an alternative to the direct computer solution. Before attempting to solve these DE numerically, one should have an awareness of the perturbation approach. An example of this occurs in boundary layer problems where there are regions of rapid change of quantities such as fluid velocity, temperature or concentration. Appropriate scaling of the boundary layer dimension is required before a numerical solution can be generated which will capture the behavior in the rapidly changing region.
When a large or small parameter occurs in a mathematical model of a process there are various methods of constructing perturbation expansions for the solution of the governing equations. Often the terms in the perturbation expansions are governed by simpler equations for which the exact solution techniques are available. Even if exact solutions cannot be obtained, the numerical methods used to solve the perturbation equations approximately are often easier to construct than the numerical approximation for the original governing equation.
First, we consider a model problem for which an exact solution is available against which the perturbation expansion can be compared. A feature of the perturbation expansions is that they often form divergence series. The concept of an asymptotic expansion will be introduced and the value of a truncated divergent series will be demonstrated.

Projectile motion
This example studies the effect of small damping on the motion of a particle. Consider a particle of mass M which is projected vertically upward with an initial speed U 0 . Let U denote the speed at some general time T. If air resistance is neglected then the only force acting on the particle is gravity, ÀMg, where g is the acceleration due to gravity and the minus sign occurs because the upward direction is chosen to be the positive direction. Newton's second law governs the motion of the projectile, i.e., Integrating (1), we obtain the solution U ¼ C À gT. The constant of integration is determined from the initial condition U 0 ð Þ ¼ U 0 , so that On defining the non-dimensional velocity v, and time t, by v ¼ U=U 0 and t ¼ gT=U 0 , the governing equation becomes with the solution v t ð Þ ¼ 1 À t. Taking account of the air resistance, and is included in the Newton's second law as a force dependent on the velocity in a linear way, we obtain the following linear equation where the drag constant K is the dimensions of masa/time. In the nondimensional variables, it becomes Let us denote the dimensionless drag constant by ε, then the governing equation is where ε > 0 is a "small" parameter and the disturbances are very "small". The damping constant K in (4) is small, since K has the dimensions of mass/time and a small quantity in units of kilograms per second.

Perturbation expansion
It is possible to solve (6) exactly since it is of variables separable form. Here, we solve by an iterative process, known as perturbation expansion for the solution.
Let v i ð Þ denotes the ith iterate, which is obtained from the equation The justification for this iterative scheme is that the term εv involves the small multiplying coefficient ε, and so the term itself may be expected to be small. Thus, the term εv i ð Þ which should appear on the RHS of (7) to make it exact, may be replaced by εv iÀ1 ð Þ with an error which is expected to be small. The first iterate is obtained by neglecting the perturbation, thus This is known as the unperturbed problem, and direct integration yields v 0 ð Þ ¼ 1 À t: Rearranging the terms in these iterates in ascending powers of ε, we obtain Clearly as the iteration proceeds the expressions are refined by terms which involve increasing powers of ε. These terms become progressively smaller since ε is a small parameter. This is an example of a perturbation expansion. It will often be the case that perturbation expansions involve ascending integer powers of the small parameter, i.e., ε 0 , ε 1 , ε 2 , ⋯ È É . Such a sequence is called an asymptotic sequence. Although this is the most common sequence which we shall meet, it is by no means unique. Examples of other asymptotic sequences are ε 1=2 , ε, ε 3=2 , ε 2 , ⋯ È É and ε 0 , ε 2 , ε 4 , ⋯ È É . In each case the essential feature is that subsequent terms tend to zero faster than previous terms as ε ! 0.
An alternative procedure to that of developing the expansion by iteration is to assume the form of the expansion at the outset. Thus, if we assume that the perturbation expansion involves the standard asymptotic sequence ε 0 , ε 1 , ε 2 , ⋯ È É , then the solution v, which depends on the variable t, and the parameter ε, is expressed in the form The coefficients v 0 t ð Þ, v 1 t ð Þ, ⋯ of powers of ε are functions of t only. Substituting expansion (9) in the governing Eq. (6) yields the following Thus, the coefficients of powers of ε can be equated on the left-and right-hand sides of (10): The proof of validity of this fundamental procedure can be developed by first setting ε ¼ 0 in (10) which yields the first equation of (11). This result allows the first member of the left-and right-hand side of Eq. (10) to be removed. Then, after dividing the remaining terms by ε we obtain the equation This is valid for all nonzero values of ε so that on taking the limit as ε ! 0 we obtain the second equation of (11). Repeating the procedure, we obtain the other equations.
Integrating the equations in (11), we obtain Using these values in (9), we obtain that v t; This is the same as the expansion (8) which is generated by iteration.
The IVP (6) can be solved exactly as The perturbation expansion can be obtained from (12) by replacing the exponential function by its Maclaurin expansion, i.e., v t ð Þ ¼ 1 This is the same as the expansion (12). Thus, the perturbation expansion approach is justified in this case. One can refer the books [1,2].

Asymptotics
The letters O and o are order symbols. They are used to describe the rate at which functions approach limit values. We will consider the types of limit values, namely zero, a finite number but nonzero and infinite.
If a function f x ð Þ approaches a limiting value at the same rate of another function g x ð Þ as x ! x 0 , then we write The functions are said to be of the same order as x ! x 0 . The test for this is the limit of the ratio. Thus, if lim where C is finite, then we say (15) holds. For example, we have the following functions: The expression This is a stronger assertion that the corresponding Oformula. The relation (16) implies the relation (15), as convergence implies boundedness from a certain point onwards.
We have the following functions satisfy the o-relation:

Asymptotic expansions
Consider the expansion is an asymptotic expansion as x ! ∞, if, for any N, The following expansion is used when (17) and (18) Here, lim n!∞ R N ¼ 0, for any value of N. The sequence 1, 1=x, 1=x 2 , ⋯ È É is an asymptotic sequence as x ! ∞. The characteristic feature of such sequences is that each member is dominated by the previous member. In constructing examples it is easier to deal with the limit zero than any other. Thus, for the case x ! ∞, we let ε ¼ 1=x, which for x ! x 0 , we let ε ¼ x À x 0 so that without loss of generality we may confirm our attention to the limit ε ! 0. The standard asymptotic sequence is 1, ε, ε 2 , ⋯ È É as ε ! 0. If we let δ n ε ð Þ represent members of an asymptotic sequence δ 0 ε ð Þ, δ 1 ε ð Þ, ⋯ f g as ε ! 0, then the following condition must hold Some examples of asymptotic sequences are The general expression for an asymptotic expansion of a function f ε ð Þ, in terms of an asymptotic sequence δ n ε ð Þ is where the coefficients a n are independent of ε. The expression (20) involving the symbol $, means that for all N, where If a function possesses an asymptotic expansion involving the sequence δ 0 ε ð Þ, δ 1 ε ð Þ, ⋯ f g then the coefficients a n of the expansion (21) given by the expression (24) are unique. However, another function may share the same set of coefficients. Thus, while functions have unique expansions, an expansion does not correspond to a unique function.
Consider a function f x; ε ð Þ, which depends on both an independent variable x, and a small parameter ε. Suppose that f x; ε ð Þ is expanded using an asymptotic sequence δ n ε ð Þ f g, The coefficients of the gauge functions δ n ε ð Þ are functions of x, and the remainder after N terms is a function of both x and ε. For this to be an asymptotic expansion, we require Refer [3,4] for more details. For (24) to be a uniform asymptotic expansion the ultimate proportionality between R N and δ Nþ1 must be bounded by a number independent of x, i.e., for ε in the neighborhood near zero, where K is a fixed constant. An example of a uniform asymptotic expansion is f Here, one cannot find a fixed K which satisfy |R N | ≤ K|ε Nþ1 |, because for any choice of K, x can be chosen so that x Nþ1 exceeds this value.

Nonuniformity
The expansion (27) becomes nonuniform when subsequent terms are no longer small corrections to previous terms. This occurs when subsequent terms are of the same order or of dominant order than previous terms. Subsequent terms dominate previous terms for larger x, for example, when ð Þ since then subsequent terms decrease by a factor of ε. The expansion remains valid for large x, provided x is not as large as 1=ε. For instance, the expansion is valid for The critical case is such that subsequent terms are of the same order. This determines the region of nonuniformity. In (27), the region of nonuniformity occurs when εx

Sources of nonuniformity
There are two common reasons for nonuniformities in asymptotic expansions, they are 1. Infinite domains which allow long-term effects of small perturbations to accumulate.
2. Singularities in governing equations which lead to localized regions of rapid change.
Consider the nonlinear Duffing equation Suppose the solution may be expanded using the standard asymptotic sequence On substituting this in (28) and in the initial conditions, we get Equating like of powers of ε on both sides, we get and O ε ð Þ : Solving Eqs. (30) and (31), we obtain The term t sin t ð Þ in the expansion (32) is called a secular term. It is an oscillating term of growing amplitude. All other terms are oscillating of fixed amplitude. The secular term leads to a nonuniformity for large t. The region of nonuniformity is obtained by equating the order of the first and second terms, The trigonometric functions are treated as O 1 ð Þ terms. Thus, the region of nonuniformity is t ¼ O 1=ε ð Þ, as ε ! 0. The second common source of nonuniformities is associated with the presence of singularities. Consider, the following initial-value problem: where ε > 0 is a small parameter. Suppose y has the expansion Substituting (34) in (33), we have Equating coefficients of like powers of ε on both sides, we get Clearly, y 0 cannot satisfy the boundary condition y 0 0 ð Þ ¼ 2 as no constant of integration is available because the equation determining y 0 is an algebraic equation not a differential equation, and no additional conditions are required. Thus, we have obtain the expression but the initial condition y 0 ð Þ ¼ 2 has not been satisfied. The unperturbed problem, obtained by setting ε ¼ 0 is not a DE, but an algebraic equation y ¼ e Àx . This cannot satisfy an arbitrarily imposed condition at x ¼ 0. For any nonzero value of ε, (33) becomes a first-order DE which can satisfy an initial condition. This is an example of a singular perturbation problem (SPP), where the behavior of the perturbed problem is very different from that of the unperturbed problem.
Thus, the perturbation expansion (36) is a good approximation of the exact solution away from the region x ¼ 0. To see this, let us compare (36) with the following exact solution: The perturbation expansion (36) generates the second member of (37), but not the first member. The coefficient e Àx=ε Þ is a rapidly varying function which takes the value of unity at x ¼ 0, and rapidly decays to zero for x > 0. Clearly, y 0 provides a good approximation away from the region x ¼ 0. The region near x ¼ 0 is called the boundary layer. These regions usually occur when the highest order derivative of a DE is multiplied by a small parameter. The unperturbed problem, obtained by setting ε ¼ 0 is of lower order and consequently cannot satisfy all the boundary conditions. This leads to boundary layer regions where the solution varies rapidly in order to satisfy the boundary condition.
Boundary layers are regions of nonuniformity in perturbation expansions of the form (36).

Boundary layer
Boundary layers are regions in which a rapid change occurs in the value of a variable. Some physical examples include "the fluid velocity near a solid wall", "the velocity at the edge of a jet of fluid", "the temperature of a fluid near a solid wall." Ludwig Prandtl pioneered the subject of boundary layer theory in his explanation of how a quantity as small as the viscosity of common fluids such as water and air could nevertheless play a crucial role in determining their flow. The viscosity of many fluids is very small and yet taking account of this small quantity is vital. The essential point is that the viscous term involves higher order derivatives so that its omission necessitates the loss of a boundary condition. The ideal flow solution allow slip to occur between a solid and fluid. In reality the tangential velocity of a fluid relative to a solid is zero. The fluid is brought to rest by the action of a tangential stress resulting from the viscous force.
Mathematically the occurrence of boundary layers is associated with the presence of a small parameter multiplying the highest derivative in the governing equation of a process. A straightforward perturbation expansion using an asymptotic sequence in the small parameter leads to differential equations of lower order than the original governing equation. In consequence not all of the boundary and initial conditions can be satisfied by the perturbation expansion. This is an example of what is commonly referred to as a singular perturbation problem. The technique for overcoming the difficulty is to combine the straightforward expansion, which is valid away from the layer adjacent to the boundary. The straightforward expansion is referred to as the outer expansion. The inner expansion associated with the boundary layer region is expressed in terms of a stretched variable, rather than the original independent variable, which takes due account of the scale of certain derivative terms. The inner and outer expansions are matched over a region located at the edge of the boundary layer. The technique is called the method of matched asymptotic expansions.
Consider the following two-point boundary value problem: where ε > 0 is a small parameter. If we assume that u possesses a straightforward expansion in powers of ε, then the equations associated with powers of ε leads to O ε n ð Þ : and the boundary conditions require which leads to Equation (42) require that each u n x ð Þ satisfy two boundary conditions. This is in general impossible since Eqs. (41) and (42) governing each u n are of first-order. Now the question is which one of the boundary condition has to be taken into account. We will find out that the boundary condition at x ¼ 0 must be abandoned and consequently the expansion (39) is invalid near x ¼ 0.
The general solution of (42) is u 0 x ð Þ ¼ x 2 þ x þ C, using the boundary condition u 0 1 ð Þ ¼ 4, we obtain From (42), we obtain the equations and its solutions are Therefore, the outer expansion is where 'out' label is used to indicate that the solution is valid away from the region near x ¼ 0. Clearly u out fails to satisfy the boundary condition at x ¼ 0. The reason why the outer solution is of use is that it closely follows the exact solution of the problem except in a narrow region near x ¼ 0, where the exact solution changes rapidly in order to satisfy the boundary condition.
The exact solution of the BVP (38) can be obtained as The constants A and B are determined from the boundary conditions: We know that e À1=ε ¼ o ε N À Á , as ε ! 0, for all N. This means that the exponential term tends to zero faster than any power of ε, as ε ! 0. It is called a transcendentally small term (T.S.T.) and can always be neglected since its contribution is asymptotically always less than any power of ε. Thus, (45) gives and the exact solution is after rearranging the terms in asymptotic order, we obtain Comparing the exact solution with the outer expansion shows that the terms involving e Àx=ε are absent. The effect of these terms is negligible when It is clear that as ε ! 0 the region in which the outer solution departs from the exact solution becomes arbitrarily close to x ¼ 0 with a thickness O ε ð Þ. This region is called the boundary layer. The behavior of the exact solution and the zeroth-order term of the outer expansion are plotted in Figure 1 for various values of ε. By differentiating the leading order term u ex 0 , of the exact solution, we have Outside the boundary layer, i.e., for , ∀N, so ε À1 e Àx=ε and ε À2 e Àx=ε are also transcendentally small. Within the boundary layer The order of u ex 0 and its derivatives are given below: This indicates that x is the appropriate independent variable outside the boundary layer where u ex 0 and its derivatives are of O 1 ð Þ quantities. However, within the boundary layer the appropriately scaled independent variable is s ¼ x=ε, then so that within the boundary layer The variable s ¼ x=ε is called a stretched variable. The differential equations becomes We assume a boundary layer expansion, called the inner expansion of the form The inner expansion will satisfy the boundary condition at giving v 0 0 ð Þ ¼ 1, and v n 0 ð Þ ¼ 0, n ¼ 1, 2, ⋯. Substituting (49) into the DE (48), we obtain the following set of equations: The boundary condition at x ¼ 1 cannot be used to determine the constants appearing in these solutions because the DEs (50) are only valid in the boundary layer. The constants in (51) are determined by matching the inner and outer expansions. We shall first restrict our attention to matching the leading order expansions u 0 and v 0 . The method which we shall apply is Prandtl's matching condition.
The leading order terms in the 'inner' and 'outer' expansions are to be matched at the 'edge of the boundary layer'. Of course there is no precise edge of the boundary layer, we simply know that it has thickness of order O ε ð Þ. A plausible matching procedure would be to equate u 0 and v 0 at a value of x such that the region of rapid change has passed. We might choose to equate the terms at the point x ¼ 5ε. The leading order expansions are Equating at x ¼ 5ε gives the following: If, instead we choose to match at x ¼ 6ε, then we obtain These two expressions differ in the argument of the exponential and differ algebraically with 5ε replaced by 6ε. The exponential functions are approaching transcendentally small values so that their contribution can be neglected. The algebraic difference is of O ε ð Þ. Thus, the arbitrariness in the decision of the point at which we choose to equate the expansions leads to a difference of O ε ð Þ. But we are only dealing with leading order expansions anyway. The difference between the exact solution and the leading order expansions will of O ε ð Þ so that an arbitrariness in v 0 and u 0 of O ε ð Þ is immaterial. Rather than choose between, for example, 5ε and 6ε as the value of x to evaluate u 0 we may take the value at x ¼ 0, since where the remainder is uniformly O ε ð Þ since the gradient of u 0 is O 1 ð Þ. For the inner expansion we are to ensure that the rapidly varying function has achieved its asymptotic value at the edge of the boundary layer. This means that the term e Àx=ε should be replaced by zero. This can be achieved by taking the limit s ! ∞. Thus, rather than choosing a specific point to equate the inner and outer terms er are led to the following Prandtl's matching condition: The limit s ! ∞ may appear rather dangerous since although it certainly removes the exponential term it could lead to an algebraically unbounded term. For example, if v 0 ¼ As þ 1 À A ð Þe Às , then the first member would be unbounded as s ! ∞. This possibility can be eliminated since the inner expansion must be of a form which varies rapidly for , not for s ! ∞. In practice, if the boundary layer has been properly located and the correct inner variable is used then Prandtl's matching condition is valid and elegantly avoids the need to choose an arbitrary 'edge' of the boundary layer.
Applying these conditions to the current example leads to lim x!0 which yields A ¼ 2. Thus the leading order terms in the expansion solutions are To prove that these are valid leading terms we consider u ex : We conclude that the matching condition has correctly predicted the leading order terms.

Composite expansion
As single composite expression for these leading order terms can be constructed using the combination where u match 0 is given by (52). Then, For the current example, u match 0 ¼ 2, so the composite expansion is Prandtl's matching condition can only be used for the leading order terms in the asymptotic expansions.
The outer, inner and composite expansions of the BVP (38) are presented in Figures 2 and 3 for different values of ε. From these figures, one can easily identify the need and efficiency of the composite expansion.

Boundary layer location
Consider the following linear DE The following general statements can be made about the boundary layer location and the nature of the inner expansion.
Case I. If a x ð Þ > 0 throughout x 1 , x 2 ð Þ, then the boundary layer will occur at x ¼ x 1 . The stretching transformation will be s ¼ x À x 1 ð Þ=ε, and the one-term inner expansion will satisfy The solution of this equation is The other condition to determine the constants A and B is obtained by matching with the value of the outer expansion at x ¼ x 1 .
Case II. If a x ð Þ < 0 throughout x 1 , x 2 ð Þ, then the boundary layer will occur at x ¼ x 2 . The stretching transformation will be s ¼ x 2 À x ð Þ=ε, and the one-term inner expansion will involve the rapidly decaying function e a x 2 ð Þ x 2 Àx ð Þ =ε . Case III. If a x ð Þ changes sign in the interval x 1 < x < x 2 , then a boundary layer occurs at an interior point x 0 , where a x 0 ð Þ ¼ 0 and boundary layers may also occur at both ends x 1 and x 2 .

Boundary layer thickness and the principle of least degeneracy
The boundary layers which we have met so far have all had thickness O ε ð Þ. By this we mean that a variation of O ε ð Þ in the independent variable will encompass the region of rapid change in the dependent variable. The associated stretched independent variable s, appropriate for the boundary layer is related to x by a linear transformation involving division by ε.
There are practical situations where the boundary layer thickness will be of O ε p ð Þ. This means that if the boundary layer is located at x ¼ x 0 , then the appropriate stretching transformation is s ¼ x À x 0 ð Þ=ε p . More generally, the choice of the function δ ε ð Þ to use in the stretching transformation s ¼ x À x 0 ð Þ=δ ε ð Þ is determined by the need to represent the region of rapid change correctly. We must ensure that the boundary layer solution contains rapidly varying functions. The form of the governing equation in the boundary layer region must have sufficient structure to allow such solutions.
Consider the example Since the signs of the first and second derivatives are the same, and the boundary layer will occur at x ¼ 0. We are not going to assume at the outset that the boundary layer thickness is O ε ð Þ. Our intension is to deduce that the appropriate stretching variable is s ¼ x=ε.
The one-term outer expansion u 0 satisfies du To determine the inner expansion we first wrongly assume that the boundary layer thickness is O ε 1=2 À Á . The stretching transformation s ¼ x=ε 1=2 changes the original DE (56) into the following one: The one-term composite expansion is The leading order boundary layer equation associated with the stretching transformation s ¼ x=ε, (61) involves more terms than (59), associated with s ¼ x=ε 1=2 , and (60) associated with s ¼ x=ε 2 . The extra term in (61) allows sufficient structure in the solution to produce the required boundary layer behavior. An aid for choosing the boundary layer thickness is to seek a stretching transformation which retains the largest number of terms in the dominant equation governing v 0 . This referred to as the principle of least degeneracy by Van Dyke.
The composite expansion (62) can be verified by comparing with the exact solution of (56). The general solution of (56) is , m 2 ¼ À1 À ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À 4ε p 2ε : We expand ffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À 4ε p using the binomial series, so that u ex ¼ C 1 e Àx þ C 2 e Àx=ε Á e x þ x À 1 ð ÞþO ε ð Þ: Using the boundary conditions and by neglecting the transcendentally small term e À1=ε , we have C 1 ¼ 2e, C 2 ¼ 2 1 À e ð Þ. Then, (63) becomes u ex ¼ 2e 1Àx þ 2 1 À e ð Þe Àx=ε Á e x þ x À 1 ð ÞþO ε ð Þ: There is an apparent discrepancy between (64) and the composite expansion (62) in the coefficient of the e Àx=ε term. There is an extra term only contributes in the boundary layer where x ¼ O ε ð Þ so that the coefficient e x may to leading order, be replaced by unity. Thus, the leading order composite expansion and the leading order term in the exact solution are in complete agreement.

Boundary layer of thickness of O
ffiffi ffi ε p ð Þ Consider the following two-point BVP: We seek a one-term composite expansion for the above BVP. We will tentatively assume that a boundary layer occurs at x ¼ 0 although the vanishing of the coefficient of the first derivative suggests the possibility of nonstandard behavior.
While the corresponding coefficient is negative in the range À1 < x < 0 indicates a boundary layer located at the right-hand limit which is again is x ¼ 0. Thus, we are led to expr = ect two outer expansions for positive and negative x respectively and an inner expansion in the interior layer located at x ¼ 0. We denote the leading term in the outer expansion for positive x by u þ 0 , it satisfies with the solution u þ 0 x ð Þ ¼ 2e Àx . The outer expansion for negative x, u À 0 satisfies with the solution u À 0 x ð Þ ¼ e Àx . We suppose the boundary layer at x ¼ 0 has thickness O ε p ð Þ and determine the index p using the principle of least degeneracy. Let s ¼ x=ε p so that the DE becomes ε 1À2p d 2 v ds 2 þ s dv ds þ ε p sv ¼ 0: 2u 0