A singular perturbation method Part I and II

lab.

v. Scheepslouwkuncfe

ARCHIEF

.Technische Hogeschool

DellI

Division of Engineering and Applied Physics

-,Rngineering Sctences Laboratory

- Harvard University

-/

A SINGULAR PERTURBATION MT}IOD. PART I

By

N.- D. PoWkes

1967

A SINGULAR PERTURBATION METHOD. PART 1*

BY

N. D. FOWKES Harvard University

Abstract. An approach to ingular perturbation problems, introduced by Máhony [1] which arose out of the consideration of a problem mvolvmg a boundary layer is applicable to other singular perturbation problems. It lends itself particularly well to problems involving wave propagation, where "multiple scales" are involved. In this paper and the paper to follow, interest is centered around the equation

- g(x)t' = 0,

where is a small positive parameter and g(x) is a bounded function of x which vanishes along simple closed curves in the solution domain. The, one-dimensional case (the Langer turning point problem) is considered iii this paper andit will be shown that the approach leads to exactly the same results as obtained by Langer and his associates usmg a "related equation" method. .

As far as the author is aware, this is the only nontrivial case in which it can be dis. played rigorously (a Ia Langer and others) that a multiple scaling approach leads to the correct result. The apprbaàh employedhere, unlikethe "related equation" technique, is available in nonlinear and partial differential equation applications, and with this work as background one may be confidently hopeful that the approach wifi also lead to the correct resiilts.in these more complicated problems, and in particular in, the partial differential equation application of Part II. In the remaining section of this paper' teäh-mques which will be used in Part II are developed in the simpler context, again m the hope that confidence will be gamed m the vandity of the results obtained using these techniques.

1. Introduction. If the solution of a differential equatiOn with parametric dependence on a small number differs significantly from the unperturbed solution (i.e., the solutiOn obtained by equating the parameter to zero), then the problem is referred to as a singular perturbation problem In the particular case in which the small parameter appears as a coefficient of one of the highly differentiated terms of the equation, the solution may exhibit oscillatory

(eg. sin x/e) or exponential (boundary layer) type dependence

(eg. exp - Ixj /) on the small parameter, and entirely different solution procedures have been developed to cope with these two cases. In both these cases' a straightforward perturbation approach1 fails to yield an interpretable result because the small parameter plays the role of a second solution scale m the problemnot just a solution perturbation parameter. The boundary layer case is essentially the degenerate case in which the second scale is felt only near the boundaries. It will be seen that one approach recently developed by Mahony enables one to deal with both these cases.

*Received February 22, 1967. This research was supported in part by the Office of Naval Research under Contract Nonr-1866(20), and by the Division of Engineering and Applied Physics, Harvard

University.

'y(z, ) y,(x) + e yi(x) ±

58 N. D. FOWKES [Vol. XXVI, No. 1 The one-dimensional wave propagation problems represented by the equation

g(x, )y = 0,

(1.1)

where is a small positive parameter, form a class of oscifiatory type problems fOr which suitable techniques have been developed. The methods employed, however, have leant strongly on the linearity and one-dimensiOnality of this equationrequirements which are not necessary for the application of the method to be employed here. Sturm-LiOnville theory indicates that the solutions are highly oscillatory in a domain in which g(x, ) is negative, and exponential in a domain in which g(x, e) is positive In the W K B J approach solutions of the form exp (j(x, )/3/2) are sought, and a direct substitution of t!iis expression into equation (1.1) leads to the W.K.B.J. approximations

exp g"2 dx). (1.2)

Solution difficUlties might be expected in the neighborhood of a zero ofg(x,

), which

has been referred to as a turning point, because of the radical change in the solution form across this point Evidently the solution approximations (1.2) cease to be valid in the neighborhood of such a point, because of the singularity caused by the g"4 factor-Eq. (1.1) is egi1ar at a zero of g(x, ), so the exact solution must be finite there. More generally, any point at which the W K B J approximations cease to be valid solution representations has been referred to as a transition point Methods have been devised for obtaining connectmg relations between the W K B J approximations valid on either side of a transition point. Modern work on Eq. (1.1), however, has been based on a series of papers by Langer [2], in which he expounded his "comparison" or "related" equation techmque Langer's approach is rather different from that of earlier authors Instead of trying to connect together the W K B J approximations across transition points, he sought a single asymptotic expansion valid throughout the solution domain Transforma-tions of both dependent and independent variables are sought to convert the given equation into one differing only a little from a convement "comparison" equation whose solutions are known If the' difference" between these two equations is sufficiently small, a solution of this transformed equation can be shown to be asymptotically approximated by a solution of: the comparison equation.

In the field of hydrodynamics, nonlinear partial differential equations are encountered, and the solutions exhibit boundary layer dependence on the small parameter Prandtl recognized the need fOr a stretched co-ordinate solution description in the boundary layer. Kaplun and Lagerstrom [3], Proudman and Pearson [4], first devised a suitable systematic method for connecting. the boundary layer (or inner) expansion to the ex-pansion valid throughout the remainder of the field (the outer exex-pansion) This they did

by rearranging and comparing the terms of the boundary layer and outer expansions a limiting procedure being employed to determine the relative orders of magnitude of these terms By removing from the boundary layer expansion the terms in common with the outer expansion, and adding the resulting expression to this outer expansion, a "composite" expansion describmg the solution behaviour throughout the solution domain may be obtained. The fact that the matching technique employs a limiting process to compare the relative magintudes of terms, restricts the parameter range for which the approximations obtained are useful. For example, a term of order exp (h/) is judged small in comparison with a term of order , when the limit as -, 0 is used as a basis

'Often referred to asthe "multiple timing" method.

1968] A SINGULAR PERTURBATION METHOD-PART I 59

for comparison. Numerically, this may not be so for values of for which it is desired to apply the solution expansionespecially if h is rather small. Since the limiting pro-cedure must be discarded to overcome this difficulty, the concept of the two expansion regions must also be discarded, and a single expansion valid .throughout the complete solution domain must be sought. The terms of this expansion will also be required to decrease for small not only decrease in the limit as 0 as before. With this in mind Mahony [1] introduced a general singular perturbation approach which enables one to obtain an apparently uniformly valid expansion and applied this approach to a problem arising from the theory of large deflections of an elastic plate. Firstly, the composite expansion is obtained. An attempt is then made to extend the domain of validity of

this local expansion. This is done by heuristically "summing" terms of the expansionwhich would otherwise grow and thus give rise to trouble if the range of interest were extended. Considerations of this kind suggest a suitable expansion modification. The expansion form thus suggested is formally substituted into the equations and, if necessary, the "summation" process is repeated and further modifications made in this way until the complete solution domain is covered. In appropriate cases the summation process indicates the necessity for the inclusion of further scales in the solution description; in other cases more subtle changes are indicated. Cochran [13], and Cole and Kevorkian [14] independently (in a different context) advanced the "multiple scale" idea and arrived at essentially the same "multiple scale"2 method as that thrown up by Mahony as a biproduct of this approach. It will be shown here that this approach leads to useful

solution expansions in oscifiatory type problems.

Both the Langer and the Mahony approaches lead to uniformly valid asymptotic expansions and therefore depend either directly or indirectly on the finding of a suitable dominating differential equation. In the Langer approach, the dominating equation is guesseda sometimes difficult task even in ordinary differential equations, a virtully impossible task in other cases. In Langer's works [15] "The formalism must be pursued inventively, whereas the lines of the concluding rigorous analysis have been pretty well

laid down." The approach used in this paper enables one to deduce the dominating equa-tion. Furthermore, this approach is available in nonlinear and partial differential equation applications. In the simple turning point problem it will be shown, in Sec. 2 of this paper, that this approach leads to the same dominating equation and subsequently to the same result, as obtained by Langer's method. In the remaining section (Sec. 3) the multiple transition point problem is considered briefly and an approach which is needed in the extension work of Part II is discussed. Also included in this final section is a brief discussion of a problem which arose out of a stability question, and which is also in-teresting here because it is related to the "helium molecule ion problem". The main reason for including a brief discussion of this problem here is to indicate the ease and quickness with which one can determine the dominating equation using the above ap-proach.

