ABSTRACT NONLINEAR PERTURBATIONS OF SPECIAL-FUNCTION OPERATORS BY Larry Charles Andrews A systematic perturbation theory is presented for the analysis of nonlinear forms of boundary-value problems. In particular, those equations are considered whose unperturbed form belongs to the class of linear special-function equations. The nonlinear terms are then regarded as perturbations of the special-function operators. A matrix representation is obtained for the perturbed operator utilizing the coordinate (Schrodinger) representation of quantum mechanics, trunca- tion and diagonalization of which will determine the perturbed eigen- values and eigenvectors. In most cases the method is applicable even when the perturbation term is of the same order of magnitude as the remaining terms, or perhaps even larger. To illustrate this point the nonlinear Legendre-like equation (d/dx)(l - x2)(du/dx) + Au + ox2u2==0 is solved for those cases when a = l and a = 5. Other examples include the Hartree equations for the helium atom where a qualitative compar- ison of the ground state energy is made with experimental data, and a detailed analysis of the van der Pol equation for a = 0.5 and a = l. NONLINEAR PERTURBATIONS OF SPECIAL-FUNCTION OPERATORS By Larry Charles Andrews A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Metallurgy, Mechanics and Materials Science 1970 ,4- // f 4 / 5 _ \T'J’ -' U“? M. -~2 M‘ :5 ’70 ACKNOWLEDGEMENTS I should like very much to thank my thesis director, Dr. T. Triffet, for his encouragement and many helpful suggestions provided throughout this thesis. In addition, may I express my gratitude to the Magnavox Company for their cooperation, and to Mrs. Jean Emig for her excellence in typing. ii TABLE 1. INTRODUCTION . . . . . . . 2. GENERAL DEVELOPMENT . . . Formulation . . . . . The Operator B . . . The Operators F and A . The Operator Z . . . an> 0. OF CONTENTS 3. A NONLINEAR LEGENDRE-LIKE EQUATION . Formulation . . . . . . The Operator B . . . . Matrix Representation . Matrix Representation . Summary of Results . mcnw> for A . for Y . 4. GENERAL NONLINEAR SPECIAL-FUNCTION OPERATORS . A. A Nonlinear Hypergeometric-like Equation B, Other Forms of Nonlinearity . . . . . . . 5. OTHER EXAMPLES . . . . . . A. The Hartree Equations for the B. The van der Pol Equation . . 6. GENERAL REMARKS . . . . . A. Higher-Order Perturbations . B. Analytical Methods . BIBLIOGRAPHY . . . . . . . . . APPENDIX Computer Program . . . . iii Helium Atom Page 10 10 11 20 22 26 26 27 31 33 35 4O 40 44 49 49 SS 62 62 65 69 71 Table 1° 2. LIST OF TABLES Page A comparison of eigenvalues for Equation (3.1) obtained by standard and matrix perturbation methods with a = l . 37 Eigenvalues for the matrix given in Equation (3.58) with a = 5 O O 0 O I O O O 0 O O O O 0 O 6 O O O O O O O I O 0 39 Approximations of u(t) and u'(t) for the van der Pol equation obtained by numerical and matrix techniques . . 60 iv 1. INTRODUCTION From the theory of representations it is known that a linear operator L that transforms a Hilbert space into itself gives rise to a matrix representation of that operator.1 The matrix elements can be defined by the inner product, ij = (¢j. L CAR). Cl-l) where {¢k} represents a complete set of orthonormal vectors over the space. We shall extend such representations to certain nonlinear forms Of special-function Operators, the matrix elements of which can be computed by algebraic manipulations instead of utilizing ex- pressions like (1.1). The type of nonlinear Operator to be considered is represented by the hypergeometric-like form 8 = .9. (1 - x2) _51_ - 2(u + vx) _gl_+ax2uk, (1.2) dx dx dx from which we obtain the nonlinear eigenequation (S + Ak)uk = O. (1.3) The term axzuk is regarded as a perturbation of the linear hyper- geometric Operator L=_£l.(1-x2)_§_-2(u+vx)_‘_l.. (1.4) dx dx dx Thus a will serve as the perturbation parameter, but will be treated as a positive constant of arbitrary magnitude in contrast to the usual restriction that<1 remain small.2 The Operator S can now be written S = L + N, (1.5) where N =(1x2uk. Note that we can further represent N as the product of a linear perturbationax2 with a nonlinear perturbation uk. For a = 0, Equation (1.3) becomes a hypergeometric equation with known eigenvalues and eigenvectors.3 We shall seek an eigenvalue-eigen- vector solution for Equation (1.3) with a f 0 that will reduce to the unperturbed solution when we again let a = O. The perturbation method is of fundamental importance in quantum theory. The mathematical description of noninteracting particles is governed by the (nonrelativistic) SchrOdinger equation or by the (relativistic) Klein-Gordon equation.4 The separation of variables technique applied to these equations often leads to extended forms of special-function equations. Because direct solution of these equations is usually not possible, approximation methods such as perturbation techniques have proved to be most valuable. In general, though, any basic field equation of the elementary particles cannot be linear since it must define interactions between the particles. Heisenberg suggests that the mathematical form Of such an equation should be given by i av §x_.+ 22 0V: x(x* 0v x): = 0, (1.6) BXV where x(x) is a local field Operator, the 0V are the conventional Pauli matrices, and 2 is an arbitrary constant with the dimension of S . a length. (The dots :: in the second term refer to the definition 3 of a product of three field operators at the same space-time point.) The anharmonic oscillator defined by the equation of motion O + aQ3 = 0, where each dot designates a time derivative, is a simpli- fied representative of (1.6). When equations such as that for the anharmonic oscillator define a conservative holonomic system with time-independent constraints which possesses a kinetic energy that is homogeneous-quadratic, then the system may equally well be described by a linear SchrOdinger equation through an application of the Hamilton equations Of motion. De Broglie has also suggested the existence of a nonlinear wave equation like (1.6) in accordance with his double solution theory.6 Thus it seems natural to consider perturbation methods for approximating solutions Of these nonlinear equations. However, the nonlinear terms are Ordinarily of the same order of magnitude as the linear terms, so that conventional perturbation theory is not appropriate. It is for this reason, and because of the possible area of application cited, that we present our perturbation method from a quantum mechanical point of view, although this is not essential. It is most convenient to write differential Operators in terms of dynamical variables (observables such as position, momentum, and energy), but these variables must be treated as operators which are generally noncommutative.7 The commutations of such operators are proportional to the classical Poisson brackets, i.e., [Q: R] = QR " RQ = j-{QaR} 3 (1°7) where Q and R represent two quantum mechanical operators and { , } denotes the Poisson brackets.8 4 We wish to write the differential equation (1.3) entirely in terms of dynamic variables. This involves the transformation of the functions uk(x) of the differential equation to the vectors w(k)(q) of the Operator equation. Such a transformation can be represented by9 _. ' (k) ”k(X) = w(0) e1px 1 (q). (1 8) where x is treated as an unrestricted coordinate variable and where p and q satisfy qp - pq = i (with‘fi = 1). In addition to (1.8) we shall also find it useful to establish the two relations x URCX) = TAO) eipx qw(k). (1 9) d (0) - i a.)_(.uk(x) = $' eipx pw(k). (1.10) To verify (1.9) we first note that since x e1px = eipx q - q eipx, (1.11) it follows that E{O) x eiPx W(k) x uk(x) _ ETD) eipx qw(k) (1.12) 1 where it is assumed that qwco) = 0. We verify (1.10) merely by differentiating (1.8). Thus, using (1.8), (1.9), and (1.10) it is easily deduced that (k) Tm) e11px pa - q21p w , (1.13) d _ 2 a§'(1 x ) 9-0- >1: 201 + vx) $13 = "(7(0) eiPx 2101 + vq1p W), (1.14) X 1kuk = ITO) eiPX 1(k) ((k). (1.15) The nonlinear term in (1.3) requires further analysis in its trans- formation to Operational form. We write N l (XX u -<1x2 $10) eiPx ¢(k) ETC) eiPx (Ck) (k) -(0) ' W 1px e aq2 W(q,w ) w(k) (1.16) where W(q, wck)) = u(k) 6(0) eipx is in general some operator function of q and u(k). For most of the applications we have in mind, however, (k) we can express this operator as V(q) w where, as will be shown in Section 2, V(q) is automatically absorbed in the process of forming the matrix elements of the nonlinear operator. Hence, (1.16) becomes _. . .. 2 axzui = w(0) e1px <1q2 V(q) w(k) ; (1.17) and combining our results we may now rewrite (1.3) in the operational form (A - 1(k3)w(k) = o, (1.18) where A -- pa - q2)p + 2101 + vq)p - mq2 V(q) 10‘). (1.19) In general, the Operator A may be taken to have the form A = I‘(P, q)B(p. q) + Mp. q)ch. q. W”). (1.20) where B represents an unperturbed self-adjoint linear operator, T and A represent non-singular linear perturbations, and Z is a non- linear perturbation. For the above case B = p(1 - q2)p + ZiCu + vqlp. r = 1, A = - aqz, 2=meWL (Lu) Thw30ther examples of differential equations with their Oper- ational forms are: 2 (l) §_%k + (Ak - x2)uk + axuufi = 0, x _ 2 [p2 + q2 - an ch) w(k) - A‘k)] w(k)= o. B = P2 + qz. r = 1, A = - aq“, Z = V(q) v(k)2. _ 2 41_ _ 2 du _ 2 £23 = (2) (l ax ) dx (1 x ) E;k + 2a(l x )xuk dx +Akuk 0, [(1 - aq2)p(1 - qzlp - ia(1 - qzlqp VIq) w(k) - A(k)]w(k)==0, B = p(1 -_q2)p. T = l - aqz, A = - iacl - q21qp, Z = VIq) v(k). The unperturbed form (a = O) of Equation (1.18) is given by k k (B - bc ))U( ) = o, (1.22) where bck) and u(k) denote the unperturbed eigenvalues and eigen- vectors, respectively. A factorization method for solving linear eigenequations was suggested by SchrOdinger, developed in the differ- ential equation form by Infeld and Hull,10 and given its more recent operational form by Green and Triffet.9 The method is also a power- ful tool for studying recurrence relations obeyed by the special functions.11 In the factorization method it is generally assumed that for some self-adjoint Operator M, whose eigenvalues m we wish to determine, there exist mutually adjoint (adjoint as used here refers to the complex conjugate transpose) linear Operators J+, J-, satisfying [M, J*] = MJ+ - J+M = J*, (1.23) [M, J‘] = - J”. (1.24) If, in addition, they satisfy certain other conditions we say the Operator B admits a factorization. In particular, if a transformation can be found such that B has the form B=p2+wm. man then Equation (1.22) can be factorized into the two equations O) ' . [K(q, M + 1) - 1p] um = [b‘33 _ ,(m , 1)]a U(J) m+1 (1.26) u 0 8 o [K(q. M) + 1p] ”#3) = [bcj3 - a(m)] Ufiii (1.27) with 3* = K(q, M + 1) - ip, (1.28) J- = K(q. M) + ip- (1.29) We should point out that Equations (1.26) and (1.27) were originally presented in the differential equation form m+1 , 1.30 j ( ) [K(x, m + l) 5;) uj [bj a(m 1)] u _d_ m- _ 1: m-1 [1((x, m) + dx] uj — [bj a(m)] uj . (1.31) It has been shown that a necessary and sufficient condition for factorization is that §_and a_satisfy the equation Ji.[K(x, m + 1) + ch, m)] + K2(X. m + 1) dx - x2(x, m) + a(m + 1) — a(m) = o. (1.32) Miller has also shown that Riccati equations similar to (1.32) are sufficient to determine the 4-dimensional Lie algebras£3(a, b), whose representations correspond to a study of the special functions of . ll hypergeometrlc type. The factorizations (1.26) and (1.27) now permit us to obtain the eigenvalues and eigenvectors of the Operator B in a simple and + - . elegant manner. However, the step Operators J , J are not unique, as more than one factorization is often possible. In addition, Coulson and Joseph have exploited the use of self-adjoint step 12,13 operators. Green and Triffet have recently introduced a more systematic procedure for determining the step Operators.9 They have shown that construction of a sequence of linearly independent Operators, which do not commute with the unperturbed Operator B, will ultimately lead to an apprOpriate form of the Operators. Step Operators determined in this fashion will usually not be mutually adjoint. However, they re- tain their most useful properties and are more general than those defined by the Infeld and Hull factorization technique. Also featured is the basic notational scheme utilized here and an algebraic method for finding matrix representations of linear operators which has been adapted to the present purpose. 2. GENERAL DEVELOPMENT A. Formulation A nonlinear differential equation of the form (L + 0N + Ak)uk = 0, (2.1) where L is a linear differential Operator, N is a nonlinear differ- ential operator, and a is a constant of arbitrary magnitude, can be treated as a perturbation problem. The boundary conditions will be generally based on physical considerations and referred to as "physical boundary conditions" as opposed to the more common arti- ficially imposed boundary conditions related to the Sturm-Liouville equation.14 For the special case when a = 0 Equation (2.1) becomes linear, so that we Shall regard a as the perturbation parameter. Writing (2.1) in the Operational form displaying the quantum mechanical conjugate variables p and q we arrive at (A - 100M“) = (1‘3 + AZ - 100w“) = o; (2.2) thus we make the associations PB + L and AZ + N. The linear operator L consists of a self-adjoint operator, designated by B, multiplied by a perturbation T (ordinarily unity in the examples to be considered). The nonlinear Operator N is decomposed so as to separate out the (k) dependency on the eigenvectors w (k) function of w . Introducing the distinct Operators A and Z is not ; consequently, only Z will be a necessary but is a convenience for calculating the matrix elements, since A may then be interpreted as a linear perturbation. 10 11 We seek a matrix representation for A, A. = Z (T. B Jk n jn nk + Ajnznk)’ (2'3) diagonalization of which by numerical methods will determine the perturbed eigenvalues. Hence, we find separate matrix represen- tations for each of the operators T, B, A, and 2. For each choice of a basis which defines our Hilbert space we expect a different, but equivalent, representation. The most convenient basis to select is the set of normalized eigenvectors belonging to the unperturbed Operator B; thus the representation for B is the diagonal matrix - (j) Bjk — b sjk, (2.4) where the ch) are the eigenvalues of B and 6jk = 0 for j # k while 5jk = l for j = k. We then proceed to determine the matrix repre- sentations for F, A, and 2 with respect to this basis. B. The Operator B The Operator 8 need be self-adjoint only in the sense that 3': nB*n'1 = B, where 8* denotes the Hermitian conjugate and n is a positive-definite linear operator. We shall assume that the eigen- values of B are distinct, nondegenerate and bounded below, but not above, with b(1) as the first eigenvalue. Certain properties are characteristic of such an Operator, and these we now list; proofs 15 can be found in the references. ’16 12 (1) A linear Operator of the second order in p is self-adjoint if and only if it has the form 8 = pF(q)p + Gm)- (2) For a finite domain the eigenvalues form a denumerably infinite sequence. (3) The eigenvalues are real. (4) There exists a complete system of associated eigenvectors satisfying the orthonormal relations mm. W”) = 6. Jk’ We wish to define the step operators JCI), ch) associated with B that give the recurrence formulas J(1)u(k) = uCk+1), (2.5) J(2)U(k) = u(k-l), (2.6) and obey the commutation rule [3, J(n)] = A(n)J(n), n = 1, 2 (2.7) where A(n) designates operators to be determined. A set of linearly independent Operators J1, J2, ., can be considered as a basis for some linear vector space. Hence, any arbitrary operator defined on the space can be represented as a linear combination of the base vectors, and the Operators J(n) will have the expansions (n) = (n) = J 2kagk , n 1, 2 (2.3) 13 (n) (n) k (2.8) will define the step Operators once the Jk and sin) are where the 5 are components of J in the given basis. Equation determined. Employing (2.8) we note that [3, .100] = 330‘) - 3003 (n) £k(BJk - JkB)£k . (2.9) Since the Jk form a basis for a linear vector space they must satisfy BJ — J B = X.J.c , 2.10 k R J J jk ( ) where the cjk are either numerical constants or Operators which commute with B. This is a consequence of the condition that the Operator [B, Jk] must be representable as a linear combination of the base vectors. Substituting (2.10) into (2.9) yields (n) (n) B, J = X E J, , . 2.11 Comparing (2.11) with (2.7) indicates that the Operators Cjk must comply with the relation 0‘) (n) (n) = A 2.12 (n) (n) which identifies A as the eigenvalues of the matrix cjk and g as the right-eigenvectors. Accordingly, Equation (2.11) becomes (11) = (n) (n) [B,J ] A 23.33.53. = 10‘) 3‘"). (2.13) 14 Usually there will only be two nonvanishing eigenvalues A(n), a ,0) positive one designated by and a negative one indicated by AC2). Equation (2.10) is the defining relation for the sequence of linearly independent Operators J1, J Experience indicates 2’... that J1 should be chosen as an elementary function of the coordinate q and then, forming the commutations 8J1 - J18 = chll + J2c21’ BJ2 - J28 ch12 + ch22 + J3c32, BJ - J B = Z.J.C. k k J J Jk’ (2.14) successive Operators J2,..., can be determined. In the examples presented here, where the Operator B is a special-function operator, such a sequence will always terminate. In many cases it is desirable to introduce an Operator M, (k) where B = B(M), whose eigenvalues m are separated by unity, i.e., m(k) = m(1) + k - 1. Such a device usually simplifies the algebra, as the Operator M satisfies the relations MJ(1) = J(1)(M + 1), (2.15) MJ(2) = J(2)(M - 1). (2.16) When B is a special-function Operator it will be at most a quadratic function in M, 3 = a0 + alM + a2M2, (2.17) 15 where a1 and a2 are numerical constants that may be determined from the relations (1) [a1 + a2(2M + 1)]J(1), (2.18) F—. as b C.) hi ll - [a1 + a2(2M - 1)]J(2), (2.19) and (2.7). These give > ll a1 + a2(2M + 1), (2.20) > A II - [a1 + a2(2M - 1)], (2.21) which can be solved for a1 and a2. Usually the constant a may be 0 placed equal to zero, but in other cases may be determined from the form of the basic operator B. Suppose for eXample B = p2 + qz. Selecting Jl = q we have [8. J1] = pzq - qu = - 2ip = 2J2, [B, J2] = - iqu + ipq2 = Zq = 2J1, l where J2 = - ip. The eigenvalues of the matrix cjk are A( ) = 2 and (2) . . (1) (2) A = - 2; the right-eigenvectors are E, = (l, 1) and E! = (l, -l). . (1) . (2) _ . Thus we obtain the step operators J = q - 1p and J — q + 1p. (2) Substituting A(1) and A into Equations (2.20) and (2.21) suggests that a1 = 2 and a2 = 0. Here we have a0 = 0 so that B = 2M. 16 (1) (2) Having the desired step Operators J and J , we next apply them to the calculation of the eigenvalues and eigenvectors of B. If u(l) denotes the first eigenvector of B then multiplying by J(2) must give J(2)U(1) = 0, (2.22) which follows from Equation (2.6). Multiplying by JCI) yields J(1)J(2)U(1) = 0. (2.23) (1) (2) . . . The operator J J 15 express1ble as a funct1on of B (or M), J(1)J(2) = f(B). (2-24) Therefore the first eigenvalue of B or M is determined from (2.23), which gives f(b(1)) = 0. (2.25) Successive eigenvalues of B can be determined from the relations m(k) = m(1) + k - l and B = B(M). To illustrate the above for the operator B = p2 + q2 we observe that f(3) = J(1)J(2) = q2 + iqp - ipq + 32 = B - 1 = 2M - 1. (2.26) Hence, mCl) = l/2 with m(k) = k - 1/2. It follows that the 17 eigenvalues of B are (k) b = 2m(k) 2k - 1. (2.27) We will also be interested in the matrix elements for the (1) (2) operators J and J . These can most easily be obtained by writing them as unit step Operators multiplied by their respective magnitudes. This is analogous to the prOperty that a vector in Euclidean space can always be written as a unit vector in the direction of the given vector multiplied by its magnitude. The unit step Operators will be designated by 3+ and 6 with matrix elements (8+)jk = 63' k+1, (9-)jk = 6j+1 k (2.28) Th . (1.) (2) . e magnitude of J and J can be calculated from 1 l -l |J( )I = (J‘ )J( 5%. (2.29) |J(2)l = (3(2)J(2))15, (2.30) where 3(n) denotes the adjoint of J(n), e.g., the adjoint of J = a(q) + ib(q)p is 3': a(q) - ipb(q). It can be shown that a _. . 9 . relation between J(1) and JCZ) must exist, thus we conjecture 3(1) = (2) J g(B). (2.31) Note that 3(1) and J(2) are mutually adjoint only if g(B) = 1. This distinguishes these Operators from the similar ones referred to in Section 1 where it was assumed such operators are mutually adjoint. 18 Now J(1)3{1) = f(B)g(B) = [11mm - 1)]2, (2.32) and 395(2) = f(B)/g(B) = [h(2)(M - 1)]2, (2.33) where h(1) and h(2) are introduced here for notational convenience. We may therefore express the step operators in the form 3(1) = h(1)(M — 1)e+ = 6+h(1)(M), (2.34) (2) 31(2) J a h(2)(M)e_. (2.35) (M-l) The matrix elements of J(1), J(2) are now clearly given by (1) 1 k (J )jk 6 k+1 h( )(m( )), (2.36) i (J(2)) 6 h(2)(m(j)). (2.37) jk j+1 k For the Operator B = p2 + q2 we have already determined that f(B) = 3 - 1. Also, 3(1) = q + ip = J(2) (2.38) so that g(B) = 1, thus Jcl) and J(2) are mutually adjoint in this particular case. We next resolve 19 h(1)(M - 1) = h(2)(M - 1) % (2M - 1) , (2.39) or replacing M by M + l h(1)(M) = h(2)(M) = (2M + 1)%. (2 40) Consequently, JCI) = 6+(2M + 1)%, (2.41) 3(2) = (2M + 1);2 e_, (2.42) so explicitly displaying the matrix elements we have 0 /fi' 0 o ... J(1) = , (2.43) o o 45 o o o o /8 20 J(2) = . (2.44) C. The Operators F and A The Operators F and A are not required to be self-adjoint although we will presume that F is positive-definite. The matrix representations of both operators can be handled similarly so that we will only display A in the following discussion. We denote by A(0) the diagonal matrix whose nonzero elements are identical with the center diagonal elements of A. Similarly, we denote by ACr)ef and B:A('r) the matrices whose elements are identical with those of A in the r-th diagonal above and below the center diagonal respectively, and zero elsewhere. Hence A = A(0)+ Z (9:A('r) + A(r)6f), (2'45) r=l Mr) = A“) (M). with matrix elements 21 . 6 A(‘r)6 ACT)6 , 2.46 3k 1 1k i til ( k j k+r * j j+r k) ( ) Air) = A(r3(m(k)). (1) Also we suppose that A can be represented as a function of J , J(2), and M, + E (J(1)r3('r) + ECT)J(2)T), (2.47) g(r) = 50‘) (M). From Equations (2.34) and (2.35) it follows that J(1)r = 6:h(1)(M)h(1)(M + 1) ... h(1l(M + r - 1), (2.48) flnr=flmmm9MM+n.HhQMM+r-nq, 94% so that upon comparison of (2.45) with (2.47) we identify Am) = 3(0) (M), (2.50) A(’r) = B(‘r)(M)h(1)(M)h(1)(M + 1) ... h(1)(M + r ~ 1), (2.51) a(r) = 3(r)(M)h(2)(M)h(2)(M + 1) ... h(2)(M + r - 1). (2.52) From these relations the matrix representation of appropriate pertur- bation terms can now be obtained. This is explicitly illustrated in Section 3. 22 D. The Operator Z There is an inherent difficulty present in the form of the operator Z that does not occur for the linear Operators, viz., it is a function of the unknown perturbed eigenvectors w(k). However, this difficulty can be alleviated by applying the same techniques character- (k) istic of the standard theory. It is assumed that w is a continuous function of the parameter a (when such a parameter does not explicitly appear in the eigenequation it can be introduced and later set equal to one) so that we can form the Taylor expansion about a = 0, k k k w( )(a) = ¢( ) + 40f ) + ... , (2.53) k k where wf ) is the derivative of w( ) with respect to 0, evaluated at k a = 0, and ¢( ) is a normalized eigenvector of the Operator B. Assuming that Z is analytic we may also expand it about a = 0, (k) 2(6, 0 )= Y + 021+ , (2.54) where we define (R) Y = 2(09 ¢ )1 (2'55) and 32 (k) 21 = (k) w) . (2.56) am 0:0 For first-order perturbations we shall henceforth approximate the operator Z by Y. The validity of such an approximation will, of course, depend on the convergence properties of (2.54). From (2.55) we note the dependence of Y on the index k. 23 Since this index will refer to the k-th column in the matrix represen- tation, we must calculate the matrix elements of Y by columns, i.e., for each column, Y is essentially a different operator. It is this column selection for the matrix elements of Y that constitutes the primary function of the operator V(q) defined in the Introduction. For present purposes we have absorbed this operator in the expansion (2.54). The Operator Y can be represented in a manner analogous to that for A. We shall denote by Y(O), 6:Y(-r), and YCr)ef matrices similar to those defined by the corresponding expressions containing the Oper- ator A. Thus we write Y = Y(O) + Z (0:Y('r) + Y(r)ef), (2.57) with matrix elements Also if we assume that Y, like A, can be expressed as a function of 2 J(1), J( ), and M, then Y=E(0)+ Z (J r l (1)rE(-r) + B(r)J(2)T) (2.59) 24 which, with the aid of Equations (2.48) and (2.49), allows us to make the associations 3“” = Emcm. (2.60) Y(’r) = 3("r)(M)h(1)(M)h(1)(M + 1) ... h(1)(M + r - 1), (2.61) Y(r) = B(r) (M)h(2) (3011(2) (M + 1) ... 11(2) (M + r - 1). (2.62) Because of the operator V(q), the above representation must be repeated for each column of the matrix Y, but only those elements that occur in the given column need be calculated. For example, suppose that we wish to calculate the matrix elements in the third column of Y 3 ¢( ). Utilizing the above Equations (2.57) - (3) (2.62) we would Obtain the matrix expression for ¢ , that correspond to Y = P (3) (3) (3) (3) T ¢11 12 13 ¢14 '°° (3) (3) (3) (3) 3 ¢ ¢ ¢ ¢ ... ¢( ) = 21 22 23 24 . (2.63) (3) (3) (3) (3) 31 ¢32 33 ¢34 "' (3) (3) (3) (3) $41 ¢42 ¢43 44 "' However, out of this matrix representation we only select those elements in the third column for the representation for Y. Therefore, it becomes unnecessary to calculate the remaining elements. To help eliminate other unnecessary calculations in computing the matrix elements of Y we note that whenever the coordinate q has a 2S symmetric matrix representation then E(r) (2) (2) h(2) (M)h (M)h (M + 1) ... (M + r - 1) (-r) (1) (1) h(1) = E (M)h (M)h (M + l) (M + r - 1). (2.64) This property means that we need not calculate any of the coefficients (2 r of J ) in Equation (2.59), as utilization of the coefficients of (1)r J will determine all matrix elements in a given column. 3. A NONLINEAR LEGENDRE-LIKE EQUATION A. Formulation For the special case of Equation (1.3) when u = v = 0 we Obtain the nonlinear Legendre-like equation -£L 1 - x2 ii.+ axzu + A u = 0 3.1 which in Operational form featuring p and q becomes [3(1 - q2)p - aqzVIq)w(k) - A(k)]w(k) = 0. (3.2) Expanding w(k) in a Taylor series about a = 0 yields w(k) = ¢(k) + awgk) + ... (3.3) where 6(k) again represents the normalized eigenvectors of the un- perturbed Operator. Hence we may write [p(1 - qzlp - aqzVIq)(¢(k) + awfk) + ...) - A(k)]w(k) = 0, (3 4) and, upon retention of linear terms only in a, we finally Obtain [3(1 - q2)p - aqzVIq)¢(k) - A(k)]v(k) = 0. (3 5) Thus we can identify the Operators B = 3(1 - qzlp. r = 1, A = - aqz. Y = V003“). (3.6) 26 27 Next we must find the eigenvalues and eigenvectors of the un— perturbed Operator B, then the matrix representations of the operators A and Y can be constructed. B. The Operator B The unperturbed form of Equation (3.5) is given by [p(l — q2)p - b(k)]U(k) = 0. (3.7) To determine the associated step operators J(1), J(2) it is convenient to select J1 = q; then [8. J1] [32. q] - q2[p2, q] + Ziq[p. q] - 21(1 - q2)p + 2q 2(J1 + J2), (3.8) {- ip(1 - q2) + 1(1 - q2)p} p(1 - q2)p [8. J2] i[p. q2]p(1 - q2)p ZJlB, (3-9) where J2 - 1(1 - q2)p. The eigenvalues of the matrix cjk are A(1) = l + (l + 4B)%, A(2) = 1 - (l + 4B)%. To avoid the square root M(M + 1) so that A(1) = 2(M + 1) and A(2) = - 2M. The we set B right-eigenvectors of cjk are 2(1) = (M + 1, 1), 2(2) = (M, - 1), hence J(1) J(2) q(M + 1) - 1(1 - q2)p. (3 10) QM + 1(1 - q2)p. (3 11) where we have made use of Equation (2.8). 28 Now f(B) = J(1)J(2) [q(M + l) - i(l - q2)p]J(2) qJ(23M - 1(1 - q2)pJ(2) szz + iqCI - q2)pM - 1(1 - q2)qu + (1 - q2)p(1 - q2)p M2 - (1 - q2)M(M + 1) + (1 - q2)p(1 - q2)p M2. (3.12) from which we conclude the first eigenvalue of M is m(1) = 0. Therefore, m(k) = k - l and by applying the relation b(k)= m(k)(m(k)+1) we then have b(k) = k(k - 1). (3.13) The first eigenvector u(l) satisfies J(2)U(1) = 0. (3.14) or in differential equation form (1 - x2)9§l-= 0, (3.15) where v1 corresponds to the unperturbed eigenfunction. Therefore, v1: constant is a solution of (3.15). Successive eigenvectors are de- termined from U(k+1) = J(1)ku(1), (3.16) 29 which in differential equation form becomes vk+1 = [(k + 1)x + (x2 - 1) j§.]kv1 k = Ck _9_k_ (x2 - 1)]" (3.17) d x where the C 's are constants. Normalization of the eigenfunctions Vk k yields _ 1 y - Vk(x) " (k ' 3).? pk_1(X), (3 18) 1 dk k where Pk(x) = -———--——-(x2 - l) are the well-known Legendre poly- k . k 2 k. dx nomials.17 The normalized eigenvectors Of (3.7) are then given by 1 8(k) = (k - l.)é Pk_1(q). (3.19) 2 The function g(B) is easily determined from the relation 3T1) = J(2)g(B). (3.20) We have 3(1) = (M + 1)q + ip(1 - qz) (M + 1)q + 1(1 - q2)p - 2q (M - 1)q + 1(1 - q2)p. (3.21) From (3.10) and (3.11) it follows that q = (3(1) + J(2))/(2M + 1) (3.22) and - 1(1 - q2)p = [J(l)M - J(2)(M + 1)]/(2M + 1). (3.23) 30 Inserting (3.22) and (3.23) into Equation (3.21) gives 3(1) = [J(1)M + J(2)(M - 2) - J(1)M + J(2)(M + 1)]/(2M + 1) = J(2)(2M - 1)/(2M + 1), from which we see that g(B) = (2M — l)/(2M + 1). Thus h(1)(M - 1) [f(B)g(B)]* M(2M - 1)x/(2M + 1);5 and h(2)(M - 1) = [fan/g(B)];é M(2M + l)%/(2M - 1)%. Summarizing results for the Operator B we have: (‘0- B = P(1- q2)p. n =1 (b). B = M(M + 1) (°)- J, = q = (J(1) + J(2))/(2M + 1) J2 =-i(1 - qzlp [J(1)M - J(Z)(M + 1)]/(2M + 1) (d). f(8) = M2. g(B) = (2M - 1)/(2M + 1) (M -1)(J(1) + J(2))/(2M + 1) - [J(1)M — J(2)(M + 1)]/(2M + 1) (3.24) (3.25) (3.26) (3.27) 31 (e). m(k) = k - 1, 6(k3 = k(k - 1) (k) = ' - l (f). 4 (k 2);: Pk_1(q) (g). h(1)(M) = (M + 1)(2M + 1);‘/(2M + 3);‘ h(2)(M) (M + 1)(2M + 3);5/(2M + 1);” (h). J(2)J(1) = h(2)(M)h(1)(M) = (M + 1)2 C. Matrix Representation for A The matrix representation for A = - aq2 can readily be obtained from the expression q2 = M2[(2M + 1)(2M - 1)]'1 + (M + 1)2[(2M + 1)(2M + 3)]'1 + J(1)2[(2M + 1)(2M + 3)]'1 + [(2M + 3)(2M + 5)]‘1J(2)2, (3.28) derived by employing (c), (d), and (h) from the above list of properties. From Equation (2.47) A . E(0) . E (J(1)rE(-r) . B(r)J(2)r) (3.29) r=1 so that we may also identify 3(0) (M) = - a(M2[(2M + 1)(2M - 1)]-1+ (M + 1)2[(2M + 1)(2M + 3)]'1} (3.30) E(’2)(M) = - a[(2M + 1)(2M + 3)]‘1, (3.31) B(2)(M) = - O[(2M + 3) (2M + 5)]‘1, (3.32) 32 while all other Ecr)(M) are zero. Hence, from Equations (2.50) -(2.52) and h(1)(M) = (M + l)[(2M + 1)/(2M + 3)]%, (3.33) h(2)(M) = (M + l)[(2M + 3)/(2M + 1)]%, (3.34) we have upon simplification A(O)(M) = - a(2M2 + 2M - 1)[(2M - 1)(2M + 3)]‘1, (3.35) A(-2)(M) = 3(2)(M) = - a(M + 1)(M + 2)[(2M + 1)(2M + 3)2(2M + 5)]‘5. (3.36) Since m(k) = k - 1, we obtain for the matrix elements A. = A(0)(m(j))6. + A(‘2)(m(k))5. + A(2)(m(j))6 (3.37) Jk 3k 3 k+2 j+2 k or finally D II - 2. 2 _ _ l_ + l _ 3 . kk 2 (k k 2)/[(k E-)(k i9]. (3 38) Ak k+2 ‘ Ak+2 k - 2.1(k + 1)/[(k + l)2(k - 1)(k + §)]%, (3.39) 4 2 2 2 and all other Ajk = 0. 33 D. Matrix Representation for Y The results of the two preceding sections are among those tabu- lated for linear operators by Green and Triffet.9 However, by making use of the relation q = (J(1) + J(2))/(2M + 1) we can now extend these results to include certain nonlinear operators by formulating the 2) normalized eigenvectors of the Operator B as functions of Jcl), J( , and M, i.e., <1“) = (1/2)? (3.40) 4‘” = (3/2)3°-q = (3/2)%(J(1) + J(2))/(2M + 1), (3.41) 4“) = % (5/2);§(3c12 - 1) .1 (5/2)%{3[(J(1) + 3(2))/(2M + 1)]2 - 1}, (3.42) 2 Keeping in mind that the matrix elements for Y must be calcu- lated by columns as determined by V(q), we have for the first column (Y = 3(1)) B(O)(M) = (1/2)%. (3.43) Therefore 5“” (0) (1/2)$5 (3.44) and the only nonzero matrix element in this column is 34 Y11 = (1/2)%. (3.45) For the second column (Y = 3(2)) E(0)(M) = 0, E('1)(M) = (3/2)%/(2M + 1), (3.46) so that by Equations (2.57) (2.62) the only nonzero matrix elements must be = 3(4) (0)h(1) (0) 12 = (1/2)%, (3.47) Y22 = 3(°)(1) = 0, (3.48) v32 = 3(‘1)(1)h(1)(1) g-(5/2)%. (3 49) For the third column = §.(5/2)%{M2[(2M + 1)(2M - 1)]"1 + (M + 1)2[(2M + 1)(2M + 3)]‘1- %. + J(1)2[(2M + 1)(2M + 3)]-1 + [(2M + 3)(2M + 5)]'1J(2)2}. (3.50) Thus E(0)(M) = %.(5/2)%{M2[(2M + 1)(2M - 1)]”1 35 + (M + l)2[(2M + 1)(2M + 3)]‘1 - é), 3(‘2)(M) = ;.(5/2)%[(2M + 1)(2M + 3)]‘1, (3.51) and the matrix elements become Y13 = E('2)(0)h(1)(0)h(1)(1) = (1/2)%, (3.52) v33 = 3(0)(2) = .3. (5/2)? (3.53) Y53 = 3('2)(2)h(1)(2)h(1)(3) = g (1/2);§. (3.54) Again, the elements not listed in a given column are to be taken as zero. Following the same procedure for ¢(4)3 ¢(S)3 ..., will de- termine as many columns of Y as may be needed. E. Summary of Results The first few matrix elements calculated for B, A, and Y as they appear in the nonlinear Legendre-like equation (3.1) are dis- played below: .1 F3 0 0 0 0 1 o ... 3 = 0 6 o ... . (3-55) 0 0 12 __ ............... J) 36 . 1/3 0 2/3/5' ...q A = - o 0 3/5 0 , (3.56) 2/3/5 0 11/21 ’1/72 1/72 1/72 ...1 0 0 0 Y = 0 2//1_0 2/577/‘2' . (3.57) 0 0 0 0 0 6/75 Combining these results gives the desired matrix representation for A, Ajk = b(J)5jk + ZnAannk. (3.58) A FORTRAN program, described in the Appendix, was used to find the eigenvalues and eigenvectors of the resulting nonsymmetric matrix. The eigenvalues corresponding to a = 1 were computed by truncating A to 10 rows and columns. These are listed in Table 1, along with comparable values computed by the standard perturbation method as presented, for example, by Saaty and Bram. The reason that every other value agrees exactly is that the nondiagonal elements of A are zero for every other row, hence the diagonal element is itself an eigenvalue. This is a peculiarity of the particular matrix being considered. The difference in the other 37 Table 1. A comparison Of eigenvalues for Equation (3.1) obtained by standard and matrix perturbation methods with a = 1. Standard Method Matrix Method k 10‘) k 10‘) 1 - 0.23570 1 - 0.24990 2 2.0000 2 2.0000 3 5.3977 3 5.3943 4 12.000 4 12.000 5 19.631 5 19.638 6 30.000 6 30.000 7 41.703 7 41.709 8 56.000 8 56.000 9 71.744 9 71.750 10 90.000 10 90.000 eigenvalues can be accounted for by the fact that standard pertur- bation theory presumes that the matrix A is essentially a diagonal matrix from the beginning, i.e., all nondiagonal elements are either zero or negligible in comparison with the diagonal elements. For such a matrix the eigenvalues can be approximated by the diagonal elements themselves. To assure that the nondiagonal elements are small conventional theory requires that the perturbation parameter a be small, thus forcing the desired condition. When a = l the nondiagonal elements of the matrix (3.58) are relatively small so that the two methods of computing eigenvalues give nearly the same results. But, as 6 increases in magnitude a great deal more error is introduced by neglecting these nondiagonal elements of A. To present a qualitative comparison of the two methods we shall examine how accurately the eigenvalues and eigenfunctions satisfy the original equation. The first eigenvalue and eigenfunction (k = 1) 38 computed by the standard theory are A1 = - 0.23570 (3.59) and u1(x) = v1(x) + 0.0351 v3(x), (3.60) where vk(x) is defined by Equation (3.18). Substituting these ex- pressions into Equation (3.1) yields the inequality |Ji.(1 - x2) 921 + Alul + x2u§| :_0.069, [x] :_1, (3.61) dx dx whereas an exact solution would reduce the right-hand side to zero. The corresponding eigenvalue and eigenfunction for the matrix method are A1 = - 0.24990 (3.62) and u1(x) 2 vl(x) + 0.0374 v3(x) + 0.0008 v5(x) + 0.0001 v7(x). (3.63) Substituting these values into Equation (3.1) yields |j%.(1 — x2) g?) + Alul + x2u§| :_0.004, |x| :_1. (3.64) Equation (3.60) has been limited to the sum of two unperturbed eigen- functions because that is all the standard theory predicts. On the other hand, the matrix method generates an infinite series of terms, the first four of which have been used for this calculation. To provide an example where the perturbation term dominates the remaining linear termswe have chosen the case where a = 5. Table 2 contains the eigenvalues obtained by truncating A to 25 rows and columns. 39 computed into Equation (3.1) gives d dx X I——-(l - x2) £21 + Alul + 5x2u§l :_0.44, le 5.1. Substituting the first eigenvalue and eigenfunction (3.65) This represents substantially more error than for the case when a = l, but it still may be considered a fair approximation for a first-order perturbation. NO attempt at comparison with standard perturbation theory is made, since the latter is not expected to give even a qualitative solution for such a large perturbation. In general the amount of error will increase with an increase in the magnitude of a, Additionally, for the nonlinear theory a significant amount of this error can be introduced when we approxi- mate equations like (3.2) with equations like (3.5). No similar approximation occurs in the linear perturbation theory, and it is this point that most clearly sets the one apart from the other. Table 2. given in Equation (3.58) with o = 5. Eigenvalues for the matrix k 1(k) k 1(k) 1 - 1.6896 14 182.00 2 2.0000 15 209.04 3 3.0686 16 240.00 4 12.000 17 271.08 5 18.309 18 306.00 6 30.000 19 341.07 7 40.591 20 380.00 8 56.000 21 418.97 9 70.758 22 462.00 10 90.000 23 504.78 11 108.88 24 552.00 12 132.00 25 598.51 13 154.97 4. GENERAL NONLINEAR SPECIAL-FUNCTION OPERATORS A more generalized form of the nonlinear Legendre-like equation, as well as other types of nonlinearities, may appropriately be treated by the perturbation method. It is not essential to re- strict considerations to equations Of the special-function type, but when other kinds of equations are featured expressions like (1.1) generally must be used for computing the matrix elements. This can sometimes result in a computational disadvantage. A. A Nonlinear Hypergeometric-like Equation Again consider the nonlinear hypergeometric-like equation introduced in Section 1, d 2 d d 2 —_. l - x .__ - 2 + vx __.+ ax u + A u = 0. 4.1 [dx ( ) dx (U ) dX k k] k ( ) For first-order perturbations the Operational form becomes [p(1 - q2)p + 2101 + 14011 - aq2V(q)¢(k) - 1m] 10‘) = o (4.2) with B = 3(1 - qzlp + 21(u + vq)p. r = 1, A = - 9612. y = V(q)¢(k). (4.3) Since most Of the special-function equations can be Obtained from the hypergeometric equation we feature it as a general example, the results of which may be applied to many other equations through appropriate transformations. 40 41 Following procedures similar to those of the last section we can obtain for the unperturbed Operator B the relations B = 3(1 - q2)p + 21(u + vq)p. n = (1 - q)'(“*v)(1 + q)'(V'“). 3 = (M - v)(M + v + 1), J = q = (J(1) + J(2))/(2M + 1) - uv/[M(M + 1)]. J2 = - i(1 - q2)p =[J(1)(M - 0) - J(2)(M + \) + 1)]/(2M + 1) + “(M - V)(M + 6 + l)/[M(M + 1)], f(B) = M2 + u262/M2 - 02 - 62, 8(3) = (2M - 1)/(2M + 1), 6(k) = (m(1) + k + v)(m(1) + k - v - 1), where m(1)2 + 3262/m(132 = “2 + V2. Details Of computing the matrix elements of A and Y will be omitted since they are analogous to those given in Section 3. To represent A = - aqz, however, we note that q2 = R1(M) + J(1)2R2(M) + R2(M + l)J(2)2 1. J(1)33(M) + R,(M)J(2) (4.4) where the following have been defined for notational convenience: 31(M) R2(M) Rch) R.(M) Thus we By setting M = m element 42 f(B)/[(2M + 1)(2M - 1)] + h(2)(M)h(1)(M)/[(2M + 1)(2M + 3)] + 32v 2/[M2(M + 1)2], [(2M + 1)(2M + 3)]-1, - uv{[M(M + 1)]”1 + [(M + 1)(M + 2)]'1}/(2M + 1), - uv{[M(M + 1)]"1 + [(M + 1)(M + 2)]‘1}/(2M + 3). identify 3(0) = — GR1(M), A(1) . A(-1) = - 0R3(M)h(1)(M), A(2) . A(~2) (k) m(k) = m(1) , k _ 1, 9 5 _ (k) Akk -"aR1(m ), A = A k k+l k+1 k 01R:‘3(m(k))[m(k+'1)2 + pzvz/mowl)2 - “2 - v2]$2 x [(2ka) + l)/(2m(k) + 3)]5 -OR2(M)h(1)(M)h(1)(M + 1). (4. (4 (4. (4 (4. (4 (4. we calculate the matrix (4. (4. 5) .6) 7) .8) 9) .10) 11) 12) 13) 43 A A k k+2 k+2 k _ OR2(m(k))[m(k+1)2 + u2v2/m(k+1)2 _ u2 _ v2]% x [m(k+2)2 + quZ/m(k+2)2 _ u2 _ vzfg x [(2m(k) + l)/(2m(k) + 5)]%. (4.14) All elements not listed are to be taken as zero. For the Operator Y we start by assuming that the normalized .i‘ “t—rfi W- IND-tn.” i eigenvectors of B are of the polynomial form N ((1‘) = Z Cmq“. (4.15) n=0 n where the Cik) are constants. When this is not the case the eigen- vectors can always be approximated by a polynomial through the use of a truncated Taylor series. For the first column (Y = ¢(1))we have E(O)(M) = 031) (4.16) so that the only nonzero element in the first column is _ (1) Y11 _ c0 . (4.17) The second column has Y .__ ,(2) C52) + C£2)q cgz) + cf2){- uv/[M(M + 1)] + J(1)(2M + 1)‘1 + (2M + 3)'1J(2)}, _ (4.18) 44 from which we identify E(0)(M) - C52) - Cf2){uv/[M(M + 1)]}, E(’1)(M) = Cf2)(2M + 1)'1. (4.19) Thus the matrix elements are Y12 ' = Cf2)[m(2)2 + quZ/m(2)2 _ u2 _ 021% x [(2m‘1) + 1)(2m(1) + 3)]“2. (4.20) YZZ = E(0)(m(2)) = ng) - Cf2)[uv/m(2)m(3)], (4.21) Y32 = 5"1)(m(2))h(1)(m(2)) (31(2) [m(3)2 + 1J2V2/m(3)2 _ U2 _ v2];§ x [(2m(2) + 1)(2m(2) + 3)]'%. (4.22) Subsequent columns can be computed in the same manner. B. Other Forms of Nonlinearity Theoretically, matrix representations can be found for any nonlinear perturbation for which a Taylor series expansion like (2.54) exists. Listed below are more examples of nonlinear special-function operators illustrating certain other forms of nonlinearities. In addition to the first few matrix elements of the perturbations, we have tabulated the operator relationships essential for obtaining the matrix elements. Details of the calculations are omitted since they are analogous to those already presented. 45 1. Nonlinear Legendre-like Harmonics d 2 du 2 du -—- l - x + 20 l - x u + A u = 0, dx ( ) dkk ( ) Eik k k k _. (k p(1 - 42)p - 14(1 - q2)pch)¢ ), >* H a! H p(1 - q2)p = M(M + 1), F = 1, - 10(1 ' q2)P3 D) H V(q)¢(k), ..< H + q = (J(1) + 1(2))/(2M 1), - i(1 - q2)p = [J(1)M - J(Z3(M + 1)]/(2M + 1), 3(3) = M2. g(B) = (2M - 1)/(2M + 1). b(k) = k(k - 1), m(k) = k - 1, (k) _ _ )_3 ¢ ‘ (k 2) Pk_1(Q)3 where Pj(q) is the j-th Legendre polynomial, ak(k - 1)/[(2k - 1)(2k + 1)]%, Ak+1 k __ _ 16 AR k+l — ak(k + l)/[(2k 1)(2k + 1)] . The matrix elements of Y are exactly those given in Section 3 for the nonlinear Legendre-like equation discussed there. 46 2. Nonlinear Hermite-like Harmonics 2 d u 2x 22k 2k dx dx + OXUE + Akuk = 0, (k)2 2 aqVIq)¢ . A p + Ziqp 2 C17 ll 2M, 3 + 2199 r = 1, A='aq: Y = VIQ)¢(k)23 ip = ch), f(B) 2M, g(B) = 1, (k) (k) b 2(k - 1), n) = k — 1, k- -9 4(k) [2 11550 - 1):) zuk_1(q). where Hj(q) is the j-th Hermite polynomial, Ak+1 k = Ak 3+1 L = " “(k/2) 2, _ ;2 - Y11 - (l/n) , le — 0, ... '8 Y21 = 0, Y22 = (9/n) , ... Y31= 0, Y32 = 0, 000 47 The above unperturbed Operator B = p2 + Ziqp is related by a similarity transformation to the familiar harmonic oscillator Operator B' - p2 + q2, which we have featured as an illustrative example in Section 2. 3. Nonlinear Laguerre-like Harmonics x $211 + (1 -x) 533;). + 01qu sin (uk) + Akuk = 0, A = 932 - i(1 - q)p - <1qu sin (41m). " :. B=qu-i(1-q)p=M, r = 1, A = - 99. Y = V(q) sin (40”), q = (3(1) + 3(2)) + 2M + 1, {(3) = M, g(B) = 1. 13m =k-1,m(k)=k-l, ¢(k) = Lk_1(q), where L. (q) is the j-th Laguerre polynomial, J Akk = - 01(2k - 1), Ak+1 k = Ak k+1 = - “1‘35, 0.0, 48 -1- sin (1), 2 0, l - -2- cos (l), .00, 5. OTHER EXAMPLES A. The Hartree Equations for the Helium Atom Thelkutree equations for the helium atom are classic examples of nonlinear equations arising because of the interactions between particles. They follow from applying the theory of the self-con— . .18 . . Sistent field to the Ham1lton1an H=— V2+V2 -.i-i+__za Sol (1 91.11.21,”. () where};1 and}:2 are the coordinates of the two electrons relative to the nucleus, r12 = r21 represents the distance between the electrons, and a is a parameter related to the strength of the Coulomb repul- sion between the electrons. Spin-dependent forces are neglected. The central-field postulate is the general principle under- lying the Hartree selfeconsistent field. It is assumed that each electron moves independently of the other in a central potential representing the Coulomb potential of the nucleus and the averaged Thus the time-independent potential field of the remaining electron. wave function of the helium atom is written as 1165,. 5,) = u(1)(;,)u(23(3,), (5.2) where (j) (J') (J') - _ u (£J.)~Rm (rjnrm (ej. 4j). 3 - 1. 2 (5.3) designates the one-electron orbitals. Here r, e, (p are the standard 49 50 spherical coordinates, the Y3!) are the spherical harmonic functions 18 and the R0) are radial functions. n2 Applying the principles of the calculus of variations to the total energy of the atom we obtain the two nonlinear radial equations u: 7' .‘b. 1L (. . . [— d2 +££LLll - i + 39 S(r , n, 1)] 3n? = 33))13), (5.4) O??- 1‘? I'. I‘. i J 3 J J E for j - l, 2, (L.- (J') (3') Pn1 (rJ) rjan (rj), (5.5) and 1'. - 0° . r. 8(rj. n. 1) = ((3123 (rinzeu-i + 1[(13%) (111)]2 (i) dri. (3.6) J where i = l, 2, but i 75 j. When a = 0 Equations (5.4) both become hydrogenic radial equations with the Operational form [p2 + 2(2 + 1)q'2 - 4q'1 - b(k)]U(k) = 0. (5.7) The operator B = p2 + 2(2. + l)q’2 - 4q"1 does not lend itself to the calculation of the step operators J0) and J(2) as we have defined them. This is a consequence of the fact that the discrete eigenvalues of B are bounded above, as will soon be made clear. We consider the transtrmation k _ ()3. _ _ (10) a p (M) , 51 Q? - P0 = i, which transforms the operator B to 3' = ()32 + Q + 2(2 + 1)Q'1, (5-9) with the associated eigenequation 3'u(k) = 4(- (300)-1/2 u(k). (5.10) Hence, the eigenvectors of B' in the coordinate Q are the same as for B in the coordinate q. Following the procedure of Section 2 the step operators of B' are determined to be J(1) Q - iQP - M, (5.11) J(2) Q + iQP - M. (5.12) Utilizing these we may obtain the relations 3' = 2M, J1: Q = 10(1) + J(2)) + M, 2 J2 = - iQP = %.(J(1) - J(2)), f(B) = M(M - 1) - 2(9. + l), g(B) = l, b'(k) = 4(- b(k))’% = 2(k + 1), (k) 1+1 -Q 2£+1 U = 2 , 5.13 (29) e 1M (Q) ( ) m . 17 where n = k + 2 and the L.(x) are generalized Laguerre polynomials. Thus, the eigenvalues and normalized eigenvectors of B are given by i.‘ 52 b(k) = - 4/“2, n = 1, 2, 000, (5.14) 2+1 (k)__ 2(n - R. - 1)! 1(1) -2q/n 2£+1 " L 4 . 5.15 ¢ “2““ 4' 2):]3 n e n+2. ( q/n) ( ) The perturbing operators A and Y for Equations (5.4) also have similar forms defined by A = zoq'l, (5.16) and Y = V(q)S(q. 11. 1). (5.17) Because the discrete eigenvalues b(k) have the upper limit (k), of the operator B, do not form zero the discrete eigenvectors 6 Therefore, as suggested a complete system of orthonormal vectors. above, they do not provide a suitable basis for a matrix representa- tion. The eigenvectors associated with the operator B' do form a complete orthonormal system, however, and it is this set that we shall utilize for our matrix representations, even though the Operator B will no longer be diagonal. To calculate the ground state energy of the helium atom we seek to diagonalize the matrix A = B + AY, (5.18) where B represents the energy of one electron due only to the Coulomb potential of the nucleus and AY expresses the energy of interaction of the electrons. To the first eigenvalue Of A must be added the energy 53 (1) b of the other electron due to the potential of the nucleus alone. Our matrix representation will be based upon the condition that 2 = 0 even though this is not entirely correct whenever n > 1; there is a degeneracy in the 2 quantum number but the major contribution to the energy occurs for 2 = 0. Also for 2 = 0, the operator Y can be approxi- mated by18 Y = 1 — (1 + 2q/n)e-uq/n. (5.19) Because the form of the perturbation does not lend itself to the operator method of computing matrix elements in this case we have utilized an analytical technique, a general discussion of which is presented in Section 6. The matrix representations for B and AY can then be defined by G) , k 3. = f ¢(J)B¢( ) dq (5.20) Jk o and ” (') -1 (k) (AY)jk = 26 I ¢ J q S(q, k, 0)¢ dq, (5.21) o k) . . where 9 represents the norma11zed eigenvectors of the operator B'. To be explicit, the first matrix elements of (5.20) and (5.21) are given by (a = l), w p—a H I - - 4 f e'2q(q2 + 2q) dq O = - 3, (5.22) 54 (AY)11 = 8 J [e-qu - e-6Q(2q2 + q)] dq 0 at 1.6296. (5.23) Continued calculations will finally yield the representations -3.0000 1.1550 3 = , (5.24) 1.1550 -l.6667 and 1 1.6296 -0.3512 AY = . (5.25) -0.6842 1.0000 The first eigenvalue of A = B + AY is E(1) 10 = - 1.724 Rydbergs, (5.26) so that the total ground state energy becomes (1) (1) E10 = E10 + b = - 5.724 Rydbergs. (5.27) Compared with the experimental value of - 5.807 Rydbergs, this rep- resents an error of 0.083 Rydbergs. If larger matrix representations are desired, then the approximation (5.19) should probably not be used, as its validity beyond n = 2 is questionable. The calculation of Y from expressions like (5.6), however, is extremely tedious. Millman and Keller19 have computed the ground state energy for the helium atom by applying a modified form of the standard pertur- bation method to the Hartree equations. They obtained a value of 55 - 5.500 Rydbergs, representing an error of 0.307 Rydbergs with the experimental value. In addition, they made a comparison of results with those achieved by applying conventional perturbation techniques to the SchrOdinger equation and showed that the eigenvalues agree up to first-order terms, but that the eigenvectors agree only up to zero-order terms. ’ E? 15.3.2,- By the present method, however, it can easily be shown that the matrix representation for the (radial) SchrOdinger equation, whose F b - Hamiltonian is given by (5.1), becomes Hr = 23 + AY, (5.28) where B, A, and Y are the same as defined for the Hartree equations and 1 = 0. Clearly, the eigenvalues of (5.28) will not yield the same values as for the Hartree equations since the sum of the eigen- values of two matrices is not the same, in general, as the eigen- values Of the sum of the same two matrices. B. The van der Pol Equation An equation of frequent occurrence in the field of nonlinear oscillations is the van der Pol equation. A general form of this equation is given by H u - a(l - u2)u' + Au = 0, A > 0, (5.29) where the primes indicate derivatives with respect to time. Up to this point we have concentrated entirely on the problem of finding those values of the parameter A for which nontrivial solutions of the 56 differential equation exist, subject to certain boundary conditions. In nonlinear oscillation theory equations such as (5.29) are treated as initial-value problems, where one seeks periodic solutions for a fixed value of A. Because of this distinction, the matrix perturba- tion method featured here must be applied in a slightly different way. To be precise, we look for solutions of (5.29) with A = l satisfying u(t + 2n/w) = u(t) (5.30) and subject to the initial conditions u(O) = 2, u'(0)= 0. (5.31) The existence of such solutions is guaranteed by the theorem of Liénard;20 and stability requirements indicate that the amplitude of u(t) varies between +2 and -2. The unperturbed solution of (5.29) subject to the initial con- ditions (5.31) is given by uo(t) = 2 cos wot, (5.32) where ”o = /A'= l is the angular frequency. The effect of the non- linear term in the van der Pol equation is to change the angular frequency “0 to a new value m, but we can calculate w by expanding A in a Taylor series about a = 0, A=wg + 0101+ ... , (5.33) then setting A = l to obtain 01=V1 - cpl, (5.34) ”I. m 57 where 01p1 represents the first-order perturbation in the eigenvalues of (5.29). In order to obtain a complete set of eigenfunctions, from which appropriate matrix representations can be derived, we reformulate (5.29) as a boundary-value problem. Choosing 0 :_t :_n as the funda- mental domain Of the Hilbert space and imposing the boundary condi- tions u(O) = 2, u(n/w) = - 2, (5.35) which are consistent with stability requirements and the initial conditions (5.31), then the eigenvalues and normalized eigenfunctions of the unperturbed equation are 13k = k2 (5.36) and vk(t) = C cos kt, C = VZ/n , (5.37) respectively. The nonlinear Operator Of‘(5.29) is identified as 2) .9. (5.38) N0 = (l - v kdt to a first-order perturbation. To compute the matrix elements of N0 it is again most convenient to apply the method discussed in Section 6. Thus, we write 71’ Ngk = - kC2 f [cos jt(l - C2 coszkt) sin kt] dt, (5.39) 0 which upon evaluation yields 0 o _ Nkk = N3k k - 0. (5.40) I‘- n 58 and for j # k or 3k, 0 k2C2 i+k 4 - 02 302 JR 4 32 .. k2 9k2 _ j2 where C2 = 2/n. Diagonalization of the matrix A_ = b,6_ + ON? will produce Jk J JR Jk an infinite set of eigenvalues A(k) and associated eigenvectors, from which we can extract solutions to the van der Pol equation for all cases where the fixed parameter A = k2. The first eigenvalues are ,(1) 1.0539 (5.42) and 1(1) = 1.2161 (5.43) for a = 0.5 and o = 1, respectively. From these, the first-order perturbations are determined to be op1 0.0539 (5.44) and 0.2161; (5.45) 061 hence, by applying (5.34) the perturbed angular frequencies become 0.9726 (5.46) E II and 0.8854 (5.47) E II for a = 0.5 and a = 1. Ta. . ‘mn— 59 Since in this case it has been specified that A = l, the desired solution can be expressed in the form u(t) = D0 + D1u(1)(t), (5.48) l where u( )(t) is the perturbed eigenfunction belonging to the first (1) eigenvalue A , while D0 and D1 are two constants to be evaluated from the boundary conditions (5.35). Solutions corresponding to a = 0.5 and a = 1.0, respectively, are given by u(t) 3 0.1599 + 1.9921 cos wt - 0.1585 cos 2wt + 0.0095 cos 3wt + 0.0006 cos 4wt + 0.0013 cos Smt - 0.0004 cos 6wt (5.49) and u(t) 2 0.3121 + 1.9582 cos wt - 0.3095 cos Zwt + 0.0371 cos 3wt - 0.0005 cos 4wt + 0.0053 cos Smt ~ 0.0014 cos 6mt - 0.0005 cos 7wt - 0.0005 cos 8mt, (5.50) again retaining only the first few terms of the infinite series. These solutions give values of u(t) and u'(t) which are comparable with those obtained through the use of time-consuming numerical techniques.20 This comparison is presented for half a period in Table 3; by periodic extension the entire solution can be constructed. Attempts were made to obtain solutions for larger values of the parameter a, but it was discovered that the error in the perturbed angular frequency w grew very rapidly for o > 1. However, if one extended to second-order perturbations, qualitative results for larger 60 Table 3. Approximations of u(t) and u'(t) for the van der Pol equation obtained by numerical and matrix techniques. w ‘“ . ' “r" ‘‘‘‘ =.===-a a = 0.5 a = 1.0 Numerical Matrix Numerical Matrix Method Method Method Method t u(t) u'(t) u(t) u'(t) t u(t) u'(t) u(t) u'(t) 0.0 2.00 0.00 2.00 0.00 0.0 2.00 0.00 2.00 0.00 0.1 1.99 -0.19 1.99 -0.14 0.1 1.99 -0.17 1.99; -0.09 0.2 1.96 -0.34 1.97 -0.28 0.2 1.97 -0.30 1.98 -0.18 0.3 1.92 -0.48 1.93 -0.42 0.3 1 1.93 -0.40 1.96. -0.26 0.4 1.87 -0.60 1.89 —0.55 0.4 1.89 -0.47 1.93 -0.35 0.5 1.80 -0.71 1.82 -0.68 0.5 1.84 -0.53 1.89 —0.43 0.6 1.73 -0.80 1.75 -0.81 0.6 1.78 -0.59 1.85 -0.51 0.7 1.64 -0.89 1.66 -0.94 0.7 1.72 -0.64 1.79 -0.59 0.8 1.55 —0.98 1.56 -1.07 0.8 1.65 -0.68 1.73 -0.67 0.9 1.45 -1.07 1.45 -1.19 0.9 1.58 -0.73 1 1.66; -0.75 i 1.0 1.33 -1.15 1.33 -1.31 1.0 1.51 -0.78 1.58 -0.84 1.1 1.21 -1.24 1.19 -1.43 1.1 1.43 -0.83 1.49 -0.92 1.2 1.09 —1.34 1.04 -1.54 1.2 1.34 -0.89 1.39 -1.02 1.3 0.95 —1.44 0.88 -l.65 1.3 1.25 -0.96 1.28 -1.12 1.4 0.80 -1.54 0.71 -1.75 1.4 1.15 -1.04 1.17 -1.24 1.5 0.64 -1.65 0.53 -l.83 1.5 1.04 -1.12 1.04 -1.35 1.6 0.47 —1.77 0.35 -1.91 1.6 0.92 —1.23 0.89! -1.47 1.7 0.28 -1.88 0.15 -l.96 1.7 0.80 -1.35 0.74 -l.58 1.8 0.09 -2.00 -0.04 -2.00 1.8 0.65 -1.49 0.58 -1.69 1.9 -0.11 -2.10 -0.24 -2.01 1.9 0.50 -l.65 0.40 -1.78 2.0 -0.33 -2.18 -0.44 -2.00 ‘2.0 0.32 -l.83 0.22 -l.86 2.1 -0.55 -2.22 -0.64 -1.97 2.1 0.13 -2.04 0.04 -1.91 2.2 -0.77 -2.22 -0.84 -1.91 2.2 -0.08 -2.25 -0.15 —1.94 2.3 -0.99 -2.15 -1.03 -l.83 32.3 -0.32 -2.46 -0.35 -1.95 2.4 -1.20 —2.02 -1.20 -1.72 2.4 -0.58 42.62 -0.54 -1.94 2.5 -1.39 -l.83 -1.37 -1.59 2.5 -0.84 -2.68 -0.74 -1.91 2.6 -l.56 -l.58 -1.52 -1.43 2.6 -1.11 -2.59 -0.93 -1.85 2.7 -1.71 -1.29 -1.65 -1.25 2.7 -1.35 -2.34 -1.11 -1.77 2.8 -l.82 -1.00 -1.77 -1.05 2.8 -1.57 -1.95 -l.28 ~l.67 2.9 -1.91 -o.71 -l.86 -0.82 2.9 -1.74 «1.48 -1.45 -1.54 -l.96 -0.43 -1.93 -0.59 -1.87 -1.02 -1.59 -l.38 -1.99 -0.19 -1.98 -0.33 -1.95 -0.62 —1.72 -1.19 -2.00 0.02 -2.00 -0.08 -1.99 -0.29 -l.83 -0.97 -1.99 0.20 -1.99 0.18 -2.00 -0.04 -1.91 -0.71 -1.96 0.36 -1.96 0.44 -2.00 0.14 -1.97 -0.44 -1.92 0.50 -1.91 0.68 -1.98 0.28 -2.00 -0.14 lllllll'il'lll ..... 61 values of a could be obtained. Higher-order perturbations are dis- cussed in Section 6, but extensive calculations would be required to yield the additional matrix elements. As a final comment it could also be mentioned that the results given here for the van der Pol equation were calculated by treating only first-order perturbations, whereas for the standard method the second-order perturbation must be included before a change in angular frequency is detected. “airliner, 6. GENERAL REMARKS A. Higher-Order Perturbations So far we have concerned ourselves only with first-order per- turbation effects. For most of the applications we have in mind, where the nonlinear term is of the same order of magnitude as the linear terms (a = l), first-order corrections obtained by this particular method are probably sufficient. However, when the perturbing term be- comes large in comparison to the other terms (a > 1), a second-order correction to the unperturbed solution may be necessary in the interest of accuracy. This second-order correction will appear in the form of another operator added to A, A=B+A(Y+Q), (6.1) where B still designates the unperturbed linear operator, AY is the first-order nonlinear perturbation of B, and A0 represents the second- Order perturbation. Here F is assumed to be unity, since it is not essential to the following discussion. To obtain the explicit form of the operator 9 from which its matrix elements can be calculated, we reconsider the general nonlinear equation (B + AZ - A(k))w(k) = 0. (6.2) As before in Section 2(D) we expand the operator Z in a Taylor series about a = 0, 62 IB-Oh' a.“ ~M.“ Fit. I (‘41-, it 63 N ll Y + 021 + 0222 + one , (603) where again _ 3’}. (k) 21 — 30m ()1 a , , (6.4) II C and insert this into Equation (6.2), (k) (k)] d) = 0. (6.5) [B + A(Y + OZ1 + 012Z2 + ...) - A Retaining not only Y but also the term 0Z1, Equation (6.5) becomes [3 + A(Y + 0121) - 1m] ((1‘) = 0, (6.6) so that 0 = OZ Assuming the first-order perturbation problem to be (M solved, we will know 61 , hence 21, and can therefore write 1. (k) (A) 32 0 = Z .3—(k) , (6.7) 2#k 2 W a=0 where yik) represents the vector components (normalized in such a way that yfik) = l) of the eigenvectors u(k). In Obtaining (6.7) we have used the expansion 1(k) = ¢(k) , awfk) = ¢(k) + Z Y(10,01) (6.8) 2#k 1 for the first-order terms. A matrix representation for Q can now be (Obtained in a manner similar to that defined for Y. Hence, A can be eXpressed as jk jk + nAjn(Ynk + an)‘ (6'9) 64 There is, of course, one form of Z which simplifies some of the computations required to obtain the matrix elements of O. This occurs when Z = V(q)w(k) so that _. 2 Q = 2 Y£(k)V(Q)¢( ); (6°10) I R.#k l in particular, 91k = O, A. (1) . Qj1=Yj //2—9]#19 (k) = . 6.11 Qkk 2§k Y9. Yik ( ) Some of the matrix elements of 0 for the nonlinear Legendre- like equation discussed in Section 3 are displayed below; for the case when O = 1 we obtain: 0 O 0 one " O 00. Q - 0 0.0592 (6.12) - 0.0264 0 -0.0303 ... n- d) To compare the results of first-order to second-order perturba- tions we will treat 3 x 3 matrix representations for the nonlinear Legendre-like operator. The first eigenvalue and eigenfunction are A1 = - 0.2495 and u1(x) = v1(x) + 0.0373 v3(x), Ill! iii 65 for the first-order perturbation, while for the second-order pertur- bation Obtained by diagonalizing (6.9). A1 = - 0.2579 and u1(x) = v1(x) + 0.0396 v3(x). Substituting these values into Equation (3.1) gives (1 2 du 22 ... l - x ——J + A u + x u ldx ( ) dx 1 1 1| IA 0.042, IxI_<_ 1 and .39; (1 - x2) 3.31 + Alul + x2u§| _<_ 0.018, |x| _<_1 for first-order and second-order perturbations, respectively. Following similar procedures one can in principle define third-, fourth-, and higher-order perturbations. However, obtaining the matrix elements in such cases may well take too long to be practical. B. Analytical Methods As we have seen there are certain differential equations for which the operator technique of computing matrix elements is not con- venient. In these cases the usual analytical method has been applied. Perhaps a few comments on the latter approach and its relation to the Operator technique would be helpful. Again consider the nonlinear equation (2.1), (L + ON + A)u = 0. (6.13) 66 Assuming that u = u(x) is also a continuous function of O, we form the Taylor expansion about O = 0, L1 = U0 + GT1+ 0.. (6.14) where uo represents the unperturbed solution of Equation (6.13) and OT1 is the first-order perturbation. Thus it follows that Mu,u,..J =N0.agflrl+ ”,, (613 1.1 where we define N0 = N(x, uo, ...). Substituting (6.15) into (6.13) and retaining only terms linear in O we have (L + ONO + A)u = 0. (6,16) Suppose that the unperturbed form (O = 0) of Equation (6.16) has eigenvalues and normalized eigenfunctions given by )0 = bk, (6.17) u0 = vk(x), (6.18) respectively. Expanding the perturbed eigenfunctions u(x) in terms of the unperturbed eigenfunctions vk(x) yields u = kckvk. (6-19) Substitution of (6.19) into (6.16) then gives 0 _ Zk(L + on + A)ckvk- 0. (6.20) Multiplying (6.20) by vj and forming an inner product over the (l I III I III I ll lullool i I. l l l 'Il.ll '. . 67 fundamental domain we get 0 . - N = A 6 .2 zk[(vj, ka) + O(vj, vk)]ck ck jk (6 l) or — L N0 = A 6.22 Zk( jk + a jk)ck cj. ( ) where ij - (vj, ka) (6.23) and o Njk — (vj, Novk). (6.24) However, L has the diagonal representation ij = - bjsjk’ (6.25) so that (6.22) becomes 0 = Zk(bj8jk - onk)ck ch. (6.26) Thus diagonalization of the matrix A, = b ONO will 3k )‘Sjk' jk yield the perturbed eigenvalues A. The matrix A referred to here is the same as that Obtained by the operator method. It is the form of the differential equation that will generally determine which method should be adopted for the computation of the matrix elements. A comparison with the standard perturbation theory can easily be made at this point. If, instead of multiplying Equation (6.20) by vj, we multiply by Vk and form the inner product then (6.21) becomes 68 Ek[(vk’ ka) + a(vk, Novk) + A]ck = 0 (6.27) or {kwk — onik - A]ck = 0. (6.28) Thus A = 6k — onik, (6.29) which clearly defines the perturbed eigenvalues A, for the standard perturbation method, as the diagonal elements Of the matrix A. BIBLIOGRAPHY 12. 13. 14. 15. BIBLIOGRAPHY J. D. Talman, Special Functions (W. A. Benjamin, Inc., New York, 1968) Math. Phys. Monograph Series T. L. Saaty and J. Bram, Nonlinear Mathematics (McGraw-Hill Book Co., Inc., New York, 1964) International Series in Pure and Applied Math. L. Brand, Differential and Difference Equations (John Wiley and Sons, Inc., New York, 1966) A. Messiah, Quantum Mechanics (John Wiley and Sons, Inc., New York, 1961) W. Heisenberg, Introduction to the Unified Field Theory Of Elementary Particles (Interscience Publishers, Inc., New York, 1966) L. de Broglie, Introduction to the Vigier Theory of Elementary Particles (Elsevier Publishing Co., NewFYOIk, 1963) Engli§h trans. by A. J. Knodel P. A. M. Dirac, The Principles of Quantum Mechanics (Oxford Univer- sity Press, London, 1958) H. Goldstein, Classical Mechanics (Addison-Wesley Publishing Co., Inc., Mass., 1950) H. S. Green and T. Triffet, J. of Math. Phys. 10, 1069 (1969) L. Infeld and T. E. Hull, Rev. Mod. Phys. 23, 21 (1951) W. Miller, Jr., Lie Theory and Special Functions (Academic Press, New York, 1968) A. Joseph, Rev. Mod. Phys. 39, 829 (1967) C. A. Coulson and A. Joseph, Rev. Mod. Phys. 39, 838 (1967) E. C. Kemble, The Fundamental Principles of Quantum Mechanics (Dover Publications, Inc., New York) R. Courant and D. Hilbert, Methods Of Mathematical Physics, Vol. I (Interscience Publishers, Inc., New York, 1953) 69 16. 17. 18. 19. 20. 70 H. Sagan, Boundary and Eigenvalue Problems in Mathematical Physics (John Wiley and'Sons, Inc., New York, 1966) P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill Book Co., New York, 1960) International Series in Pure and Applied Physics J. C. Slater, Quantum Theory of Atomic Structure (McGraw-Hill Book Co., New York, 1960) International Series in Pure and Applied Physics M. H. Millman and J. B. Keller, J. of Math. Phys. 10, 342 (1969) H. T. Davis, Introduction to Nonlinear Differential and Integral Equations (Dover PublicatiOns, Inc., New York, 1962) APPENDIX APPENDIX Computer Program Listed below is a modification of a FORTRAN prOgram originally developed by Green and Triffet for computing the eigenvalues and eigen- vectors of a real matrix.9 Either symmetric or nonsymmetric matrices may be entered; once the matrix has been introduced, a check is made to determine whether it is symmetric or not, then the appropriate branch of the program is automatically selected. The original program was written for a CDC 6400 computer, and some difficulties with accuracy arose when the program was run on an IBM 360 computer. To compensate for this, double precision calculations were introduced. Computation time was lengthened, but a 25>(25 matrix still only requires about 5 minutes. A critical point in the program is the estimate of the largest eigenvalue. This appears to be even more important when multiple eigenvalues exist. An estimator technique, based on an iteration method for computing the magnitude of the largest eigenvalue, is built into the program. Before printout, a check is made of the largest eigenvalue by means of the estimator. If these values are not suffi- ciently close, the largest eigenvalue becomes the estimator and the procedure is repeated until convergence is reached - a process which may be accelerated, if necessary, by reducing the tolerance in FUNCTION DLEIG. 71 72 Amom.mcm9_m3 Aa._vauhe._c<>520 oz_z_am ;o_m_owaa msmaoo OPO\ oauws_e.oo moex 73 z.Hua_ Nam co we._e 2 use m. 305mm >sme<_cm:z_..xfiace..xmfivwico_a..xe.m_.\.orva<:mou saw Aa.ec<.efisfim.ecu9_mz 2.Hue mHm on 2». Hfim Ham.HHm.oomn2_ca_ can Hom.ofim.oflmhau_va_ 1-2». mom mom.~om.momhz-_va_ H+_u_ mam _e.mm-A<_.ac<.emuh<_.ev< “om A<_.63<.mm+_e.cmuaq_._c< A<_._cau_< 75 Au.uau.ma3 c_dgo Lcau m0h<:_Hmm .VHmze x..a:ou Lee m.a.<>z.c.m .. .m..xaa.quaoa Nos .Nos.e.we.a3 .......e_.caz .eaoc<.o i.ez use 30m .....x~fi.\..meaaza.m ..< :e.3 x.aem 2 en... ma» .0 .aroc<_o H 3.42 O.» :o 3.3aaa >axe .ama .m< mm:.<>zae.w are .. .N..XHH.pe am a macmq >.oea.o.;s. .H..XHHcpmkma cz< «OH<:_Hmm th m< u344>zwo_w owhq4304zuc_u kmucz<4 >44qc_zwo_w xw.a~ ace..x.H.\..moa aquaa< mwaae_zc.zo .x.a»m z a. zoowm are 2. mzz:.oc >m o.em.. me< mmoecm>zmc_w wze .s..xHHOhmmaa .c.ma.mmmmae-c.ma.mmmmaa .mmwzu>zoo kgsmmm 4_H23 mhqadh_ oz< < >m m >40_k432 oo.Hu ..Om 2.Hu. m on oono.Hu.OH .z.<.c.w.c zo.eez:a zo.m.omua m.mzoo mm cm ma ma HH OH 000 82 ozm zmapmm mH.mH.H .x<2.- .... O221+ _ u _ Hmmm\.ea..< u Haa_.< 6+.uaa. z.Hua cH on e.mH.e .Hmmmc a. ...ea..<.me