THE U N I V E R S I T Y O F M I C H I G A N COLLEGE OF ENGINEERING Department of Meteorology and Oceanography Technical Report ON THE STABILITY PROBLEMS OF THE HELMHOLTZ-KELVIN-RAYLE IGH TYPE C.L. Dolph and F. Stenger Department of Mathematics Aksel C. Wiin,-Nielsen Project Director ORA Project 08759 Supported by: NATIONAL SCIENCE FOUNDATION Grant Nos. GA-841 and GP-6600 Washington, D.C. Administered through: OFFICE OF RESEARCH ADMINISTRATION, ANN ARBOR June 1968

. "'/, ~-:(.b <... 3j;,h, X;,,s

ABSTRACT Differential equations of the form 0" - k 0 = q(y, c, k, 8) 0 subject to the boundary conditions 0(Y1, C) = 0(yZ, c) = 0 constitute a large class of "eigenvalue problems" arising from stability investigations. New methods, a posteriori, in character are developed in this paper which renders a large class of such problems tractable for all values of the parameters k and 8. In Part I, an adaption of the method due to Drukarev which systematically converts an obviously associated Fredholm integral equation is used to obtain a simple expression for the denominator of the associated Fredholm resolvent via a Volterra integral equation. By use of the Weierstrass approximation theorem and/or the Laplace transformation, an elementary solution is shown to be always possible to an arbitrary order of accuracy (involving only monomials of exponentials) provided only that the potential q above be continuous. This will certainly be valid for the regime of instability. Part II is devoted to the asymptotic case when k is large. In it error bounds are given, we believe, for the first time, for the theory based on the corresponding Ricatti equation equivalent to that given above. Part III contains an error analysis for the secular equation and several examples illustrating the methods of Parts I and II. Part IV is based on a modification of Numerov's and Newton's neumerical method and applied to a problem in baroclinic instability. It is the authors hope and belief that the tools developed herein will make many new problems of fluid mechanics tractable since they do not depend on the particular special functions defined by the first equation for any given choice of the potential q

-2Intr odu c ti on The stability problems of the title were initiated by Helmholtz in (1868), Kelvin (1871) and Rayleigh (1894). Here a somewhat wider class, to be referred to, somewhat loosely, as of the Rayleigh type, will be considered. This class arises, as those actually considered by Rayleigh, from elimination after linearization of a system of partial differential equations about a standard known situation, often a velocity or wind profile in a given channel associated with a steady state. In the simplest situations, a solution to the linearized partial differential equation is sought in the form (1) D = exp [ik(x-ct)] 0(y) with the result that the main difficulty is reduced to solving a second-order boundary value problem in order to determine the value of the unknown (usually complex) parameter c which enters in a non-linear way in contrast to that of the Sturm-Liouville theory.. Here it will be assumed that the resulting ordinary differential equation takes the form (2) d 0(y)/dy -k (y) = q(y,c,k, 8) 0(y) which is to be solved subject to the boundary conditions