2. The single, first order turning point problem. Suppose g(x) vanishes at one point, a, a simple zero, in the solution domain. Also forthe present g(x) will beassumed ana1rtic

in the solution domain which will be assumed finite. Both these conditions may be relaxed considerably in the final result.

60 N. D. FOWKES [Vol. XXVI, No. 1

is a special point in the solution domain. It seems natural, therefore, to commence by looking for solutions in the neighborhood of the transition point a. Now g(z) is analytic and vanishes at a, and therefore can be expanded in the form

g(z)

E g(x

h (2.1)

about x = a. For convenience it is assumed here that g1 = 1, a condition which can be arranged by a suitable choice Of scales.

A balance is achieved between the two terms of Eq. (1.1) in the neighborhood of x = a if significant variations of y occur in a zange of z of order . This suggests the introduction of the stretched co-ordinate i = Cr cr)/, in terms of which the expansion

(2.1) becomes

g(x) = E tngk

Equation (1.1), when expressed in terms of becomes

- ni = y E

kIgknk _(2.2)

The expansion form,

Y

Yo('7) + y() + ...

(2.3)

suggested by Eq. (2.2) certalnly cannot be expected to be useful for nonsmall values

of (x a) because the expansion (2.1) for g(x) is useless for such values, however, such a solution expansion will have a small range of validity and, therefore, will be an ap-pràpriate one to commence an investigation. An examination of this expansion will sug-gest a modification which will enable an extended domain to be covered:

The first few expansion terms of the "exponentially small" solution are given by, = [AoAi()] + [(g2Ao/5)(2Ai'(17) - 17Ai()) + A1Ai(17)] + 2[.

..]

± ...

(2.4) where Ai(17) is the Airy Bessel function which is xponentia1ly small for large positive ,. Higher order terms in the expansion contain still higher powers of i It is clear from, (2.4) that the expansion (2.3) ceases to be useful for values of greater in order than (1/12) (corresponding to a range of Cx

- a) of order

1/2) In order to obtain a first term, Yot which would dominate the solution behaviour for values of (x - a) of unit order it would be necessary to add to the Yo above, terms occurring later in the present expansion for y which, although unimportant for small values of (z - a), assume importance for values of (x -'- a) of unit order. To obtain some idea of the type of adjustment required the sum of the first two terms' of the present expansion (2.4) is considered.

Now, because overestimation of the order of a term is not serious whereas under-estimation is, 172Ai'(17) is better represented in the form

2(x - a)2Ai'() over a range

of (x a) of unit order. Thus for values of (x - a) of unit order, the sum (2.4) is better represented in the rearranged form

Yo + Ey = AoAi() ±

)2

Ai'(i7)'-

(x - a)4i() + EA1Ai(fl).

(2.5)

19681 A SINGULAR PERTURBATION METHOD-PART I 61 This expression can be written in the "summed" form

Yo + Yi = {A0 2Ao(

_{a) + O(x}

-2}

J.(xa)

_Az

(za)2(x)81)

+g2

C 56 6 _JJ (2.6) i.e. Yo + Yi = A(x)Ai(u(x)/t) + O(), where

u(x) = (x - a) + g,[(x - a)2/5J +

and

A(x) = A0 - (g2A0/5)(x - a) +

as can be seen by expanding the right hand side of the above expression. This suggests that an extension of the solution domain so far covered by the expansion form (2.3) may be achieved if a solution of the form y = y(, x, 6), where = u(x)/6, (u(x) being suitable chosen) is sought. It should be noted that close to a,

(x - a)/, andtherefore

is expected to differ inessentially from the stretched co-ordinate, , there.

NOTE. The above argument is purely heuristic but its great, advantage lies in the fact that it offers a method for seeking out solution forms which may lead to more useful perturbation expansions. The "sum" (2.6) is not the only one suggested by (2.5); how-ever, since many equivalent useful descriptions of the solution behaviour (corresponding to different "sums") are available, and we only seek one of these here, this need not be of concern.

2.1. The Formal Solution. In the hope that the above argument has suggested a useful solution modification, solutions of the form,

y = y(x,

, ),

(2.7)

where

= u(z)/,

(2.8)

a function to be determined in the solution process, are sought. The function u(x) will be referred to in the work to follow as the "boundary layer" function, to correspond with the terminology used by Mahony.

Now, in terms of the assumed form

y'(x, 6) = E'yr(x,

-, )u + y(x,

,e) and

y"(x, ) = 2y(x,

,

)u -I- '(2y(x,

,

t)u + y(x,

,

6)u,) + y(x,

,6) (2.9)

where' denotes total differentiation with respect to x, and the derivatives on' the right hand side are partial derivatives. On substituting (2.9) in Eq. (1.1) the following equation for y(x, , 6) is obtained;

y -

1 (2z)y

_{= - { (2y1./u + yu/u) + e2y/u}

Objectionab]e because the expansiOn, y = Yo + yj etc., is not asymptotic for uonsmafl values of

(z -

cx)because of the presence of these terms.

