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.