-3(3) 0(y~) =(Y) 0 In the above differential equation, k is usually a wave number, e stands for one or more other parameters, and c, the so-called "eigenvalue" is to be determined. Most work has been done for the actual Rayleigh type in which q is assumed to have the special form (4) q =Cd~i~~l [ wd 2w(y, c) [d w(y'c) ]/ [w(yc) - cj, dy the dependence on k being understood. Many examples of this type can be found in Drazin and Howard (1), an excellent survey article to which we shall frequently refer as DH, while examples of q's not of this type can be found, for example in Derome and Wiin-Nielsen (2). The'"flow" is said to be unstable if Im c > 0, neutral for Im c = 0 and stable for Im c < 0. Since c enters non-linearly in eigenvalue problems of this type, very few can be solved exactly. In general the ordinary differential equation (2) is not reducible to one having a solution expressible in known or tabulated functions. It is therefore necessary to resort to approximate methods which may be purely numerical, literal, or a combination of both. * See, however, the introduction to Part II -* All references not explicitly given can be found in (DH).

-4The approach here will stress analytic methods although one purely numerical method and example will be discussed in Part IV. In all the other cases numerical steps may be preferred for the solution of the derived secular equations but,even if this is the casethey will occur at the last step. In fact, this paper originated in an attempt to develop methods which would: (1) be significant for such areas as meterology; (2) not required detailed special function theory; and, (3) would be easily understood and easy to apply. These apparently incompatable goals - similar to those of most applied research papers - can be met surprisingly well in two ways. The first way which constitutes Part I of this paper adapts an observation due originally to the Russian physicist Drukarev (3), embellished somewhat more by the physicist Brysk (4), and then discussed rigorously by Aalto (5), (6). We will consequently refer to it as the DBA method although our adaption will be novel in that,we believe,it is the first tirme this set of ideas has been applied to other than a standard eigenvalue problem. In our adaption, in common with the DBA method, a natural Fredholm integral equation is converted into a Volterra integral equation in a systematic way. The next step consists in deducing a closely related Volterra integral equation whose solution when evaluated at the other end-point of the channel yields the secular equation. That is, it is in a real sense, a Volterra integral equation for the denominator of the resolvent of the original Fredholm integral equation. There are of course other ways of reducing boundary value problems to

-5Volterra integral equations. In fact Tricomi (7) for example, attributes one such method to Fubini and illustrates another for the transverse oscillations of a bar. Aalto (5) has subsequently treated this last example by the DBA method although, of course, the resulting eigenvalue problem is still linear in the eigenvalue. Common to both, as well as here, is that the zeros of a transcendental equation have to be found for the determination of the corresponding eigenvalue. Aalto also notes that the DBA method is somewhat related to what Henrici (8), pp. 345-46 calls the "shooting method". To our way of thinking, such methods of search which, in the end, rely on the experience and judgement of the investigator, should be grouped under the general title of biblical methods after the sermon on the mount - "seek and ye shall find" (9). In contrast,our adaptions and modifications yield a method which not only furnishes a roadmap but an actual railroad leading to the solution to within a desired accuracy, and thus we deem it superior to other methods using Volterra integral equations, In Part II our aims are achieved by asymptotic methods for large values of!kI. Its novel feature includes the, we believe, first derivation of error bounds resulting from methods based on the Riccati differential equation. This part has been considerably expanded as a result of Professor K.M. Case's observation that the Olver method as originally presented by the first author could lead to inconsistencies for these non-standard eigenvalue problems. In Note, however, Part IV.

-6Part III, which consists of a discussion on error bounds for the secular equation and examples, we hope that we have answered Professor Case's criticism by presenting a series of results which are sufficient to guarantee the validity of the asymptotic methods for this class of problems. In Part IV a numerical method is given for solving the problem (2, 3), and as an example, the eigenvalues are obtained numerically for an equation related to a baroclinic quasi-geostrophic modelof the atmosphere. Admittedly it has not been possible to obtain the objectives stated above without some penalities. To accomplish them we have had to develop a theory which has an a posteriori character. The main assumption involved is as follows: The potential q in equation (2) is assumed to be continuous throughout the channel. This will certainly be true when, as must be verified a posteriori, the eigenvalues have an imaginary part which is non-zero. If this condition is not satisfied then it may still be true,as it is in the Rayleigh case,if the real part of the eigenvalues lie outside the interval defined by [w in' w ]. [ This is the counter positive statement to Rayleigh's result that if c. = Im c t 0 then c = Re c 1 r must lie in the interval win < c < w ]. rain r max In addition, to obtain uniform estimates,, it will usually be necessary in most instances to uniformly bound I c away from zero. Since c is a function of k, 9 there will be a region in the complex c-plane in which this is true, but again this region can only be found a posteriori.

-7From a mathematical point of view, it is not necessary to impose the above assumption on q since Langer-Olver theory involving turning points and singularities could be used. We have avoided this for we felt that it would prevent two of our three over-all objectives from being achieved. From a practical point of view, our assumption on q on the other hand is probably too restrictive except for the jet-flows considered in DH. All of our examples will involve potentials q which, as long as Im c 0, are analytic in y, c and all other parameters e. Historically, Friedrichs seems to have been the first (1941) to successfully relate Rayleigh type problems with q given in the form (4) to Sturm-Liouville problems. Subsequent work by Tollmien (1935) and Lin (1945) has led to a method,which we shall call the FTL method,which is capable of giving much information about the neutral curve and the unstable solutions in its neighborhood. We will use some of the results of this theorem in Part III and refer to DH for a full discussion. It may prove useful to others to note that recently Derome and Dolph (10) have found another way to relate eigenvalue problems of the type of this paper and to non-selfadjoint classical eigenvalue problems having a "zero" eigenvalue in the classical sense. This has permitted the development of a perturbation theory which has succeeded in clarifying certain aspects of the Eady model of the atmosphere, since it takes into account North-South effects in the wavestructure and has produced predictions more in agreement with observation.' M. E. McIntyre in a preprint entitled "On the non-separable baroclinic parallel flow instability problem" apparently, made the same observation for a siniilar purpose.

-8In contrast, the DBA method when applied to standard eigenvalue problems seems to offer mainly pedagogical advantages in contrast to its application here where it yields a secular equation,to an arbitrary order of accuracy,which involves nothing except exponentials and monomials. However, it should not be forgott~enthat it is, in reality, at the very least, a "pedagogical gem" and the authors are most happy to join Anselone and his student Aalto in their effort to get it more widely known. [ The authors find it hard to believe that this observation dates merely from 1957 and not from the time of Euler or Gauss and find it very surprising,'indeed, that it is so little known five years after the appearance of Brysk's article in the Journal of Mathematical Physics. Could it be that there is also a communication gap?] With the aid of the DBA method there would seem to be no reason why the Fredholm theory, in adequate generality for most engineering applications, should not be taught at the Junior level. More specifically, in addition to the second order case considered here Aalto has shown what modifications are necessary for an n-th order differential equation possessing a Green's function or for a system of linear matrix differential equations which are, at least,as general as that system described by problem 16, page 204 of Coddington and Levinson (11). There is even the possibility of extending the method to the semi-infinite interval if due care is taken with the singular problem which would then arise and perhaps even to the doubly infinite interval for potentials of compact support. Certainly its implication and usefulness for higher-order instability problems such as that of the Orr

-9Sommerfeld equation are worthy of investigation The rigorous mathematical theory provided by Aalto in the Banach space of continuous functions is just one of many possible frames of reference to which the DBA method could be used. As an alternate and one possibly more familiar, to the undergraduate, it would be possible to work in L2[yl, y2] and use Tricomi's almost uniformly convergent series in order to discuss the convergence of the successive approximations necessary to the solution of the Volterra integral equation. In fact,since the Neumann series always converges in any structure reasonable for applications, one might as well assume that W is a world of Leibnitz in which a Voltaire space V exists and then proceed with the method. Such a procedure seems preferable to us since we feel that any mature mathematician will easily be able to furnish an appropriate framework for his needs and by avoiding the detailed considerations of any one of the methods, the procedure will be more accessible to others. Before beginning, what might whimsically be called, a little backward excursion into the 19th Century, it should be noted that care must be exercised if spurious roots are to be avoided. In common with the series method of Arnason (12) or the digital method of Haltiner (13)care must be used that one, after an initial approximation, follows only the roots of interest. In general since non-linear effects will take over and destroy the validity of any linear theory shortly after instability has set in, it will usually be the case that one is only interested in the first complex root that occurs, usually as a result of a quadratic approximation. A However, when in the frequent occuring.case when q(.....c) is analytic

-10 - in c, the Rouche-Horwitz theorem as discussed in Part III can often be useful, while for the Rayleigh form (4), one can use the Howard (1961) circle theorem and the FTL method as also stated in these references.

-11Part I An Adaption of the DBA Method. I. 1. The integral Equation for the "Denominator!" of the Fr edholm Resolvent. As indicated in the introduction, the boundary value problem (2), (3) is readily converted to an integral equation via the Green's function for the operator of the left-hand side of (2). By standard methods this is readily found to be G(y, s;k) = [ sinh k(y2 -y>) sinh k(y< -y1)] /k sinh k(yz -y1) The corresponding Fredholm integral equation can be written successively as Y2 (6) O(y,c, Ok) = f G(y, s;k) q(s;c, k, 9) 0(s;c..) ds Yl + fJ G(y, s;. ) q(s;c,..) (s;c,.) ds Y [ sinh k(y - y2) sinh k(s -yl) - sinh k(y - yl) sinh k(s - y)] q( s, c,. ) ds Y1 k sinh k(y2-yl) Y2 + sinh k(y-yl) S [sinh k(s-yY2)]/[k sinh k(y2-y1) ]q(s, c,.) 0(s)ds Yi

-12 - fY [ sinh k(y - s)/ k] q(s;... ) 0(s.. )ds + sinh k(y - yl) F(0) yl where Yz (7) F(O) - [ sinh k(s -y2) q(s;... )) 0(s)ds] / [k sinh k(y2 -Y1)] yl In (6) the transition from the first equal sign to the second is valid whenever the appropriate integrals exist and have the property that Y2 Y Y2 f' ( )ds - f ( )ds - f ( )ds Y1. Y1 Y1 It is the invalidity of this relation which seems to prevent an easy extension to the doubly infinite interval since at least one of the integrals would fail to exist in general. If q, were, however, of compact support this would not interfere formally but a singular problem might result. The transition from the second to the third equal sign in (6) involves the computation of the difference of two Green's functions and is obtained by two applications of the standard addition or subtraction formulas for hyperbolic functions. The end result is hardly surprising since its kernel would be the natural one to obtain if (2) had had an inhomogeneous term and were to be converted to a Volterra integral equation directly. Assuming as mentioned, in the introduction, that one works in a mathematical structure in which the inverse to (6) exists, we proceed to

-13 - formally derived the equation for the' denominator' of the Fredholm resolvent of the first line of (6). Setting V = fY [ sinh k(y - s)/k] q(s;... )( ) ds Yl equation (6) can be written as (8) (I - V)0 - sinh k(y - y) F(0) -l In any structure in which (I - V) exists this can be written as 0 = (I - V)- [sinh k(y-yl)] F() Application of F to this yields [I - F(I - V) sinh k(y- Y)] F(0) 0- As first asserted by Brysk (4) and first rigorously established by Aalto ( ) the quantity d(c;k, 8) = { I - F[ (I - V) sinh k(y - Y1)] }

-14is the denominator of the resolvent of the Fredholm integral equation (6) We proceed to examine its structure, Let -1'(Y) = (I - V) sinh k(y -Y1) then p(y) will satisfy the Volterra integral equation (Y) - JY [ sinh k(y- s)] / k] q(s;... ) i(x) ds = sinh k(y - xl) Y1 and from (7) F[ (I-V) l sinh k(y-Y1)] = F [I] f [ sinh k(s -y2) q(s;.. ) i(s)ds]/ [k sinh k(y2 -yl)] 1 Setting (9) c(y) = fY sinh k(y- s) q(s;.. ) i(s) ds Yl it follows that (10) d(c;k, 8) = [I-F(p)] = I + o(YZ)/ [sinh k(Y2 - Y1)]

-15 - However, since ~L(y) satisfies (8), it follows by multiplying (8) by [sinh k(y -y')/k] q(y';. ), integrating on y' from Yl to y and using the definition (9) that cr(y) satisfies (11) o(y) = fY [ sinh k(y - s)/ k] q(s,... ) cr(s)ds+fY [ sinh k(y - s)/ k]. Y1 Yl q(s,o..) sinh k(s -yl)ds In summary then, to obtain a Volterra integral equation for the "denominator" of the resolvent of (6) it is merely necessary to solve equation (11) for o(y) and then set y = Y2 in (10). Alternately stated if one wishes to define D(y;c,k, 6) = I + (y)/[sinh k(y2-y1)] so that D(y2;c, k, ) = d(c, k, 8) then the quantity sinh k(Y2 -Y1) [D(y;c, k, 8) - I] would satisfy (11) - the integral equation for the "denominator".

-16 - While any of the above can be used better error bounds result when (11) is replaced by an equivalent equation having a kernel bounded by 1/ Jkj. Setting o(y) = eky W(y), equation (11) can be replaced by (12) W(y)=:Y [(!ek(Y s))/2k] q(s) [W(s) + (1 - e )/2k] ds Y1 From this it follows that in place of (10) 9 c can be determined as a root of the following: (13) (Y2' c) = 2W(y2) + l -e Zk(y2Y1) 2 i, 2. The effective determinations of the''denominator " of the Fredholm resolvent. For any given value of c, the Volterra integral equations (11) or (12) can be, in principle, solved to any order of accuracy by the Neumann series. To obtain these solutions, even in principle, one must have a bound on the potential q independent of the yet unknown eigenvalue "c". This can be achieved by restricting attention to a region of the complex c-plane. In typical cases such as those treated by Rayleigh this will involve limiting "1c1t to a region in which its absolute value is uniformly bounded away from

-17 - zero. That is since IqI q = Iw(y)I / Iw(y)-cl M i ci where Iw" (y)| < M to obtain a uniform bound it is necessary to select an n. such that Icl > Icil' >.. The explicit region c = c(k) [ c(k, 9) in general] in which this assumption will be valid can only be verified a prosteriori just as the assumption that q is continuous must also be verified. Assuming that q is continuous and has been uniformly bounded as a result of some such restriction as that given above for "c" in the Rayleigh case, it is possible to replace the difficult task of computing the terms in the Neumann series for (11) or (12) by replacing these equations, to an arbitrary high order of accuracy,by equations for which the corresponding series involves nothing more complicated than successive integration by parts of exponentials and monomials. This is accomplished by use of the Weierstrass theorem,since under the above assumptions q(y,c,..) can be uniformly approximated in both y and c by polynomials. In the sequel we shall assume that an N(r) has been determined so that for

-18 - N P (y;c, 9) -= q.(c, e) y N j=O J one has Iq(y;c, e) - PN(y;c, 9)1 < (e) The assumptions on q will imply the existence of a bound for PN independent of "c" but the actual nature of this bound will depend on what method is used to determine P. For example one may use the proof of Lebesgue in which polygonal appro'>ilalatiOls are used (cf. 14, p. 219), the method given in Courant-Hilbert ( 15) after a linear transformation of the interval, the method of trigometric polynomials as it is described in Widder (16 ), the method of approximation based on the fundamental solution of the heat equation as it is described in Epstein ( 17 ), or a method of interpolation. If q(y, c) is analytic in y then Taylor's series about (Y2 - Y1) /2 may be the most convenient. In any event the continuity of q enables us to uniformly approximate the integral equation (12) by the equation U(y) fY K(y- s) PN(s) [ U(s) + K(s)] ds Y1 where for simplicity we have denoted the kernel of (12) by K(y). This follows from the difference

-19 - [W(y) - U(y)] = fY K(y - s) { [q(s) W(s) - PN(s) U(s)] + K(s)[q(s)- P (s)]} ds Yl fY K(y-s) q(s) [W(s) - U(s)]ds + fY K(y-s)[q(s) - PN(s)] Y1 Yl [ U(s) -+ K(s)] ds This is a Volterra integral equation for e(y) = W(y) - U(y) to which the generalization of the Bellman-Gronwall inequality can be applied. Since we will use this inequality many times we will sketch a new proof (due to Stenger (18)), For notational simplicity we designate the last term on the right of (16) by fy f(s)ds. It follows from (15) that y1 le(y)| < fY [(1/k) |q(s)l |e(s)| + If(s)l ] ds Y1 since iK(y-s)i < l/k. If p(y) denotes a suitable non-negative function' this may be written as an equality in the form ]e(y) = - Y [(l/k) Iq(s)| le(s)l + If(s)| -p(s)] ds Yl Since this equality is equivalent toa first-order linear differential equation, it follows by differentiation and use of an integrating factor that

-20(e(y)I = SY exp [ fY Iq(s')/(k)ds' ] [f(s) - p(s)] ds Y1 s < exp [ JSY Iq(s')/(k)ds'l] fY If(s)| ds Y1 Yl Since If(s)} can be estimated with the help of (14) it is clear that Ie(y) I can be made uniformly small in y and c by taking N sufficiently large, if the well-known apriori bound for U(y) is used. [For the sequel it should be noted that the same argument will apply equally well if IK(y - s) < F(y) G(s). ] It is therefore clear that if q is continuous it is sufficient to consider (15). Starting with the usual iteration scheme (17) U +1(y) - JY K(y -s) PN(s) [Ur(s) + K(s)] ds (18) U (y)= K(y) PN(y) also it islclear that the Neumann series can be explicitly calculated by integration by parts of terms consisting of monomials and exponentials in y. Thus not only the usual Volterra estimate holds but the Neumann series can be explicitly computed. For example,if one chooses Leibnitz W to be terra firma and the Voltaire V to be either L,(y1, y2) or the Banach space C on the interval [Yl,yZ] reference to either Aalto (5) or Tricomi (7) will >' We are in the space age, da?

-21provide the necessary estimates. Since we are primarily interested in the eigenvalue equation for c, we will restrict ourselves to the latter. Suppose then that a number M has been determined such that after M iterations of (17) one has that IU(y) - UM(Y)I < 1(0) We are now in a position to make an estimate on the secular equation for c which, as will be seen in Part Ill usually implies a corresponding estimate on the eigenvalues themselves. Since S(y2, c, 0) is defined by (13) it is convenient to denote the corresponding expression by SN(Y2, C, 0) when W(y2) is replaced by U(y2) and by SMN(Y2, c. 9) the same expression when W is replaced by UM((Y2) Since IS(y;c, ) - SMN(y2;c, )I < I S(yz;c, ) - SN(YZ- C, ) + IS M(y2;c0) - S MN(y; C, 0) the use of the Weierstrass theorem and the usual Volterra estimates make it clear that the left-hand side of the above can be made uniformly small if q is continuous. Using the same subscript notation as above for S, SN and SMN to the corresponding d's, we note for example that in the uniform operator norm the estimates of Aalto apply equally well to the

-22 - interval I [ Y, Y2]. Thus if q is assumed continuous for complex c in this interval, i" t" denotes the uniform norm, and if one defines: M - max [ sinh k(y- s)/k] q(y, c, k, ) I; the function f by f = sinh k(y-y1); the operator K by the relation KO = fY [ sinh k(y - s)/k- q(s;...) 0(s) ds, Yl and the approximate Fredholm denominator by the relation n d (cI,...) - F[0] [E KPf], then n one can easily derive the following estimate: (17) d(c, )-d(c, )[ < OFO OfO [eM ( MP/p/ )] Here the norms are to be understood as above. The requisite estimate now follows at once from the definition (13) I. 3. An alternate effective method for arbitrary k based on a Sinh calculus. While the method of the last section involves nothing more difficult than successive integration by parts of products of monomials and exponentials,

-23it appears simpler to us to replace this process by purely algebraic operations via the use of the Laplace transform, followed by the inversion process this may take one or several forms depending upon the algebraic scheme preferred by the user. In one form the operational calculus becomes automatic quickly while in others it is easier to establish bookkeeping rules. Consequently we will present several closely related versions which can be stated in terms of operations which apply equally well to (11) or (12) when q(y.,.. ) is replaced by its polynomial approximation. No loss of generality is involved if the interval [ yl, Y2] is replaced by the interval [0,1]. Having accomplished this by a suitable linear transformation,let the Laplace transform of an arbitrary function f(y) be as defined by the relation f(p) = S- e PY f(y) dy 0 Then equations (11) and (15) take the respective forms 2 2 (18) P l(p) = [l/(p -k )] PN(d/dp) [ p (p) + l/(p k)] (19) U 1 (p) = (1/2k) [l/p-l/(p+2k)] PN(d/dp)[Um(p) + l/p(p+2k)], while P = PN(d/dp) [l/(p - k) and Uo(P) = PN(d/dp) [ l/p(p+2k)]

-24respectively. In the above p(p) has been used to denote the unknown in (11) when q has been replaced by PN In view of the simple identities (20) 1/p2_k2) = (20) lip -k) = -' [l/(p-k) - l/(p+k)] (21) l/p(p +2k) = (1/2k) [l/p -l/(p+k)] it is clear that a simultaneous operational calculus for both (18) and (19) can be developed by defining X = (1/2k) and interpreting (20) and (21) as (22) ab = Xa - kb where a = l/(p-k), b = l/(p+k) if (18) is to be solved and a = l/p, b = 1/(p+2k) is (19) is to be solved. As written equations (22) clearly indicates that an operational calculus may be based either on the lefthand side of (22), the right-hand side or both. If the left-hand side is to be used, it should be noticed that it will involve evaluation of terms of the form p L(f) this would seem to indicate that delta functions and their derivatives would occur. The use of an operational calculus based on the right-hand side of (22) however makes it clear that in fact this does not happen and so the usual operational rule of the Laplace transform reduces

-25to the simple expression (23) piL(f) = L[f (y)]. That is, the identity (22) implies that all functions which occur in either (18) or (17) and all their derivates up to the indicated order must have zero initial values. The operator PN(d/dp) which occurs in both (18) and (19) is given explicitly by N PN(d/dp) = Z (-1)j qj dj/dpj N j= since each term in either (18) or (19) is a convolution of the kernel with the product of a monomial with a function. Before systematically developing these operational calculi,two remarks are in order. The first results from the fact that the transform of e.g. (31) leads to a differential equation of the form 2 2 d 1 PN( )p _ (p k )P = PN(d 2 p -k in which all coefficients except that of the derivative free term of the lefthand side are constant. This apparent simplicity is, however, deceptive, since one wants solutions for p corresponding to small y and thus one is,

-26 - in reality, faced with an irregular singular point at infinity. For N greater than one this appears to cause precisely the complications the method here hopes to avoid; namely reliance on special transcendental functions. The second follows from the fact that equation (18) is clearly equivalent to the differential equation,, 2 sinh ky p k-k p - PN(y)p PN(y) sinh ky p(O) = p (0) 0 This suggests the use of the transformation ky -ky p = e u + e v which will be sufficient to lead to a solution if u, v satisfy respectively the equations u + 2ku - PN(y)u = 1/2 PN(u) v - 2kv' - PN(y)v = 1/2 PN(v) subject to the initial conditions u(O) = v(O) = u'(0) = v (O) = 0. Power series solutions toeachof these linear equations will clearly exist and converge everywhere but the relation to the solution known to be correct and unique as obtained by the successive approximation process (22) is

-27 - far from clear. Thus, for example, if PN = q0 + qly a series solution to the fifth power is: u = q0y /4 - 2 + kqy + (q /4 + 3 k q1 + p k2 q)y 2 2 + [-q ql/4 + k(q1 6q -4q + k (12ql) - 8 q0k y... It appears that the convergence would be very slow at best, except for small k To resume the development of the operational calculi, note that the one appropriate for (19) follows from that for (18) if p is replaced by (p - k). Since there is more symmetry in the method when applied to (18) as well as the fact that it involves hyperbolic functions thereby justifying the name of this section, (18) will be used as a basis. Moreover, since d[ (p2 -k2)] -n /dp = -n(p2 - k)(n+) (2p) it follows that the left-hand side of (18) can be written in any approximation in terms of products of positive powers of p and l/(p -k ). Thus once the number M has been reached it is merely necessary to have a way of finding the inverse transform of these. Following up a suggestion of our colleague, Professor J. Ullman, it is convenient to introduce a generating function a la Rainville;

-28 - 00 n 2 2 -n 2p 1 n(p -k) = y/[p (k + y)] 1 Since this is always meaningful if the real part of p is sufficiently large, it follows that L [ z y (p - k ) n] = [ysinh (k + y)]/ [k + ] 1 2 2 -j Thus to find the inverse transform of (p -k ) it is merely necessary to find the coefficient of the j-th power of the right-hand side. This may be computed easily by setting w = k + y 2 and noting that since d(w )/dy = 1, the Taylor series for the right-hand side of the above identity becomes oo 1 d [ sinh (w2)1/2 x / (1/2) j+l j=) dx(w ) y=O Thus the inverse transform of [(s - k ) ] corresponds to the coefficient of yj in this series Once this has been found and denoted by, say, h(x) the inverse transform of pJL(h) is readily found by simplification of the standard rule given by (23).

-29These operations, while tedious, are all elementary, and the work can be arranged systematically once one notes that the equation for (p) [ -M — M-l i p (p) [ (fQ) + I (fQ) +... +(fQ)] f 2 M-1 [ M(fQ) + (M -l1) (f Q) +.. 2(f Q)M + 2(fQ)M] f [ (fQ)M + D (M-k+l) (fQ) ] f k=l It is also convenient to introduce the product functions 1- = a fn m n so that io,0 = f and dH /dp mTT 1 n -n if mm m -1, n m, n+l This implies that (f Q)k f = of Q) 1 Q [ TTo k+l and hence some simplification above.

-30As noted earlier the procedure for solving ( I8) is based on successive use of the identity (22). For this purpose it is convenient to let a = I/(p -k), b = l/(p +k), n-i n-l = 1/(2 k) and to introduce two new combinations r = X a and n m-i rn-1 t =X a b where rl =a, r2 = a-b, t1 = b and t2 = a-b Since d aj/dp = -j aJl d bj/dp = -j aj+l the successive approximation process based on iTT = ab will involve mT n the eventual inversion of a b until the number M is reached. This is to be by accomplishedq use of the identity ab = Xa - Xb in order to keep the successive derivative of a and b separated can be simplified by noting that n-2 n-i r = r -X b n n-l (24) m-2 m-1 t X a - t m m-l n+m-2 m n m-2 m-l n-2 n-i r t = X a b =+X a r +X b t n m n-i m-1 -2r t n-l m-1 so that rt = (a +b)r t -2r t n m n-2 m-2 n- rn-i

-31The authors are indebted to their colleague,Professor T. Storer, for derivation of the last recurrence relationship. Successive use of the above will result eventually in the necessity of inverting a sum containing coefficients of (p+k)-j for different powers of j. After the number M has been reached, the inversion of (25) rn tm = (1/2k)n+m-2 [ (k) (p+k) -n can also be accomplished by noting that n+m -2 (2k) r t = [residue of (r t ePY) at p = k] + n m n m [residue of (r t ePY) at p = -k] n m The first of these is given by (26) 1/(m-l) dm-i [ept ( + k) -n]/dpm-l and the second by (27) 1/ (n-l)! dn-1 [ept (p-k)-m dpn-1 p= -k The amount of work necessary for this procedure will be illustrated in PartIIIfor M= 1, N= 2.

-32 - Experience suggests that a mixture of these two approaches may be the most efficient. That is, one basically uses the combination 2 2 1/ (p - k ) and its powers except when it is necessary to differentiate. It is differentiated as 1/[(p-k)(p+k)]. The identity (25) need be applied only at the very end or the inverse can be found by the standard method just given above.

-33Part II Asymptotic Methods Although the methods of the last section are easy to use, they do involve a great deal of tedious computation if either M, N or both are large. Under certain circumstances the use of asymptotic methods will yield simple estimates for the eigenvalue although as pointed out to us dramatically by Professor K. M. Case care must be used if inconsistency are to be avoided - see part III for details. In this part we shall describe two methods yielding asymptotic results for large k This is not as restrictive as it might appear at first glance since many parameters are frequently lumped into a parameter which we have denoted by k. For example for the linear quasi-geostrophic problem of baroclinic instability, the quantity which would correspond to our k after the application of the Liouville-Green Transformation is actually 2 2 2 k Co p0/ f where or is a measure of static stability, p a standard value of surface pressure, f a standard value of the Coriolis parameter, and k the wave number. [For further details see Derome-Wiin-Nielsen (2)]. Thus there are many different regimes which would correspond to large lIkI in equation (12). For an excellent account of the relationship between the methods of Picard of Part I and the asymptotic methods associated with equations (11) or (12) the interested is referred to Erdelyei (19).

-34II. 1. The Method of Olver. There appears to be little advantage in dealing with (11) or (12) in contrast with dealing with (2) directly. For the latter, one seeks a solution of the form m -1 m -1 (2 8) 0(y) = ( aj/k + + e ( b./k + ) j-0 J 1 j=02 in which the a's are defined recursively by aO= 1/2k a+l -1 [a + Y q(s, c, )a. ds] aj+l 2 b. (1)j+l a. J J Now Olver has established the following: Theorem(II. 1): There exists a solution of the form (28) where Eil < exp [ fJY q(s, c, )/ (2k) ds ] Y la/ (k) ds a. a 1. i 1 provided Re k > 0. Here k =k if k> 0, and k= -k if Im k O a =0 a2 = 1. The following form of Olver's theorem is more convenient for our purposes. We omit the proof which is similar to that of the above theorem.

-35Theorem (II.1'; There exists a solution of (2) of the form ky ky m- k (28 ) 0(y) = ekY [ (a /kp) + E ] + e Z b / kp p=O P pO P where i[ < 2 exp fY Iq(s, c, 9)/ (2k')ds I fY la/ (km) ds 0 0 1 where k'- k if k > 0, and k'= k if I m k 0O To use this method one chooses an m in the above form (28') and then sets 0 (1) = 0, where 0 (y) denotes that part of (28) which remains when one sets E = 0 II. 2. An asymptotic method based on the Riccati equation. In this section we consider the equivalent equation (29) wI + 2 kf (y, c, k) w + k g(y, c, k) w = 0 since the method is unaffected by the presence of a first derivative term, and thus is sometimes more convenient.

-36The substitution [cf. e. g., Cesari (20)] (30) w = exp [k f u dy] reduces (29) to the first order Riccati differential equation (31) u' + 2k f u + ku2 + kg = 0. It is now assumed that for each integer n in the range 0 < n < m + 2 the functions f and g admit the representations n-i ~~n - 1 ~-n (32) f = z f.(y,c) k J + f (y,c,k) k j=0 o n n- 1 (33) g = z g.(y,c) k + gn(y,c,k) k j=0' In the asymptotic developments of this type it is usual to assume that f (y, c, k) and g (y, c, k) are uniformly bounded functions of k for all y on I, where I denotes the interval 0 < y < 1. With such assumptions the solutions we shall obtain below will have the usual asymptotic properties. However, since we also develop error bounds on approximations to solutions of (31), we do not require that f (y, c, k) and gn(y, c, k) be uniformly bounded functions of k0. In fact we may allow f.(y, c) and gj(y, c) to depend on k. This arbitraryness in f. and introduces an additional deree of freedom for the user who and gj introduces an additional degree of freedom for the user who

-37may want to concentrate: (i) on choosing expansions of f and g which minimize the error bounds in a certain sense, or (ii) for purposes of solving the eigenvalue problem (2, 3) to choose an expansion of f and g which, while not necessarily yielding the sharpest bounds, makes tractable the integration of the first few coefficients in the resulting approximation. The equation (31) may be solved formally by a (not necessarily convergent) series of the form oO (34) u = z u.(y)/kj j-=O J From (32), (33) and (34) we have o00 u' = Z u' /kj j=-O kg k g/kJ = gj+l/k 2 00 j+l oo j+l ku = Z Z u. u. j=-l i=O 1 oo j+l kfu = Z Z f. u-i /kj ji=-l i=O () At this stage it is perfectly all right to ignore the fact that (32) and (33) may not have convergent expansions.

-38When these are inserted into (31) it becomes oo j+l E [ul + { (2f. +u.) u } + g /k- 0 j —1 J i-=O i j+l-i gjl ] /k Defining v by the relation 2 1/2 v - + (f - o) one obtains the relations (35) u0 -f0 + v ui+l = (-1/2v)[ui + 2fi+l u + gi+l + (2f + t) i+l-t i t=l t t i+l-t provided v f 0 in I. As stated in the introduction,the error bounds we are about to derive are not only new but appear to be the first that have been derived for the differential equation given in (31) Letting U stand for the sum of the first m-l terms in (34) (where m > 2), we can define E by the relation (36) u = U + E m which, when inserted into (31) yields the following equation for ~ ~

-39(37) + 2k(f +U )E + kE 2 r m m where r = -[U' t k(2 U + U+ + g)] m m m m 2 v u m -1 (38) m-l [ g (y c, k)+ Zuf (y, c, k) rn-i1m m [ i m+l-i k k i=O i +-i 2m -3 m -1 m -juu +z km-j u u].] i j+l-i j=m i=j+2 -m +l The transition from the first equality to the second follows from the expansion of (38) in connection with (35) and (32, 33). It may be carried out in a fashion similar to that leading to (35). We omit the details. In short, we have that r = 2v u /k-1 + terms from which 1/k can be factored. m m Notice now that

-40 - f + U m Of + f/k + f (y, c, k)/k2 + U ul/k rn 0 1 20 1 + (U - -ul/k) (f0 + u) + (fl + )/k +[f (y, c, k) + k(U - - Ul/k)] /k2 0 0 1 1 2 (Urn 0 = v + [2vf1 - u0 - 2 fl(v -f0) - gl]/vk [f2(Y, c, k) + k(U - u - ul/k)] /k2 - - {u0 - g flf0 + gl)/2 v k + [ f(Y, c, k) + k (U - uO - Ul/k)]kZ where (32) and (35) have been used to obtain the last equality. Thus if one defines ji(y, c, k) - vk[1 - (f +gl + v - 2fl f/kv ] and T (y,c,k) = 2Zk(f+U m)- i(y, c, k) equation (17) can be written as

-41-,r~~~ Z~~2 + L(y, c, k) E = r - k ~ - T (y,, k) ~ from which it follows that (39) E(y) Y exp [ J (s, c, )ds] [r (t) - k s] [ (t) - it(t, c, k) E (t)]dt Y0 Y where y0 is an end-point of I We will now assume that the following is valid.: Assumption (II. 1): For all points i, q in the order y,, 71, y the relation Re L (Q, c, k) < Re [t(r1, c, k) holds. Consider that portion of (39) given by (40) p (y,ck) jY exp[ I(s, c,k) ds] r (t) dt. yo Y Now if g(y) is a point on I in the order yo0, t(), y such that |1 -expS p dt 1 exp [ (Y) p dt] > 1 exp [ ] Y Y for all t on I in the order yO, t, y, and provided that r (t)/.t(t, c, k) is absolutely continuous on I, integration by parts of (40) yields

-42YO P (y, c, k) = [ 1 - exp Sf [ dt] r (yO)/I(Yo0 c, k) Y + fY[ 1 - exp ft fdt'] dt (r /) dt Yo Y so that if one notes by assumption (II. 1) that 1- exp[ f dt'] < 2 y one obtains the estimate 1(y) (41) Pm(Y' c, k) < Pm (yc,k)= 1-exp jf dtl rm(Yo)/l(Y, c, k) Yo It should be noted that p (y c, k) = 0 and that p (y, c, k) increases as y traverses I from y0 to y1 on I. If r /p. is absolutely continuous there exists a non-negative Lebesque integrable function pm (y9 C, k) such that Pm (y, c,k) = Y Pm (tc,k) Idt YO Under these circumstances let us make the further assumption that for all y in I,

-43 - (42) SY Ik E (t) dt | < a f Y p (t, c, k) Idt| y m Yo Y0 where a is a constant which will be specified below. Under Assumption(II. 1) it follows that I exp [ f t dt ] I < 1 so that the estimate (41) implies the y following: (43) E (y) I < Y [ T(t, c,k k) I E(t) + (l+a) pm (t,c,k)] IdtI Yo The Bellman-Gronwall lemma implies that (44) i E(y)| < (l+a) exp [ Y I T (t, c, k) dtI ] Y pK (t, c, k) IdtI yo Yo Inserting (44) into (42) yields fY Ik ~ (t) dtI < (1+a)2 jY exp [ 2T (s, c,k) Ids. [p (t,c,k)] |dtl YO Yo YO < [k| y-y I (1+cv) exp [ Y I 2ZT(t, c,k) dt ~ [pm(yk)] Yo

-44so that (42) will be satisfied if (45) a/(1+a)2 > Ikj y-y0I exp [ fJY 2 T(t, c,k)dtl ] p* (y, c,k) Y0 for all y on I Now if 13 denotes the maximum with respect to y of the right-hand side of (45), it follows that 3 < - and so 4 a > (1++a) 1 That is 2 1 1 a -2a( -1) + 1< 0, o < 13< 4 This last inequality will be satisfied for all a such that i 1 ) 1/2 1 1 + 1 2 1/2 21 23 - 2P1 21 We summarize the error bounds in the following theorems: Theorem (II.2): Let r /p. be absolutely continuous on I, where r is defined by (38). Let p (y, c, k) be defined by (41), and let!=|epm |Ztck t|J |< 1=Iklexp f 12T(t,c,k)dt I f pot) I dt < I I

-45 - If the Assumption (II. 1) is satisfied then the differential equation (31) has a solution given by (36) where I|(y)I < ( + a) exp [ fY IT(t, c,k) dt| ] p (y) m Yo with 1 1 2 1/2 1 - [ ( - 1) -1] 2P 2P If we denote by U(t) the expression (46) U(y) k Y [U + E] dt k fy U dt + rl(y).. m y where y' is on I, then the substitution (30) and Theorem (II.1) imply the following: Theorem (II. 2): Let the conditions of Theorem (II. 1) be satisfied. Then the differential equation (29) has a solution (47) w(y) = exp U(y) where U(y) is given by (46) and where

-46 - l(Y) < Ikl (1+a) exp f IT (t, c, k) dtl P (t, c,k) [dt. _'(t~c~k)Y Let us now describe the application of the above for solving the problem (29) together with w(O) = w(l) = 0. To this end we construct two approximate solutions U(1) and U() corresponding to y = 0 and y0 = 1 respectively m m so that two solutions of (29) are given by wl(y) = exp [k fY (U(1) + )) dt], w(y) = exp [ k fY(U() +()) dt ] m2 rn Then clearly W(y) = w1(y) - w2(y) is also a solution of (29) which satisfies W(O) = 0. Upon setting W(1) = 0 we obtain the equation S I[ (U1) (U) + (1) E(2) dy Z n r i + + f [(U~ - U ) +- ~ -E dy n +1, +2,. 0 m mr k- 2,. In order to solve the eigenvalue problem we would replace E and (2) E by zero in this equation and thus find c = c(k). We would then check to see that the conditions of Theorems (II. 1) and (II. 2) are satisfied, and if so, obtain a bound on the error of the approximate eigenvalue. The procedure for obtaining such a bound by use of the above theorems is de s cribed in Section III.

-47 III. Error Analysis and Examples We shall now study the relationship between the error in an approximate solution and the error in an approximate eigenvalue and we shall illustrate the preceeding theory with a number of examples. III. 1 Error Estimates for the Secular Equation Our method involves the application of either Theorem(III. 1) or Theorem(III. 1'). Theorem (III. 1) is Roche's theorem, which we state in the following form. Theorem (III. 1): Let F(c) and f(c) be analytic functions of c in D = {c: Ic -c < ro0 r > 0 }, and let f(c ) = 0. Let m(r0) = inf CJ(c) l > 0, a = a, Ic-c I=r I a 0 and let 6 = sup IF(c) - e(c) l < m(r0). Then F(ct) = 0 for some ct E D, and so tEDot Ict c I < r ct a = 0 Theorem (III. 1'): Let F(c) and b(c) be differentiable and real on an interval I. Let F(ct) = 0, f(c ) = 0 where ct and c are on I. If m = inf I'(c) j A 0 then Ict - c < 6/m, where 6 = sup IF(c) - ~(c) I cEI c EI Theorem (III. 1') was discovered independently by Mr. H. W. Hethcote. Its proof follows by an application of the Mean Value Theorem. Let us return to the equations (2, 3) and let us assume that l =0, y2 =1 and that q(y, c,.) is an analytic function of c. It then follows that a solution of (2) is also an analytic function of c. The theorems which follow are thus tailored to Theorem (III. 1). Similar statements of theorems tailored to Theorem (III. 1') can obviously also be made, but these are omitted. Theorem (III. 2): Let a (y, c, k) be an approximate solution of (2) such that a (0, c,k) = 0, p (1, c,k) = 0 and such that pa(1, c,k) is ananalytic function of c for |c - c I < r. Let q(y,c,k) be analytic in Ic - c < r

-48and let there exist a solution 0 (y, c, k) of (2) such that JI (y, c, k) - a (y, c, k) I < 6 where 0 < y <, c - ca < r If inf ( a(, c, k) Ic - c I = r} > 6 then there exists a function O(y, ct, k) which satisfies (2) and (3), and for which Ica - C t < r This theorem applies directly to the Sinh calculus method and also to Olver's method. The DBA method also provides a simple criterion on the error of an eigenvalue. In the notation of Theorem 3.1 we can take F(O) given by the right of (13) and P(c) = d (c,...). A bound 6 on F(c) - <(c) is given by (17). If c = c exists such that d (c,...) = O, such that d (c,.. ) is an analytic function of c in Ic - c I < rO and if m(rO) = inf d (c,...).> 6 c-c-C I rn then the conditions of Theorem (III. 1) are satisfied. Hence we have Theorem (III. 3): For the DBA method, let d (c,... ) be an analytic function n of c and for some r0 > 0 and let m(r0) = inf {d (c,,..) Ic - c I = r } 0 0 n n 0 If the bound on the right of (17) does not exceed m(r0) then there exists a solution {(y, ct, k) of (2) which satisfies (3) and for which Ic - c I < r n t = 0 Theorem (III. 3) can be applied in the Sinh calculus method and also to Olver's procedure. We emphasize that the existence of an approximate eigenvalue which satisfies the conditions of Theorem (III. 2) or (III. 3) simultaneously illustrates the existence of a solution to (2) as well as yields a bound on the error of

-49 - the approximate eigenvalue. Thus, once we have a crude approximate solution which satisfies the conditions of either of the above corollaries we can proceed with assurance to determine a more accurate solution and a fortiori to determine c as accurately as we please. Let us next illustrate the application of the Riccati differential equation method to the equation (2). To this end, we list the first few approximations applicable to this equation in Table I. u r m m m 1 u = +1 q/k 0 2 uO = +1 q/k 3 u (1 + q/Zk2) -q' /2k2 - q2 /4k 4 (I + /2k2)' -(-q"+ q )/4k + U0qq'/4k4 - q'Z/16k Table I It would have been simpler to carry out the analysis of Section II. 2 without the presence of a first derivative term in (29). We shall now take advantage of our allowance of the presence of this term. We observe that the transformation (48) 4 exp I F(y) dy

-50does not alter the solution of the eigenvalue problem (2), (3) provided that F is continuous on the interval [y, Y2] =0, 1] This fact will now be used to considerably sharpen the error bounds for the approximate solution of the Riccati equation corresponding to the secular equation. The procedure is to choose an F in (48) which elimates a suitable term of a particular r in Table I, or replaces that term by a term of higher order. Consider for example the case m = 4. The transformation (48) with F = f3/k transforms (2) into the differential equation f'i V I + 2k(f3/k3),' + kZ(l+ + 3 + -)3 =O k k Upon transforming to the Riccati equation we obtain the approximation (49) Uo(1 + q/2k ) - (q' + f q ~~4 0 ~8k qq' 1 r= q q'fq 4 4 5 4k 16k where f3 was chosen such that u4 = O, i.e., f3 = 8 [-q' + f q ] Starting with (47) with (1) = () we obtain an approximate solution c. Hence we set a (c) = [U(1) (y, c) -U (yc) dy - 03 m k

-51and define F(c) by F(c) ='(c) + ( (1)- 2)(1) (l)) (1) where ql (y) is defined by (44), y =; 0 corresponding to y0 = 0 and () to yO = 1. Hence if 6() and 6(2) denote the bounds of ( )(1) and l(2) (1) as given in Theorem (II. 2) and if these bounds are uniformly valid in the region D = {c: Ic - c I < r } then we have Theorem (III. 4): Let 1(c ) = 0, where D(c) is defined above. If for a some rO > 0 we have m(r) = inf {ID(c): Ic - c I = r} > 6() + 6(2) 0 0 a 0 then there exists a function O(y, ct, k) which satisfies (3) and for which Ic -c < r. a t = 0 III. 2 Examples 1. Examples of the "Sinh" calculus Let P3(d/ds) = q0 + =q - qd/dp + q2d /dp and M = 1. Then it follows from (34) that (50) Pl = (1/k) qo[l/(p-k) + l/(p+k) -(1/2k) l/(p-k)- (l/2k) (l/p+k)] l/(p-k)3 (p+k) = (1/2k) (l/(p-k) (l/p-k)) - (1/2k) (l/(p-k)) + (1/2k/2k)3[1/(p-k) -l/(p+k)] l/(p-k) (p+k)2 =(1/2k) [l/(p-k) - l/(p+k)] - (1/2k) (l/p+k) l/(p-k) (p+k) = (1/Zk) [l/(p-k) - l/(p+k)] - (1/Zk) [l/(p+k) -(1/Zk) (l/(p+k) ]

-52Using these in (50), taking the inverse transform, setting y = 1 it follows that upon collection of terms that (50), with y O, y I becomes Y, =O0 kp( l ) + sinh = O0 = 1 Zq0 q1 1 1 2 q sinh k(- q - + ) + coshk [2q + 1)+ ]=0 k + c os 0 + 3qZ Following the suggestion of the last sentence of section I one could more easily obtain the equivalent of (50) by writing P1 = 2/(p -k ) [q0/(p -k ) + qld[l/(p-k) (p+k)]/ds + d [l/(p-k) (p+k)]/dp ] = q/(p-k)2(p+k)2 + q [1/(p-k)3 (p+k)2 + l/(p-k) (p+k) ] + 2q [l/(p-k)4 (p+k)4 (p+k)2 + l/(p-k)3 (p+k)3+ l/(p-k)2(p-k)4] The inversion could then be accomplished either by (40), (41) or by successive uses of (39). 2. An example in baroclinic instability Consider a case of baroclinic instability in a quasi-geostraphic model studied by Derome and Wiin-Nielsen (2 ). The differential equation for this model can be written as

-532 w" - [2/(x + c-1)] w' - k w = 0 which is to be solved subject to w(O) = w(l) = O The transformation w = (x + c-1)0 reduces the above equation to 2 z (51) H" -k = /(x + c-) without affecting of the boundary conditions. Upon combining Olver's method with the DBA method, or directly using the procedure of Section (II. 1), (28) with m = 2Z yields the approximation (52) c = 1/2 + 1/2[1 - 4(coth k)/k] /2 al, 2 This turns out to be surprisingly accurate; the error in the approximation -2 being O(k ) as k -> oo. A closer examination of the method of Section II. 1 yields an explanation for this. If in Olver's method described in Section II. 1 we choose the indefinite integral form for the coefficients rather than the definite integral form we get

-541 2 -A a0 A (arbitrary) a1 2 2 x c-1 (X + C - 1) 1 2 b 1, b1 dx (x+c 2 -1 ) with this choice of al and bl we get x A I a2 = - q f q dt] A 4 2 (-2) 0 4 (x + - 13) (x + - 1) (x + c - 1) and similarly for b' In view of Theorem II. 1' Olver's method thus yields the exact solution to (51) above: ky 1 -kx 1 - Ae y [1- + 1 + e [1 + 1 x c-1 x + cA and c can be determined so that P satisfies the condition The above result brings up the following question: When can we directly apply Olver's method to the boundary value problem (2, 3) so that in computing only a0 and al, c =o(l) as(k -> co? Olver's procedure (Theorem II. 1) provides as answer in the case when 1 I [q(t,c)dtl = o(k), namely that " Here Ca is the approximate value of c obtained in using only a0 and al in equation (28).

-55(53) q(x,k) - ( q(tc)dt = 2F(x,k) Upon setting x (54) f q(t,c) dt = -2u'/u the equation (53) becomes (55) u" + F(x,k)u =. From this, we obtain Theorem (III. 5): Let F(x, k) be any absolutely continuous increasing or decreasing function of x defined on [0, 1], such that IF(1, k) - F(O, k) I = o(k2) as k -> o. Let u be any solution of (55) which does not vanish on [0,1]. Then the function q defined by (54) satisfies (53). For example, if we take F -0, (55) yields u = clx + c2Z, where cl and cZ are arbitrary constants, so that q = 2/(x + C3) where c3 is an arbitrary constant. 3. The elementary linear eigenvalue problem We consider here finding approximations to the eigenvalues of the linear problem

-56(56) + k - = f(x)I "I4(0) =i(1) = 0 where for sake of simplicity we shall assume that f(x) is independent of k. We shall use the method of Section (II. 2) based on the Riccati equation. Our bounds on the error of the approximation are the sharpest known to us. It should be noted that the Liouville-Green transfromation y x = k' f 2+q7y7 ck)dY = [-(k + q(y, c, k)]-/4If may be used to transform (2) into an equation of the form (56) above, where if k- + q(y, c, k) is a twice differentiable function of y that does not vanish in [Y1 y2] then F(x) =(k2 + q) (k +q) i often dy bounded, slowly varying function of x and k. Thus our assumption that f(x) be independent of k is not unduly unrealistic. Solution: By use of Table I we obtain ik[ ) U( )] 2ik, r1 f Hence we obtain the approximation k. Let us checik Hence we obtain the approximation k = nlTr. Let us check the error.

-57By equation (41), (57)sin) x (57) PM, l(x) - 3 i f(0) +f Jf'(t) dt] 0 < < 1 k 0 * | sin k(l - f(l) I + P, (x ) 3 l[f(1)]J+f If'(t) ldt] 0 < l< I k x Since IT(X, c, k) I = IuO - Ull 0, we have by Theorem II. 2 * 1 P1 kf p 1(x) dx < - l[f(O) +f If'(t) dt] i I k o0 Z k s p, 2(x)dx < -[If(1) + If'(t) dt] - P I k o0 Let us set P = man (3{, P2) and assume that p < 1/4. We then compute. =., -1 - /( 1 1) _2Ti 2cyhtP This can always be achieved for!kI sufficiently large

-58Since moreover 1 Pm l(x) + Pmz(x ) <- [f(O) + If(l) + j f'(t) dt] k 0 we have the following result: There exist a solution U1(x) U2(x) y = e - e of (56) where, for all k satisfying p < 1/4 1 UIl(1) - UZ(1- Zikf < (1 + c)[If(O) + If(l) + j If'(t) I dt] /k 0 Clearly, U!(1) and U2(1) are analytic functions of k in Ik - n7r < 6 for all 6 > 0. Let 6 > 0 be given, and let k0 be sufficiently large so that 1 (58) 6 > (1 + a)[f(O)I + If(l) I + fI f'(t) dt]/(k0 - 6) p(, k) 0 Then it follows by Theorem (III. 4) that if nwr > k0 + 6, there exists and eigenvalue k of (56) such that Ik - nrl < p(6,k0)

-59 - -2 Note that p(6, nTrr) = O((nr). The best bounds obtainable by Olver's method satisfy p(6, nrr) O((nr) ), and our result is the sharpest known to us. We next illustrate the use of the equations for obtaining approximation on eigenvalues of (56) that are accurate to O(l/(ntr). Upon substituting the first of (49) into (47) and setting E () = = 0 we obtain the approximate equation 2 - 2F/2k2 = 2n or k" _ (nTr)k - F/2 = 0 where I F = f f(t)dt Solving the above, we get the approxiamtion k = [1 + \mz-]] 2 1 - 2F/(nir) for large n Now by Section (II. 2) we have pB(y, c, k) = 2ik, 3 2u0f 2 T(y, c, k) = Zk[f /k + U4 - u] = k - f'/k. Hence,if r r4(t) is given by the second of (49), we may take

-60 - x p4,1(x) [1r4(0)I + fJ (t) dt]/k 0 P4,2(x) [ r4(1) +f r4(t)dt ]/k 0 and 1 = k expf 2 T(t,c,k)dtl f p4 1(t) dt I I <k et f LT(t, C, k)dt| p4 1(1) - p3.2 = k exp f 2T(t, c, k)dt S p4 Z(t)dt I I < k exp f I 2T(t, c, k)dt| p4 2(0) - P Hence with p' =man (P' ) < 1/4, we have a - ( - - 1) - ~ / I Z z2 Moreover, if k > k =0 (59) 6 > (1+ a) [p () 1(1) * P4 z(1)] exp f JT(t,c,k)dt - p(6, ko) I and if 2 [1 + J F ] > k + 6 there is an eigenvalue k of 21 - 2F/(nrr) = 0 (56) such that

-61I k Z [ +1- 2]1 < p(6, k) 2 1 - 2F/(n~rr) Note here, that p(6,k) = O(k) = O[(nTr) ]. Summing up we have Theorem (III. 6): The boundary value problem (5.6) has solutions -I =k,x) where (a) Ik- nrrJ < p(6, k0) with p(6, k0) given by (58), or (b) k- /n [1+ < ] < p(6 k0) 2 1 -2F/(nrr) 1 where F = f f(t)dt and p(6, k0) is given by (59). 0 4. Exact solution of the Riccati approximation It is possible to deduce results for the Riccati equation which are analogous to the results of Theorem (III. 5). For example, in solving (31) with go + q/k2,0, f g 0 if s > 0 we obtain v = + NJ-2 1 + +q/k U = kv [1 - v ] d v' v' 2] r2 k [dxh( v 2v ( v) ] In view of Theorem (III. 5) we now ask "under what conditions is

-621 [frz(O)/2v(0) + Ir2(1)/2v(l)I + J IJd(rZ(t)/2v(t)/dtIdt =o(k2) 0 as k -> co?" The simplest case occurs when r2(t) = o('k), that is, 1 2 2) u1 + u1 = o(k ). Proceeding as in the proof of Theorem (III. 5) we easily establish the following Theorem (III. 7): Let F(x, k) be any increasing or decreasing function of x defined on [0, 1] such that IF(O, k) I = o(k ), IF(1, k) I = o(k), and let w be any solution of w"' + F(x,k)w = 0 which does not vanish on [0,1]. If v = c2w where c2 is an arbitrary constant, then the equation above is satisfied. For example, if we set F(x, k) = 0 we find that U2 given above is the exact solution of (31) whenever q(y, c,.) = k[(cx + c2) - 1] where c and c2 are arbitrary constants. 5. Connection with oscillation theory A problem related to the one of our present paper is whether the equation w" + p(t)w = 0, where p(t) is real,is oscillatory or disconjugate on [T, co). The equation is oscillatory if all of its non-trivial solutions vanish infinitely often on

-63[T, co), and it is disconjugate if each of its non-trivial solutions has at most one zero on [T, oo). An excellent treatment of this problem and a good summary of current known pertaining results is given by Willett in (22). Willett appears to be the first to give workable criteria that are simpler than solving the above equation and then testing whether the above equation is or is not disconjugate on [T, co) We relate some of Willett's criteria to our equation (2) with (3) replaced by O(0) = 0(1) = 0. Upon making the transformation y = l/t, = /t, using Willett's results and setting t = l/y we obtain the following Theorem (III. 8): Let S be the set of all numbers c such that 2 k + q(y, c, ) is real; let P(y) be defined by Y 2 2 P(y) = -fJ (k + q(t, c, ))t dt, 0 and let S' c S be such that any one of the following conditions holds for c in S' and O0< y < 1 y 2 -2 (i) f P (s) s ds < P(s)/4, or 0 y~2 2 2 -2 2 -2 (ii) 3 (f P (t)t dt) s ds < (1 - E) P (s)s ds/4 for some c > 0 00 0 then S' contains no eigenvalues of (2, 3). The proof of this theorem is,in part,based on the following connection between the above equation in w and the equivalent Riccati equation

-64w' + p(t) +w = O. This Riccati equation has a continuous solution on [T, oo) if and only if the above second order equation in w has a disconjugate solution. There exists a wide range of possibilities between an equation that is disconjugate and one that is oscillatory on an interval. From a practical point of view it would be desirable to know workable conditions which guarantee the existence of solutions of (2, 3) that have a finite number of oscillations on [Yl' Y2 ] Let us illustrate the application of the above theorem with an example. It is possible to show via the FTL theory(22)that the equation 2 2," -= [k2 + 2/(y - c)], p(0) =O(1) = 0 does not have any eigenvalues if Im c 0, k > 0. We now show that this equation has no eigenvalues also when c is real and k > 0 Let us first consider the case k = O. In this instance we need clearly only consider the case c > 1. Substituting into the above definitions we have Y 2 oo2n+l n P(y) = 2f (1 - t /c)dt = 2 Y y /[c (2n+l)] 0 n=l

-65so that after squaring we obtain Y 0oo0 2n+2 n- 1 P (s) s ds = 4 - 1/[2m+l)(2n - 2m+l)] 0 n=2 c n(2n+l) m=l Upon examining the last two identities we find that from the point of view of comparing coefficients of equal powers of y, the condition (i) of Theorem (III. 8) is not strong enough to yield disconjugacy of our equation in 9. Upon squaring the last of the above identities, multiplying through -2 by y dy integrating from 0 to y and simplifying we obtain y t 2 2 2 2 00 2n+l n=2 (60) f (f P (s) s ds) t dt = y _0 0 n=4 c (2n+l) m=2 (2m+l) (2n -2m+l) m-l 1 n-m+l 4 1 4 1 m +l n-m+l 1 2j+1 We are now in position to apply the condition (2) above. This will be satisfied if it can be shown that n- 1 0 < (1 - s) 1/I[(2m+1) (2n - 2m+l)] m=l n- 2 m-l n-m+l 4 1 4 z {1/[(2m+l)(2n - 2m+l)]( 2i+ n-m+l 1 m=2 i=l j=l1

-66for n = 4, 5, 6,..., for some s > 0. This last inequality is certainly satisfied if m-l 4 1 m+l. 2i+1 = i=l for m = 2, 3, 4,... Now 4 1 2 2C 2i+1 < +In (2m-1) < 1 - ~ m+1 2i - m+1 i=l where c =1 - in 9/3 Hence the above equation in ~ is disconjugate on [0, 1] for k = 0 Examination of the Neumann series for a solution of this equation shows that it is disconjugate also for k > 0.

-67IV. Numerical Methods A host of numerical methods exists for obtaining discrete solutions of differential equations. Among these, a particularly powerful 3-step method suitable for solving the equation (2) is the Numerov method (see e. g. (24)). Applied to the equation y" = g(x) y + f(x) the Numerov method is 2 1 2 1 2 (61) 62 [(1 h2 g)y hY2 + f + 1 62ff 1 1 where g = g(xn) = g(xO + nh), 6g(x) = g(x + h) - g(x - h) and where h > 0 is the step size. Expanding (61) we compute Yn+l n = 1, 2, 3... from 12 52 12 2 1 2 (1 - h gn+l) =(2 + 6h )h gnn 2h gn)Yn_ + h (f+ fn) Tthe error in using this equation is approximately equal to -/ 240 6 at the ntth step. Let us assume that g(x) = g(x, c), f- 0, and that we want to determine an"eigenvalue" c such that y(0) = y(l) = O.

-68Differentiating the equations (62) y" = g(x, c)y, y = y(O) = 0, Y1 = y(h) = Ah with respect to c, we get the derived variational equation (63) y + g(x, C)y + gc(x, c)y Y() =Yc(h) = 0 An algorithm is to set Nh =1, start with an approximate value of c and compute y for this value c of c, n = 2, 3,..., N by use of the Numerov method on (62). We next also compute y, n = 2, 3,..., N using the Numerov method on (63), with the yn given by the numerical solution of (62). We then set YN C =C - a c, N to get an improved value of c. In the case of convergence, E = Ic - C I bounds the error in c, since ~ is replaced roughly by E2 after each new iteration. In the above we have combined two methods, Numerov's and Newton's method. With fixed h, Newton's method will converge to an exact eigenvalue of the discrete solution. We can test the accuracy

-69 - of the numerical method by repeating the computations with a finer mesh size h. Thus ini tlte case wheern this method produces an approximate eigenvalue we can get at the error by a posteriori means. Unfortunately Newton's method converges only if c is sufficiently close to a true eigenvalue of the discrete solution. Thus instead of "seeking in an attempt to find", it may often be simpler to obtain an approximate eigenvalue by use of one of the above analytic methods. In the case when the first derivative is not absent in a second order equation we could of course always remove it by a transformation of the form (49). However this sometimes complicates the resulting equation to the extent that it is much more efficient to apply a less accurate method directly and to combine it with Newton's method. We illustrate with the following example taken from [25]. d2 IdU-c c (64) d _ 1 y = 0 d2 U-c U-c-c c dx U-c dx r y(O) = y(l) = 0 where y is related to the vertical wind velocity, A is related to the static stability and U = U(x) is the zonal wind. The problem is to determine c given U and c. The derived variational equation is r

-70d2 1 1 dy U - c - c Yd __ dx gtC - A(1 r ) + ]Ul_ A ( d2 [U c+U -c -c Ud-c c dx r c 1 2 1 2 dy 2 r U - - c c dx AC) Taking U = Px(l - x) we approximate 2 by Y 2Y + Y1) dx dy/dx by (Y>+1 - Y1)/2h, y(h) = 0,?(h) =5h. This results is an error of order h4 at each step. We combine this approximate method with Newton's method to simultaneously solve both of the above equations. With P = 10, a = 1, c = 2.14, h = 1/36 we obtain r c =.43030, - 1. 56988 with A =0 c =.36213, - 1.57459 with A =1 c =. 37474, - 1.61654 with A = 2 The number c can be obtained another way when A = 0 since the above differential equations can be explicitly solved then. The exact values of c when A = 0 are given by c =. 42904, c = - 1. 56904 which indicates that the above results are correct to 3 significant figures.

REFERENCES 1. Drazin, P. G., Howard, L. N., "Hydrodynamic stability of parallel flow of inviscid fluid. " Adv. in App. Mech. 9, 1966, p. 1-89 Academic Press, New York. 2. Derome, J. F., and Wiin-Nielsen, A., "On the baroclinic instability of zonal flows in simple model atmosphers." Tech. Rept. No. 2, Department of Meteorology and Oceanography, Univ. of Mich., 1966. 3. Drukarev, G. F., "The theory of collisions of electrons and atoms." Soviet Physics, JETP 4, 1957, pp. 309-320. 4. Dryak, H., "Determinant solutions of the Fredholm equation with Green's function kernel." J. Math. Phys. 4, 1963, pp. 1536-38. 5. Aalto, S. A., "Reduction of Fredholm integral equations with Green's function kernels to Volterra equations." Tech. Rept. No. 26, Department of Mathematics, Oregon State University, Corvallis, Oregon 6., "Solutions of Fredholm equations with Green's matrix kernels." Math. Res. Center (USA). Tech. Summary Rept. No. 714, 1966. 7. Tricomi, F., "Integral Equations". Interscience, New York, 1957, vi + 238 pp. 8. Henrici, P., "Discrete variable methods in ordinary differential equations." Wiley, New York, 1962, + 407 p. 9. St. Matthew, Chapter 7, Verse 7, King James Version. 10. Derome, J. F., and Dolph, C. L., "Three dimensional disturbances in the non-geostrophic Eady model of the atmosphere. " 32 pages. Submitted to J. Atmos. Sci. 11. Coddington, E. A., and Levinsion, N., "Theory of ordinary differential equations." McGraw-Hill, New York, 1955, vii + 429. t2. Arnason, G., "The stability of non-geostrophic perturbations in a baroclinic zonal flow." Tellus, 15, 1963, pp. 205-29. -71

-7213. Haltiner, G. J., "Finite difference approximation to the determination of dynamic instability. " Tellus, 15, 1963, pp. 230-40. 14. Frank, P. and v. Mises, R. V., "Die differential und Integralgleichungen der Mechanic und Physik," Vol. I, M. Rosenberg, New York, vii + 919, 1943. (see pp. 353-53). 15. Courant, R., and Hilbert, D., "Methods of Mathematical Physics," Vol. I., Interscience, New York, xv + 561 pp., 1955. 16. Widder, D., "Advanced Calculus", Prentice Hall, Englewood Cliffs, New Jersey, xvi + 520 pp., 1961. 17. Epstein, B., "Partial differential equations", McGraw-Hill, New York, x + 273 pp., 1962. 18. Olver, F. W. J., "Error bounds for the Liouville-Green (or WKB) approximation," Proc. Camb. Phil. Soc., 57, 1961, pp. 790-810. 19. Stenger, F., "Error bounds for asymptotic solutions of differential equations," J. of Res. of NBS., B., vol. 70 B., 1966, pp. 167-186. 20. Erdelyi, A., "The integral equations of approximation theory. " Contained in "Asymptotic solutions of differential equations and their applications," Wiley, N. Y., 196, Ed. by C. H. Wilcox. See pp. 211-229. 21. Cesari, L., "Asymptotic behavior and stability problems in ordinary differential equations, Academic Press, New York, 1963, vii + 271 pp. 22. Willett, D., "On the oscillatory behavior of the solutions of second order linear differential equations," Preprint from Math. Dept. The Univ. of Utah, Salt Lake, Utah, 84112. 23. Lin, C. C., "The theory of hydrodynamic stability", Cambridge University Press, 1955, ix + 155 pp. 24. Hartree, D. R., "Numerical Analysis", Oxford 1958, xvi + 302 pp. 25. Wiin-Nielsen,,A., "On baroclinic instability as a function of the vertical profile of the zonal wind", Monthly Weather Review 1967, pp. 733-739. UNIVERSITY OF MICHIGAN 1111111111111111111111111111111 111 1111 I 3 9015 02082 8151

AC KNOWLEDGEMENTS The authors would like to thank their colleagues Professors K. Case, S. Jacobs, T. Storer, J. Ullman, A. Taylor and C. Yih for their interest and constructive suggestions.