62 N. D. FOWKES (Vol. XXVI, No. 1

The "boundary layer" function was introduced to account for the solution behavior in neighborhood of x= a. In this neighborhood g(x) is of order

,

and therefore (1g(x)/u)

is of unit order. An equivalent (see (2.8)) though better representation (in terms of magnitude) of this term is therefore given by ((/u)g(x)/u). Equation (2.10)

now.be-comes

- k3(x) y

= -

(e(2y/u + yu/u) + e'y/u },

(2.11)

where

k3(x)

g/(uu).

(2.12)

An asymptotic expansion for y(x, , ) of the form

= E

"y(x,

), (2.13)

suggests itself. On substituting this expression for y in Eq. (2.11), and equating to zero the coefficients of the various powers of t, the following equations are obtained;

Yorr - k3ry, = 0 (2.14a)

- 1c3y,.

= -

{(2y_1./u + y_iru/u) + y_2/u}

(for n

1) (2.14b) where y is interpreted to be zero for negative n.

Equation (2.14a) has as independent solutions, Ao(x)Ai(k) and B0(x)Bi(kr), where Ao(x) and B0(z) are, as yet, undetermined functions of x. Because the equations (2.14b) are hnear in y,. , , and Y-2 , these independent leading terms in the expansion for y give rise to independent higher order terms in the expansion for y. The solutions,

exponentially small for large positive L, corresponding to a first expansion term of the form

Yo = A,(x)Ai(k) (2.15)

will be sought. The other independent solutions based on Bi(k) can be found in the same way. On substituting (2.15) into Eq. (2.14b) with n = 1, the equation for y

be-comes

- k3y,

= -

{F(Aok)Ai'(k) + k(Aok4)2Ai(k) , (2.16)

where F('j) denotes

(2.17)

Equation (2.16) is effectively an ordinary differential equation inr,

with x as an

in-dependent parameter. Its solution is the sum of a complementary function and an independent particular integral. This particular integral contains terms of the form

Ai'(k) and Ai(k), arising out of the

Ai(k) portion of this equation, and a term

of the form rAi(k) arising out of the Ai'(k) portion of this equation. Thus, in the

solution for y' , terms in Ai'(k) and Ai(k) arise, which correspond to the objection-able4 terms, i2Ai'() and Ai(), which arose in the heuristic development of Sec. (2.1)

1968) A SINGTJLAR PERTIJRBATION METHOD-PART I 63 in the solution form, for neither u(x) nor Ao(x) are yet determined. This freedom is now utilized to remove these objectionable terms.

By equating k to zero, the rAi(k) term of Eq. (2.16) is removed and the boundary layer function u(x) is defined. For if k = 0 i.e. k = constant, 1 (say);' then from (2.12),

uu = g.

(2.18)

Now, if u is to behave like g close to x = a, (and this is suggested by the heuristic develop-ment of the previous section) then u must be the real solution of this equation. Clearly, for real solutions of Eq. (2.18), u(x) and g(x) have the same sign throughout the solution domain, and therefore u(x) vanishes at x = a, the zero of g(x). This, together with Eq. (2.18) and the requirement that u(x) have a continuous first derivative' determines the boundary layer function u(x), as

u(x) = sign (g(x))(3

r(Jg(x))'2 dx }

2/3

(2.19)

Now that the 'Ai(k) term has been removed from Eq. (2.16) the only objectionable term remaining is the Ai'(k) term and this can be removed by simply requiring that its coefficient F(A0) vanish identically? Thus,

0 U

0 = F(A0) =

- -- + -- A0

(from (2.17)).

This is a first order, linear, differential equation from which A0 (x) can be determined as A0(x)

= (2.20)

where K0 is an arbitrary constant of integration. As far as the particular solution of Eq. (1.1) is concerned, K0 is nothing more than the usual arbitrary multiplicative con-stant associated with the solution of a linear homogeneous differential equation. The evaivation of K0 is therefore to be achieved when suitable complete boundary conditions are imposed on the solution, y. Now since g(a) has a simple zero at x = a, u(a) (from (2.19)) is bounded from zero. Also, u(x) is strictly monotoriic as is seen from its integral representation (2.19) so that (from Eq. (2.18)) u(x) is bounded from zero. Thus, A0(x) and its derivatives are certainly bounded for the analytic functions, g(x), considered in this development.

Now that Ao(x) and u(x) have been determined in the manner described above, all terms on the right hand side of Eq. (2,16) vanish identically so that this equation now integrates to give

= A1(x)Ai(), (2.21)

where A1 (x) is, as yet, an arbitrary function of x. This expression for y fits into the

asymptotic scheme for values of x of unit order and indeed for all values of x in thesolution domain8 provided that A1 (x), when determined, is bounded.

'There is no loss in generality in choosing this constant unity.

is convenient, though not necessary, to remsin on the one solution branch of (2.18) in this way. 7The reader ssi1l recall that Ao(x) is as yet undetermined.

64 N. D. FOWKES [Vol. XXVI, No. 1

The equation for Y2 now reduces to

Y2

- Y2 = - tAi'(r)F(A1)

±

Ai()A'(x)/uJ.

(2.22)

The particular integral associated with Ai() is Ai'(). Thus the particular integral

associated with Ai() term of Eq. (2.22) fits into the asymptotic scheme. Again the Ai'() term is required to vanish identically which requires F(A1) = 0. This determines A1(x) as

A1(x) = Ki/IuI"2, (2;23)

which is bounded as required. The constant of integration K1, again is determined when suitable boundary conditions are imposed on the solution, y. With A1 (z) determined by (2.23), the equation for Y2 can be integrated to give

A2(x)Ai() - Ai'() (A'(z)/u).

(2.24)

The equation for y now becomes

-

Ai'()F(42) - Ai()(A/u) ± Ai"(r)F(A'/u).

(2.25)

The particular integral associated with Ai"() is (2 - )Ai'(), which, for values of

z of unit order, is better represented in magnitude in the form (2Ai'()

-The particular integral arising out of the Ai"() term of Eq. (2.25), therefore, does not fit into the asymptotic scheme. Thus, both the Ai'() term and the Ai"() term give rise to unsuitable particular integrals. The particular integral terms corresponding to these cannot cancel one another in other than a very restricted region because their functional forms are different. Thus the two terms must be coped with independently. As before the Ai'() term can be removed by requiring that its coefficient vanish identically which determines A2(x) as

A2(x) _K2/Iuj"2. (226)

The Ai"() term in Eq. (2.25) cannot be dealt with so simply however. The coefficient of Ai"() in ths equation is completely determined, and in general does not vanish identically. Now,

Ai"() = Ai()9

=

(u(x)/)Ai(),

and herein the trouble lies; for it is now clear that the Ai"() term on the right-hand side of Eq. (2.25) for y is of order e' and, therefore, should have appeared one stage earlier in the solution proceedings (i.e. in the Eq. (2.22) for y2). It is thus apparent that the formal solution process, as described by (2. 14b), has led to an unsatisfactory arrangement of terms and that an alternative arrangement must be sought. Now, since the difficult

term (Ai"()F(4o"/u)) in Eq. (2.25) for ya can be rewritten in the equivalent form

('u(x)Ai(r)F(A0"/u)) this suggests the rearrangement based on simply including

this term as an additional term in the Eq. (2.22) for y2'thereby modifying the schematic solution process as described by Eq. (2. 14b). Now the particular integral associated with

Ai() is Ai'() so that the extra term included in the equation for

Y2 wifi give rise to

an extra Ai'() term ([F(A'/u)u(x)]Ai'()) in the solution fOr Y2

. This, in turn, will give rise to an extra Ai"() term (j1(x)Ai"() (say)) in the Eq. (2.25) for y , which, as before, is rewritten in its equivalent (eT'u(x)f1(x)Ai()) form and is included in the

1968] A SINGULAR PERTURBATION METHOD-PART I 65 equation for Y2 to give rise to stifi another Ai'() term (U1(x)u(x)jAi'()) in the solution for Ya , and so on. Thus extra Ai'() terms in Y2 are generated successively and an

ex-pression for Y2(x, ) of the form,

Y2(X, ₎

A2(x)Ai() + Ai'()(F(A'/u)u(x) + fi(x)u(x) + 12(x)u(x) + ...)

(2.26a)

results. The process of simply transferring large terms (the Ai"() = c1u(x)Ai()

terms) in the equation for y back into the equation for Y2 as described above, therefore, does not lead to a useful result (the successive coefficients of Ai'() in (2.26a) do not decrease). The form of this expression (2.26a) for Y2 , however, suggests that the ap.-propriate "sum" is of the form,

y2(x,

) = A2(x)Ai() + a2(x)Ai'(),

(2.27)

where clearly a2 (x) must be determined in such a way that the successive generation of extra terms in the equation for Y2 is avoided. With Y2 so defined (where a2(x) as yet undetermined) the equation for y becomes,

Y3n -

= Ai()A'(x)/u - Ai"()F(a2),

(2.28)

(where F(A2) has been required to vanish as before) and the term now to be transferred back into the equation for Y2 is Ai"()F(a2). When this term is written in the form 1u(x)Ai()F(a2) and is included in the equation for Y2 this equation becomes,

Y2

- Y2 =

Ai()Ao"/u - Ai()(u(x)F(a2))

(2.29)

= Ai()(A'/u + u(x)F(a2))

Now Y2 as given by (2.27) satisfies this equation if a2(x) is defined to satisfy,

uF(a2) + a + A'/u = 0

(2.30)

By defining a2(x) in this way the successive generation is automatically avoided. Thus by simply modifying the schematic process as described by (2. 14b) in the way indicated. above the Ai"() difficulty is overcome.

Equation. (2.30) is a linear first order differential, equation in a2(x) whose general solution

az(x) = exp (-f p(x)

dx) q(x) exp (f p(x)dx) dx,

where p(x) = (u/u + 1/u)u/2, and

q(x)

(u/2)A'(x)/g(x),

and x0 is arbitrary. An examination of this expression shows that a2(x) has bounded derivatives of all order throughout the solution domain if, and only if, a2(x) vanishes at x = a. This property is required to ensure the higher terms on the expansion for y are bounded. Thus a2(x) is uniquely determined as

a2(X) = exp (-f p(x)

dx)

f q(x) exp (f p(x)

dx)dx. (2.31)

The uniqueness of a2 (x) is to be expected because only one arbitrary constant can arise in the expression, (2.27), for Y2 , and this has already arisen in the Ai(r) term.

66 N. D. FOWKES [Vol. XXVI, No. 1

With a2(x) determined as above, the equation for y now becomes

-

J3 =

(A'(x)/uAi(),

which, in essence, is the same as the equation for Y2 as it appeared in (2.22). The process as it has been described above can therefore be repeated formally ad infinitum since no new type of dependence on will appear on any right hand side. In the application of this method to less tractable equations new types of dependence will appear at each stage of the solution process so that the problem of determining the higher order solution terms will become successively more complicated, but in the present example since the particular integrals will involve only through Ai'() it is convenient to take advantage of this simplicity of form and derive the successive terms by formal substitution of the expression

y =

Ai()

"A,,(x) + Ai'(r)

E "a,,(x), (2.32)

0 2

in the differential Eq. (2.11) and equating the coefficients of 4i() and Ai'() to zero. This leads to the following ordinary linear differential equations for A,,, and a,,which can be solved successively;

F(A,,) + a_(x)/u

= 0

a,,(x) + uF(a,,) + A2(x)/u,, = 0 (2.33)

with a,,(a) = 0 for n 3.

The exponentially large solution is simply obtained by replacing Ai() by Bi() and A's and a's in the above expressions by B's and $'s.

This is exactly the same result as that obtained by Olver [5] who used a "related equation" technique. Olver [5], Langer [2], and others have shown that under fairly general condlitions, this is indeed an asymptotic representation of the solution of Eq. (1.1) which is exponentially small for large positive x. Thus here is a nontrivial case (i.e. not constant cOefficients) where an answer derived by a multiple scaling method has been shown (by an independent indirect mvestigation) to be a genuine asymptotic representation of the exact solution over a significant range in both scales. As far as the author is aware, this is the only nontrivial case in which one has been able to ifiustrate that "multiple

scaling" does lead to the correct result. Thus, with some justification, one may be confidently hopeful that the results detained in further "multiple scale" applications, and in particular the partial differential equation application of Part II will be correct.

The important thing to note about the approach presented here is that it enables one to deduce the correct "dominating" or "related" equation. In the "related equation" approach the "related" equation is guessed. In Langer's own words [15] "The formalism must be pursued inventively whereas the lines of the concluding rigorous analysis have been pretty well laid down." The fact that

y = 0.

with defined by Eq. (2.19), is the appropriate dominating or "related equation" follows immediately from Eq. (2.11) since none of the "leftover" terms in this equation

'°Weber functions arise in the above treatment of the two turning points case.

110ne can only hope to patch the expanthons together up to terms of finite order.

"In the exponential region Ai() is of order exp (fP") while Bi() is of order exp (+ 38l2).

1968] A SINGULAR PERTURBATION METHOD-PART I 67

are greater in order than 6. If one wished, one could at this stage pursue the problem

using the Langer integral equation approach [2], however, if one wishes to obtain an asymptotic expansion, then the formal expansion procedure above indicates clearly the stage in the asymptotic expansion at which care must be taken (e.g. the necessity of the Ai'() term is clearly indicated). The approach indicated in this paper also leads to the dominating equation and asymptotic solution expansions in partial differential equations and nonlinear differential equations.

Since the aim here is to illustrate a method rather than establish a result, it suffices here to state that (although rigorous proofs of this result are not yet available) (2.33) appears to be correct for physically interesting functions g(x) in the infinite domain except in very special cases and point out that the Mathieu equation is one such special case. Mahony [121 in a recent paper discusses this point and in particular he considers the Mathieu equation.

3. Further problems. Two problems which are primarily of interest in association with the partial differential equation extension work of Part II are now considered

briefly.

3.1. Multiple Transition Point Problems. Uniformly valid asymptotic expansions have been obtained by Langer [7] and Kazarinoff [8] using a "related equation" technique in the two turning points case. Clearly such expansions will be useless if they are not in terms of functions with known properties and which are also well tabulated. The limit of usefulness thus defined is reached at the two transition points cases.'° The Airy func-tion expansions of the previous secfunc-tion are valid representafunc-tions of the solufunc-tions of

Eq. (1.1) in a domain containing a single turning point so that the possibffity of describing the solution behaviour in terms of expansions valid in overlapping domains about each of the turning points is worth investigating. Olver [9] has considered this possibility and has derived connecting relations between the two expansions by tracing particular solutions around in the complex plane. However, this technique is not available in the partial differential equation application of Part II. An alternative simple minded ap-proach which leads to the correct results in this case and which is available in the partial differential equation application will be discussed here. With this work as background, one may be confidently hopeful that the approach will lead to the correct results in the later application. Here the question "Under what conditions will patching (matching the solution and its derivative at a point) lead to a useful result" will be answered. Since different solutions may have exactly the same asymptotic expansion in certain regions in the solution domain, it is reasonable to inquire whether there is anything to be gained by being selective in one's choice of the point of patching.

If the expansions are to be patched in the oscifiatory portion of the solution domain (which must be accomplished in the case in which g < 0 for a, < x < a2) where both Ai() and Bi() are of unit order, then errors which are algebraically small in s will be picked up in the coefficients of Ai() and Bi() due to patching.1' The Bi() error thus introduced will swamp the solution in the exponential solution domain if the solution itself is exponentially small.'2 In all other cases the patching error can be shown to be relatively unimportant [6]. The potential well problem is a difficult case in point. In this case, since any portion of a Bi(r) solution (however small) cannot be allowed, one can

"Excluding email neighborhoods of both turning pointsthe Airy function expansion about a

turning point ceases to be valid in the immediate neighborhood of a second turning point.

'4at , = ,, Ai(,) is of order exp (-2/3 ,3'2) so that an expansion of the form [IA(x, ) exp (2/3 F" - 2/3 ,8") Ai (,)} + (B (z, ) exp (2/3 F" + 2/3 Fi") Bi (2)}]is needed to patch onto the

left hand expansion. The result follows from this and patching symmetry.

one extracts out the dependence one has to face up to an equation in which poles and zeros

are closely separated.

68 N. D. FOWKES lVol. XXVI, No. 1

sidestep the above difficulty by ensuring the correct expansion behavior at the two end points ± - The eigencondition requirement,

Ai(,)Ai'(,) - Ai(2)Ai'(,) = 0(s)

where

, =

2 is the patching point, then results from the patching conditions. The independence of this eigencondition on the actual position of the point of patching follows directly from the fact that the left-hand side of the above expression is the Wrouskian of a second order ordinary differential equation. Thus, in the potential well problem the esgenconditions can be established by patching the Ar(r) expansions at any convenient point between the two turning points.13 The most convenient point to patch the expansions together is the point given by = = . The analysisis then simplified

because in this case u1(x) and u,(x) are continuous across the patching point. If the expansions are to be patched m the exponential portion of the solution domam (so that g> 0 for a, <x< a2) then clearly no hope can be held out if the solution is exponentially large because a very large error will be picked up in the Ai() (exponentially small) component, and in the oscifiatory domain where. Ai() is of unit order this error will swamp the solution itself. If, however, the solution is exponentially small then patching results m a correct answer if and only if the "centre"

₌

= (say) is chosen as the patching point. To see this it is necessary to notice first of all that only an exponentially small Bi() term is allowed and errors introduced in this term are important nowhere, and secondly that exponentially large Ai() errors will be introduced if patching is undertaken either to the left or to the right of "centre".'4

3.2. The Dominating Equation. It was pointed out at the end of Sec. 2, that after one has obtained the appropriate expansion form, by either guessing or using the heuristic approach suggested in section 1, one can quickly derive the "related" (or dominating) equation by simply assessing the various terms of the differential equation. To ifiustrate this the equation

s2y

- g(x, s)y = 0,

(3.1)

where g(x, s) /(x) - s(h(x)/x), where f(x) and g(x) are analytic and positive in the domain of interest, is considered. This problem arose out of a stability context but is of interest here because its two dimensional analogue is closely related to a problem that arises when one tries to extract the high energy states of the helium molecule ion.15 Now a heuristic discussion similar to that used in sectiOn 1 leads to the suggestion that a solution representation of the form y = y(, x, s) where = u(x)/s, an as yet unknown function of x, be sought. In terms of the equation for y becomes,

ía

-\

I 2 2

y - at -

= s2y/u + yç

-

y/u

where a = 2f''2/u2 and y = uh/(axu). Now one has to assess the various terms of this equation. The terms that are apparently of order do not, in fact, give rise to large

U

U

)

du

'1/2 1/2

= x,ej1 ax

where S is defined by 8 =' u(w) where w is the zero of g(x, s), and although the term in fact gives rise to particular integrals of order s such terms are too small to influence the dominating equation. With S and defined as above the equation for y becomes,

- (1 - S/)y = s{2y/u +

yu/t! - sy/u}

(3.3)

thus y - (1 - S/)y = 0, with

defined by (3.2) is the required related or dominating equation'° and one can use any one of a number of techniques to estimate the error term. One could, of course, proceed as in Sec. 2 and obtain an asymptotic expansion; however, in this case the process would be extremely laborious.

REFERENCES

J. J. Mahony, Aust. Math. Soc. 2, 440-463 (1962)

R. E. Langer, Trans. Amer. Math. Soc. 33, (1931.), 23 ibid 51, (1937) 669 and many others

S. Kaplun and P. A. Lagerstrom, J. Math. and Mech. 6, 585 (1957) I. Proudman and J. R. S. Pearson, J. Fluid Mech. 2, 237 (1957)

F. W. J. Olver, Phil. Trans. A. 247, 367-369 (1954)

N. D. Fowkes, Ph.D. Thesis submitted 1965, Queensland University

R. E. Langer, Trans. Amer., Math. Soc. 90, 113-142 (1959)

N. D. Kazarinoff, Arch. Rat. Mech. Anal. 2, 129 (1958)

F. W. J. Olver, J. Res. Nat. Bur. Standards B63, 131-169 (1959)

J. Heading, Pha.se.-integral methods (Methuen Moüograph), (1962) M. J. Lighthili, Phil. Mag. 40, 1179 (1949)

J. J. Mahony, Resonance in almost linear systems (to be published)

J. Cochran, Ph.D. Thesis Stanford University

J. D. Cole and J. Kevorkian, Nonlinear differential equations and nonlinear mechanics, Academic Press, N. Y.

R. E. Langer, Boletin de Ia Sociedad Mathematica Mexicana, 1960

leThe solutions of this equation are confluent hypergeometrie functions.

3.2)

1968] A SINGULAR PERTURBATION METHOD-PART I 69

particular integral terms compared with the first order solution if the boundary layer function is given by

ARcHIEF'

Lab1

_Technische

v. Scheepsbouwku

_Hogeschool

Division of Engineering and Applied Physics

Deift

Enaineerin Sciences Laboratory

Harvard Ufliversity

A SINGULAR PERTURBATION _{IIEThOD. PART -II.}

N.D. Fowkes

1937

A SINGULAR PERTURBATION METHOD. PART 11*

BY N. D. FOWTCES Harvard University

Abstract. This paper extends the work of Part I to the partial differential equation

- g(x)

=

0

where e is a small positive parameter and g(x) is a bounded function of z which vanishes along simple closed surfaces in the solution domain. In particular, the eigenproblem corresponding to the case in which g(x) is positive at infinity and in which the boundary conditiOn '

*

0 as lxi is imposed, is considered. One class of eigensolutions is

extracted.

1. Introduction The ideas which have been developed by consideration of linear ordinary differential equations in Part I are now applied to an eigenvalue problem as-sociated with the partial differential equation,

g(x)1' = 0, (1.1)

where x is a vector in rn-dimensional Euclidean space, and s is a small positive parameter. The case will be considered where g(x) is analytic,' positive at infinity, and negative only in a simply connected domain, D, bounded by a closed surface I'. Hereafter, the surface F will be referred to as the transition surface, to correspond with the terminology, transi-tion point, in the one-dimensional case. The boundary conditransi-tion,

' - 0 as

lxi

is imposed, and it is required to find those values of e for which a nontrivial solution to this problem exists. This problem is closely related to the problem of finding the highly excited bound states of the Schrodinger wave equation potential well problem.

There is no need to stipulate that there is a single closed transition surface for, if there are several, the solution can be considered as in the case of ordinary differential equations, 'as the union of solutions in separate regions matched across suitably deter-mined boundaries. It will become apparent to the reader when the form of the present solution is obtained that the other restrictions imposed place no essential limitation on the method. Thus, for example, if g is negative at infinity and suitable boundary con-ditions are imposed on the solution, the same method can be employed to solve the equation.

If (g(x))"2 is thought of as the variable refractive index of a medium, and(e_23)

as the propagation constant which is proportional .to the angular frequency of the field, then the solution, '(x), of Eq. (1.1), represents the standing wave pattern, set up in a medium bounded by the refiexion surface r. The propagation constant (_218) is large, so the problem is a high frequency propagation problem, and therefore, it maybe expectcd that many of the results to be obtained will be closely ]inked with the results of geometric

*Received February 22, 1967. This research was supported in part by the Office of Naval Research

under Contract Nonr-1866(20), and by the Division of Engineering and Applied Physi, Harvard University.

'This condition may be relaxed considerably.

72 N. D. FOWKES [Vol. XXVI, No. 1 optics. Indeed, it will be seen that the rays and wavefronts of geometric optics play an important role in the analysis to follow.

Although the story is quite complicated, the eventual solution technique is straight. forward and is summarized in the conclusion.

Procedural motivation. The same approach is employed here as was employed

inthe simpler one-dimensional case, described in Part I It will be seen that there is a very

close analogy to be drawn between these two cases, and for this reason, only a brief description will be included here of the ideas motivating the formal procedure as dc-scribed in the next section.

As in the one-dimensional case, a description of the solution behaviour close to the transition surface if first sought. In order to do this the equation is described in terms of a co-ordinate system centered on the transition surface. Any co-ordinate systemCr, s),

where r determines distance from, the transition surface and s specifies, uniquely, position around the transition surface in its neighbourhood, may be employed for this purpose,

because this system is introduced only in order to develop a suitable mathematical form for a uniformly valid solution. If the transition surface is suitably regular, the equation for ' in terms of this new co-ordinate system will be of the form

(4' + Va + VV + Vb ±

c) - E gk(s)r' = 0,

(2.1)

where a, b and c are analytic iii s and r close to r, and the symbol V denotes the gradient with respect to s, r remaining constant.

If nontrivial èolutions to this equation exist, they must varyrapidly in some or all of' the (r, s) directions. It seems reasonable to commence by looking for solutions slowly varying around the transition surface, r In this case, the dominating terms of the above Eq. (2.1) will be 4ç, and go(s)r', and a balance will be achieved between these two terms if & varies by unit order in a range r of order . This suggests the introduction of the stretched co-ordinate = r/s, in terms of which Eq. (2.1) for 4' becomes

- go(s)4' = terms

in 2 (2.2)

which is analogous to the equation obtained in the One-dimensional case of Part I. Again, as in the one-dimensional case, a normal perturbation determination for 4' of the form,

=

() + ep(,) ± ...

, runs into difficulty because the terms on the right hand

side of' Eq. (2.2) for 4' throw up terms in 4',, which do not fit into the asymptotic scheme, for large . To overcome this, heuristic ideas analogous to those employed in the simpler case, suggest that an asymptotic solution of the form 4' = 4'(, s, r,

), or, in

terms of x, 4' = 4'(, x, ), where

= 'u(x), be sought.

It should be noted that this solution form suggests itself, only when solutions varying slowly around the transition surface r, are being considered. By seeking solutions varying rapidly in certain tangential directions around 1', it may be possible to extend the class of solutions considered. This, however, will not be done here.

The formal solution. In this sectiOn a formal substitution of the solution form 4, 4'(, x, ), suggested by the work of theprevious section, is carried out. In terms of

the (, x) co-ordinate system,

V4'(x, ) =

'4'1.(x, ,

)Vu(x) + V4'(x,

,

)

and,

1968] A SINGULAR PERTLJRMTION METHODPART II 73

V24'(x, ) e_z{4,(x ,e)['V'u(x)]2}

+ e'{2V4'(x,

,

).Vu(x) + 4'(x,

,e)V2u(x)1 + {V'4,(x, , e)}, where V4,(x, , e) denotes the gradient of 4, with respect to x, remaining constant. A

substitution of this expression in Eq. (1.1) results in the equation for 4' in the form

- (g/[u(Vu)'])4, = 1/(Vu)'{e(2V4,rVu

+

4'tV2u) + e2V24'I, (3.1)

where, as in the one-dimensional case, is replaced by c/u.

This result is rather surprising because it suggests that the partial differential equa-tion (Li) for 4' can be dominated by the ordinary differential equaequa-tion (3d). It should be remembered,' however, that the previous heuristic development applies only to those solutions which vary slowly around the transition surface.

This equation is formally equivalent (with the derivative, did; which occurs on the ordinary differential equation case, being replaced by the vector Operator. V) to that obtained in the simpler case of Part I. One can proceed formally as in the one-dimensional case, or simply recognize that the two cases are formally equivalent and look for solutions of the form

4, = A(x, c)Ai() + 2B(x, )Ai'()

(3.2)

where the "boundary layer co-ordinate"2 u(x) satisfies

u(Vu)2 = g(x) ' ' ' ' (3.3)

with u C) along r. Direct substitution leads to the following equations for A (x, e) and B(x, e),

2(V'AVu + V2uA = 3V2B

' (3.4

and

2u(VB Vu) + uBV2u + (Vu)2B = - V2A,

(3.5)

which are again formally equivalent to the equations obtained by direct substitution in the one-dimensiOnal case of Part I..

It is clear from the earlier work that it is the real solution of the nonlinear hyperbolic boundary layer Eq. (3.3) that is required here. The useful property thatthe characteristics of this equation are the orthogonal trajectories ofthe level curves u = const follows directly

from the characteristic form of Eq. (3.3). Although the characteristic equations for u cannot be explicitly integrated, they can be integrated numerically in any particular case. The question of whether or not Eq. (3.3) has awell-behaved solution does however arise, and this question wiU be discussed later in Sec. 5 when more is known of the signif icance of the "boundary layer" function.

In the hope that the fOrm suggested by the heuristic arguments of Sec. (2) might lead to an evaluation, or at least a simplification, of Eq. (1.1) for 4', certain formal manipulations, resulting in the relationships (3.2), (3.3), (3.4), and (3.5), have been performed. The question therefore to be asked now is"Does the form(3.2) suggested by the heuristic arguments lead to a straightforward perturbation problem, (as defined by Eqs. (3.3), (3.4), and (3.5)) or have the transformations so far performed merely led to equations whose solutions are at least equally singular in theirdependence on

and

2u(VBk. Vu) + vBkV2u + (Vu)2Bk = -V24k, (4.3)

for all lc 0, where it is understood that B_1 = 0. An iterative procedure can now be employed to determine Ak and Bk for successive values of k. The equations determining Ak and Bk are all first order, linear, hyperbolic differential equations, the characteristics of which coincide with the characteristics of the equation for u. These charaderistic8

therefore form the "rays" of propagation (using the terminology of geometric optics) of the sOlution i'(x). Also, the orthogonal trajectories of the level u curves on the (x, y) plane are perpendicular to the shared characteristics, and therefore mark off on the (x, y) plane the "wavefronts" of the solution 4'(x). In order to obtain a clear picture of the solution prop-erties, it is necessary to change to a new (characteristic) co-ordinate system based on these rays and wavefronts. This new co-ordinate net is defined by

a

= u and,

$ is constant along any characteristic, and the value of $ uniquely determines this characteristic.

Thus, for the moment, the particular choice of $ is kept open and disposable at convenience. Equations (4.2) and (4.3) when referred to this co-ordinate net reduce to their simplest forms,

(AkS''2)a = +1/(28"2)7(Bk_I),

and

(Bk(aS)"2)

a = -

1/(2(aS) 112)y(Ak),

where y is a partial differential operator defined by 3 1 &,'\

3 (13w

y(w)

74 N. D. FOWKES [Vol. XXVI, No. I

?"

In order to answer this question, a close investigation of the equations for A and

B has to be undertaken.

For the sake of simplicity and ease in diagrammatic representation in the work to follow, the two-dimensional case will be considered. It will be obvious to the reader that the underlying ideas in the following are by no means restricted to this case; and that an extension of the method to higher dimensional spaces involves nothing essentially different. In the two-dimensional case the transition surface r degenerates to a line, which will be referred to as the "transition line".

4. An iterative procedure for determining A and B. The form of Eqs. ((3.4) and (3.5)) determining A and B, suggests that asymptotic expansions of the form

A =

E akAk(x) and B E 3kBk(x) (4.1)

be sought. A formal substitution of these expressions for A and B in Eqs. (3.4) and (3.5), yields

1968] _{A SINGULAR PERTURBATION METHOD-PART II}

75

S(a, 5) = h2/h1 , (4.7)

where h1 and h2 are scaling factors, corresponding to the new co-ordinate net, defined by

h da2

_{+ h d$2.}

_(4.8)

In an iterative determination of A,. and B,. Eqs. (4.4) and (4.5) become, effectively, ordinary linear differential equations in _{a, which can be integrated to give as general} solutions,

A _81/2

.10 2S12 dci,

and

B,.

-(aS)112

I

2(aS)"27(A,.) dci +

_(S)"2'

.

where F,.() and G,.(9) are, at present, arbitrary functions of _{9. Clearly, to ensure the} continuity of B,. and A,. close to the transition line r (where a vanishes (see Eq. (3.3)), it is necessary to put G,.() = _{0, as in, the ordinary differential equation case of Part I.e.}

The boundedness conditions on the derivatives of B,. and A,. with' respect to a near.

a

= _{0 are also ensured by this condition,. as can be seen by comparison with the}

discus-sion in the ordinary differential equation case. The solutions for A,. and B,. now become

+ f

da, _4.9)

and

- -

_(aS)

1/2!

1 7(A,.)da, 0 2(aS)12

for all lc 0, where B_1 is interpreted to be zero, and where the F,.(),_{for all k 0,} are still arbitrary functions of 3. In particular,

A0 =

F0($)/S'.

_(4.11)

The formal equivalence of the results obtained here and the results obtained in Part I suggests that the expansion,

= Ai(-) E

e8A,.(a, )

+

e2Ai'() E e3B,.(a,fi) (4.12)

(which follows from (3.2) and (4.1)) is indeed _{an asymptotic representation of the desired} solution except possibly in isolated regions of the solution domainwhere S vanishes or is singular. It will be seen (in Sec. 6) that it is the isolated _{points of misbehaviour of 1S in the} solution domain which determine the arbitrary functions of $ still present in the solution expansion.

NoTi. _{It is particularly significant that the approach presented above leads one} quite straightforwardly into the "natural" co-ordinate

_{net of this problemthe net}

based on the rays and wavefronts.

5. Mapping singularities.

_{It has been established that any iterative solution of}

Eqs. (3.4) and (3.5) for A and B is closely associated with _{the properties of the level}

',(A,.) and y(B,.)areregular at a 0 once this condition is imposed.

FIG. I CAUSTICS

FIG. 2

'It is clear also from the fact that the equation determining u is hyperbolis and the boundary condi-tion u = 0 is to be satisfied along a closed curve r that a unique solucondi-tion will not result.

See [2) for a more detailed discussion.

76 N. D. FOWKES [Vol. XXVI, No. I

curves (the wave fronts) of u and their orthogonal trajectories (the rays). A satisfactory determination of A and B over the whole field by this method is thus dependent on the existence of a satisfactory solution for u at every point in the field. In particular, a single-valued determination of u at each point in the field is required.

At a particular point x, the boundary layer function u(x) is evaluated solely in terms of the characteristic (i.e. "ray") passing through this point so that a multiple u coverage of the (x, y) plane can only be realized in regions of overlap of these characteristics in the (x, y) plane. The transition line r is closed, and therefore, at least within r, it is inevitable that characteristics will overlap and a multiple coverage of the (x, y) plane will result.4 There may also be regions in the (x, y) plane not crossed by characteristics within which, therefore, no u coverage has yet been achieved. Since the characteristics form the "rays" of propagation, these regions by their very nature are simply "shadow" regions so that this situation need not be of concern. A closer examination of the implica-tions of Eqs. (3.3) and (4.7) confirms that the solution (x), as given by (3.2) indeed drops to zero on passing into such a shadow region.5

Some typical overlap situations are now envisaged. The first that springs to mind is the focusing situation in which a band of characteristics converges on and passes through a single pointthe focus. e.g. in the centrally symmetric case g(x) = g(r) characteristics focus on r = 0. The characteristics may overlap to form caustics, a situation which is fflustrated in Fig. 1. Variants of this situation which require no special attention do arise. One such variant is ifiustrated in Fig. 2.

The fact that characteristics overlap so that the (x, y) plane is multiply covered does not prevent use of the (a, ) co-ordinate net in the present problem, however, if the solutiQu (3.2) is to have any meaning, a single sheet coverage must be defined. In the single focus case, the obvious way to produce a single sheet coverage of the (x, y) plane

1968] A SINGULAR PERTURBATION METHOD-PART II 77

is to stop all characteristics once they have reached the focus. Even so the mapping remains singular at the focus. It will be seen later (in Sec. 7) that it is this singularity in the mapping at the focus that restricts the solution behaviour in such a way that the tinknown F0(3) in the expansion for '(x) (see (4.11)) can be determined. In Figs. 1 and 2 above, clearly, to produce a single sheet coverage of the (x, y) plane it is necessary to cut from one solution sheet to another, somewhere between LM and LN. There are, of course, many ways of cutting from one sheet to another to produce a single-valued cover-age of the (x, y) plane. it is important, however, to notice that, although this cutting procedure can be defined in many different ways, the end-point L, of the cut is fixed in the (x, y) plane.

Aside. It may seem strange that a certain amount of arbitrariness is associated withthe positioning of this cut in the (x, y) plane. It should be remembered, however, that in the one-dimensional case of Part I, a similar situation arose. There it was found that the matching point (corresponding, to the cut, in the more general case under considera-tion here) could be placed anywhere between the two transiconsidera-tion points excluding a small neighbourhood of these points. As in the one-dimensional case, therefore, it might be expected that different procedures for introducing the cut produce simply different asymptotic representations of the same solution. This point is taken up again in Sec. 8. it is clear from the above, that by introducing suitable "cuts" in the solution domain a single valued coverage can be achieved (except along "cuts" and at foci) for any g(x) and that mani, different cutting arrangements can by employed to do just this. It is also clear from the above that, having adopted one of these "cut configurations", one at most has to cope with a number of "cuts" and "foci" in the solution domain. In order to be explicit a procedure which leads to a particular cut configuration will now be de-

-scribed, and this configuration will then be adopted for the analysis of the next two sections where the solution behaviour along characteristics emanating from cuts (Sec. 6) and foci (Sec. 7) will be examined. The question of the equivalence of the various cut configurations wifi be discussed in section 8.

Along the transition line, F, u is required to vanish. The characteristics are orthogonal to the level u curves so they are "initially" orthogonal to r. There will be a smallest value of ui at which at least two characteristics intersect (e.g. iu(L) j in Figs. 1 and 2). For values of ui less than this value the characteristic equations for u lead to u(x, y) as a single-valued function of position. For larger values there is no unique determination and one has to decide on which 'characteristics to base further evaluations (e.g. beyond L in Figs. 1 and 2 characteristics cross and there is no unique determination). The decision is made to make further determinations on the basis of characteristics other than those that have so far intersected (e.g. other than the characteristics intersecting at L) and furthermore, as one integrates forward in ui to ignorecharacteristics for further evalua-tions once they have intersected. The cut is, of course, the locus of the intersection points. The area of the solution domain left uncovered is continuously reduced as one integrates forward in ui so that eventually this process leads to a complete coverage. Notice that this particular cutting procedure leads to a co-ordinate net in which u is continuous across cuts.

Some simple illustrative examples are now discussed. In the case in which g(x, y) = a2z2 + b2y2 + c2 symmetry arguments indicate that the level u curves are concentric ellipses so that all characteristics focus on x = y = 0. A more interesting situation arises when there are two or more foci within i-e.g. in the case in which

78 N. D. FOWKES [Vol. XXVI, No. 1

b 1 1

g(x, y)

a +

-

_{- 1)2 ±}

,2]l/2

[(x + 1)2 ± 2j112 (5.1) which is cFosely related to the helium molecule ion problem, there are two foci correspond-ingto the two neuclei at x = ±ly = 0. An application of the above prescribed procedure leads to the "cut configuration" indicated in Fig. 3.

r

FIoYEz4

x

-

u const. curves - characteristics cut

FmtraE3

6. Cut analysis. The scaling factors h1 and h2, are well-behaved in a neighbourhood of the cut, except possibly at foci which are considered in the next section. The problem associated with the cut is therefore simply the problem of patching solutions across the cut and is akin to the "patching" problem encountered in the one-dimensional case of Part I. In Fig. 4, the two characteristics, C and D, meet the cut E, at F. Suppose the solution along C is

= AAi(-1) ± 2BAi'(i),

and the solutionalong D:is

= AAi(r2) + e2EAi'(2).

From the construction of the co-ordinate system, = = (say) at F. Now (x), being the solution of an elliptic partial differential equation, is required to have continuous derivatives of all order throughout the solution dOmain, and in particular along the cut. A necessary condition for this to be true is that ', and its first order derivatives, be continuous at every point F, of the cut E. If fi and are unit vectors perpendicular to,

1968) A SINGULAR PERTURBATION METHOD-PART II 79

and tangential to the cut respectively, then these continuity conditions at F require; #1(F) = L'2(F),

when

= (6.1)

and

i.e. to first order in ,

A4i( = AAi(,

AAi'() cos (n1) = AAi'(

cos (n2), (6.2)

AAi'() cos (sfl) = AAi'() cos (S2),

where (nt,) denotes the angle between l and the direction of increasing _{etc. (Use is} made here of the fact that Vu is continuous across the cut which follows from Eq. (3.3) and the fact that u is continuous across the cut.)

These relationships are to be valid for all points F of E. This requires either:

A=A=O()

or

the cut, E, is perpendicular to both families of characteristicsEqs. (6.2)

are then consistent provided the A, A coefficient determinant, 2Ai() Ai'(), vanishes

or

E is a limit line of the characteristics above and below E. No further con-sistency conditions are required in this case. A cut of this type has not been encountered in the examples considered, and it appears mlikely that such a double limit line will arise hi practice. This case wifi not, therefore, be pursued further.

Case I. Because S is finite at the cut, A and A can only vanish to order at the cut, as a result of F($) and F() vanishing to order . (This follows from (4.9) and (4.10).) This implies that A = A = O() along the characteristics emanating from the cut, which in turn 'implies that A = A = B = B = 0, to terms of order ? for all a. Thus the solu-tion (x), along characteristics emanating from a cut of this type, is zero to terms, of order (a solution which is exponentially small in may exist in this "cut region"). Thus, under no circumstances can a standing wave pattern be set up in the "cut region" which forms, therefore, a region of complete interference in the solution domain. The situation envisaged is this. If the initial value problem was to be considered and disturbanceswere

initially set up in the solution domain, then these disturbances would propagate back and forth along the characteristics, interfering with one another, finally leading to a steady statethe state being sought here. In the present case, the transition surface, which acts as a refiexion surface for the propagating waves, is such that no standing wave pattern can be set up in the "cut region", and the steady state solution vanishes due to annulment.

Case II. A nonzero solution exists in this case, provided one of the eigenconditions

given to first order by Ai() = O(), or Ai'()

= O(), is satisfied.

eFor the equation which .&(x) satisfies is elliptic.

7Later (see Sec.. 8) the further restrictions that need to be imposed on Fk()will become evident.

8A reduction in the amount of algebra is achieved by doing this for the presentgeneral discussion, but may not be desirable in a particular problem.

80 N. D. FOWKES [Vol. XXVI, No. I

be referred to as a reflexion region in the work to follow) only those which suffer in-phase refiexion at the transition surface will survive and will be present in the steady state. All others will interfere and annUl one another. The eigenconditions above may be interpreted as being simply phase conditions for these steady states. If these conditions are satisfied, then the matching relations (6.1) are consistent, and may be solved. As in the one-dimensional case, these relations ensure the continuity of all higher order derivatives along the characteristics (i.e. in the direction of Vu). The solution, 4'(x), is also required to have continuous derivatives of all order in the $ direction (i.e. in

the

u = constant direction).° To ensure this, the Fk(19) for all k 0, which occur in the coefficient terms of Ai(r) in the solution expansion (see Eq. (4.9)), must be of class 6C.T At this stage it might be noted that a focus may be thought of as a degenerate case of this type of cutthe case in which the cut degenerates to a point. It might beexpected therefore, that the focus will lead to eigenstates in much the same way.

7. Focus analysis. An infinite band of characteristics intersect at a focus which therefore maps onto a line (parallel to the 3 axis) on the (a, ) plane. Thus at a focus the scaling factor S .has a first order zero and Ao(a, 3) (=F0(j3)/S"2) is singular. Higher order terms in the expansions for. A and B (see (4.9), (4.10)) are still more singular there. The solution '(x) is required to be well behaved at the focus so the expansion (4.12)

does not truly represent the solution behaviour close to this point. An investigation of the solution behaviour close to such a mapping singularity is now undertaken.

Close to the focus, a0 say, the scaling factor S may be usefully represented by a series expansion of the form

S =

s3)( +

s(/2

+ ...)

(7.1)

where

=

(a - ao),

(7.2)

and (because the scaling factor is, by definition, always O)so(13) > O. [In the centrally symmetric case (g(x) =g(r))so(3) reduces to a constanti. Similarly since a is well behaved in the neighbourhood of the focus a power series expansion of the form

a =

a ± u1(3) + ...

, (7.3)

is assumed. The value of a at the focus is determined by the scaleof g(x), and a change in the scale of g can be effected by simply altering the value of

.

Thus there is no loss in generality in assuming that the value of a at the focus is unity (i.e. a0 i).

The expansion (4.12) for '(x) is invalid close to the focus onlybecause the coefficient

expansions (4.1) are singular there. One way of dealing with this difficulty would be simply to modify the expansions (4.1) for A and B in such a way that they coped with the situation close to the focus. This could be done, for example, byseeking expressions for A and B of the form A = A(x, w, ) and B = B(x, c., ), where c., is a suitable chosen

locally stretched co-ordinate, introduced to cope with the singular character of the Eqs. (3.4) and (3.5) for A and B close to the focus, a0. If a uniformly valid asymptotic expansion is to be sought, this is clearly an obvious way to go about it, because it is

1968] A SINGULAR PERTURBATION METHOD-PART II 81 essential to retain the Ai() dependence in '(x) in order to cope with the situation close to

r.

However, it will be remembered that in the two turning points ordinary differential equation problem of Part I a similar situation arose. There, in order to obtain a simpler description of the solution behaviour, the requirement that the solution be described in terms of a single expansion was relaxed, and two expansions were used to describe the solution behaviour. In this case the Eqs. (3.4) and (3.5) determining A and B, are quite complex, and in order to obtain equations which are at all tractable it is highly desirable, if not necessary, to relax the single representation requirement, and to seek a separate description (the "inner expansion") of the solution i(x), valid close to the focus.

The Inner Expansion. Equation (1.1), when referred to the (a, ) co-ordinate system, reduces to

31[ô (

a

a (ia1

(7.4)

(The relationship h1 = u/g which follows from the scaling factor definition (4.8) and the definition (3.3) for u, is used here.)

Consider a ö neighbourhood of the focus a0 . Within this neighbourhood, S is of order , and u is of unit order. Thus a balance is achieved between the various terms of Eq. (7.4) if ô is of order 3,'2 This suggests that Eq. (7.4) might be better described close to in terms of the stretched co-ordinate r = (ao - a)/&, where 8 In ternis of this stretched co-ordinate the equation reduces to

i[a (_a

a (ia'

Lar \ an

a

\

aj

where

and the expansions (7.3) and (7.1), for a and S, become

(7.7)

and

= rso[1 + &Sn + . . ]. (7.8)

On substituting the expansions (7.7) and (7.8) into (7.5), and neglecting terms of order 8, the following approximate equation for ' is obtained;

i a I a\

1

1 a 11 a#\

f' 0. _(7.9)

It will be remembered (see Sec. (4)) that the choice of $ was left open. Equation (7.9) indicates that a most convenient choice of $ is that defined in such a way that So ($) 1 at the focus. A 3 so defined, simplymeasures the angular displacement of characteristics at the focus. This particular choice will be employed in the work to follow. Equation (7.9) is separable, and the solutions which are bounded at r = 0 are of the form

K(â, )f(y$)J?(n) (7.10)

where f('y3) denotes

csin(-y$+d)

(7.11)

= 0,

(7.5)

82 N. D. FOWKES [Vol. XXVI, No. I where c and d are arbitrary constants of integration, and -y is the separating constant whose values are thrown up later when appropriate boundary conditions are imposed on the solution 4'(x). It will be seen that the outer solution (represented by (4.12)) on reaching the focus is not of unit order in , so that it is necessary to include the scaling factor K(o, ) in the above expression for

Matching. The required solution, '(x), is represented by the "inner expi.nsion" whose first term is (given by (710)) in the neighbourhood of the focus a0 , and by the outer expansion (4.12), away from this point. The Lagerstrom, Kaplun matching by rearrangement technique is now employed to cOnnect up these two expansions, and so evaluate the unknown quantities at present associated with each.

For large + ver,

(2/7r)''2K(', )f(yfl)/r''2 cos (r - (ylr/2) - (ir/4)),

(7.12) where use has been made of the asymptotic formulae for the Bessel functions.

Because is of order 8,'2at ao , Ai() can be approximated by its asymptotic form

i.e.

1 1 1 "2 _3/2

-

;i7

IT

cos(r

with the result that the outer expansion (4.12) close to a0 may be approximated by 1 F0(i3) (1/4\

(

ir 2ag"2

175 17r cos

r + -

-where use is made of the stretched co-ordinate relation

r = (ao - a)/ô.

Assuming that there is a nonzero domain in which both (7.12) and (7.13) are valid representations of the solutions '(*), the following matching relations are obtained:

K(o, 6)

(/4

h/2), F0(3) = f(y13)/2"2, (7.14)

and

575 (7.15)

where n is any large positive integer. Thus, the solution along characteristics leading into a focus (except for a small neighbourhood of that point) to first order is given by

- (2S)'

Ai(),

(7.16)

and the standing wave patterns represented by this expression can only be set up if the radial "phase" eigencondition (7.15) (which determines eigenvalues for ) is satisfied to order e. It is more convenient to write this eigencondition in the form;

= ± )'7T/2, where .is determined from either

Ai(-) = 0,

(7.17)

or

Ai()

Ai'() = 0.