The lanczos algorithm and spline expansions

Nuclear Physics A463 (1987) 139c - 144c North-ltolland, Amsterdam

139c

THE LANCZOS ALGORITHMAND SPLINE EXPANSIONS G.L. PAYNE Department of Physics and Astronomy, The U n i v e r s i t y of lowa, lowa C i t y , lowa 52242, U.S.A.*

1. INTRODUCTION There e x i s t many techniques to solve the Schrodinger equation for three nucleons with r e a l i s t i c two- and three-body i n t e r a c t i o n s .

The solution of the

Faddeev equations in configuration space has the advantage that the equations have a simple form, which can be programmed in a straightforward manner. difficulty

The

with t h i s method is that the standard techniques for solving the

p a r t i a l d i f f e r e n t i a l equations can lead to numerical c a l c u l a t i o n s which are p r o h i b i t i v e even on modern supercomputers.

In t h i s paper we review two tech-

niques which have been used to reduce the problem to one which can be solved e a s i l y on a CRAY or s i m i l a r computer.

T r a d i t i o n a l l y , Faddeev c a l c u l a t i o n s ,

both in configuration space and momentum space, are categorized by numbers of channels.

These channels specify the angular momentum quantum numbers of the

i n t e r a c t i n g pair of nucleons, and the corresponding quantum numbers of the remaining spectator nucleon with respect to the i n t e r a c t i n g pair.

To obtain a

converged solution to the bound-state problem, one must include enough channels to accurately represent the three Faddeev amplitudes. t i o n s I in configuration space used 5 or 18 channels.

The i n i t i a l

calcula-

Recently, i t has been

shown2 that i t is necessary to use 34 channels to obtain a converged r e s u l t for the case where there is a three-body force as well as the usual two-body force. 2. CONFIGURATIONSPACE FADDEEV EQUATIONS The Faddeev equations in c o n f i g u r a t i o n space were f i r s t f o r the three-body scattering problem.

derived by Noyes3

These equations were f i r s t

used for the

bound-state problem by the Grenoble group 4 who used f i n i t e difference techniques to solve the equations.

To derive the equations one introduces the

Jacobi variables for three equal mass p a r t i c l e s with position vectors ~ i , ~ j , and ~k : ~
;

Yi = ~ ( r j + ~k) - ~i

,

*This work was supported, in part, by the U.S. Department of Energy. 0 3 7 5 - 9 4 7 4 / 8 7 / $ 0 3 . 5 0 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

(1)

G.L. Payne / The Lanczos algorithm and spline expansions

140c

where i , j ,

k imply c y c l i c permutation.

Using these v a r i a b l e s the Schrodinger

wave f u n c t i o n is w r i t t e n as the sum of the three Faddeev amplitudes: +

= (1 + P- + P+)@(~z,Yl) =

e(~1,Yl) + ~(xZ,Y2) + ~(X3 ,Y3)

= ~i

+ ~b2 + ~3

where P- and P+ are the c y c l i c permutation operators.

(2)

,

For a system with a

Hamiltonian of the form: H = T + V(~I) + V(~2) + V(;~3) = T + VI + V2 + V3

,

(3)

where T is the k i n e t i c energy operator and V(~i) is the two-body p o t e n t i a l acting between p a r t i c l e s j and k, the Schrodinger equation can be replaced by the three coupled Faddeev equations (T - E)~ i + V ( ~ i ) ( ~ i + ~j + ~k) = 0

(4)

For the case of i d e n t i c a l p a r t i c l e s a l l three Faddeev amplitudes have the same f u n c t i o n a l form, and i t

is only necessary to solve one of the equations.

A f t e r expanding the Faddeev amplitude in a complete set of channel f u n c t i o n s , one obtains a set of coupled e l l i p t i c

partial differential

equations. The

standard numerical methods f o r solving these equations f o r the eigenvalue, E, and the e i g e n f u n c t i o n lead to p r o h i b i t i v e matrix equations to solve.

The size

of the numerical c a l c u l a t i o n can be considerably reduced by using a scheme s i m i l a r to one used in momentum space c a l c u l a t i o n s . 5

This method consists of

r e w r i t i n g the Faddeev equation in the form: [E - (T + V1)]~l : XVl(@2 + @3)

(5)

Now one assumes a value f o r E and solves t h i s equation f o r the eigenvalue x. The value of E is varied u n t i l

the eigenvalue ~ is u n i t y .

Using the Lanczos

method described below, one must i n v e r t a matrix corresponding to the l e f t - h a n d side of equation (5).

Since the operator (T + V - E) couples at most two chan-

nels ( f o r the case with a tensor f o r c e ) , the r e s u l t i n g matrix equation never i n v o l v e s more than two channels at a time.

Therefore, the computer time

required to solve the equations varies l i n e a r l y with the number of channels, r a t h e r than as the cube. 3. SPLINE EXPANSION To solve equation (5), we f i r s t

expand the Faddeev amplitude in a complete

set of channel functions

,Yl) = Z

*~(xI'Yl) xlYl

(6)

G.L. Payne / The Lanczos algorithm and spline expansions

141 c

where ~ l a b e l s the various angular momentum and s p i n - i s o s p i n f u n c t i o n s used to d e s c r i b e the three-body system.

A l s o , we have introduced the reduced channel

f u n c t i o n s @ a ( x l , y l ) which have the boundary c o n d i t i o n s t h a t the f u n c t i o n s are zero when xz or Yl is zero.

M u l t i p l y i n g t h i s equation by x l Y l and t a k i n g the

i n n e r product w i t h <~I gives a set of coupled e l l i p t i c

differential

equations.

In a d d i t i o n , we i n t r o d u c e the h y p e r s p h e r i c a l v a r i a b l e s x l = p cos o

,

This change of v a r i a b l e s s i m p l i f i e s

Yl = ~

P sin e

(7)

the p r o j e c t i o n i n t e g r a l s on the r i g h t - h a n d

side of the e q u a t i o n , as they depend only on the v a r i a b l e e and not p.

Thus,

we o b t a i n the set of coupled equations

(E - T)@a(p,O) - Z V B(p COS 8)@8(p,8) O+ V¥(p COS O) f _ Kyg(e,o')@6(p,O)de' By o

,

= ~ Z

(8)

where V~B(p cos e) is the projection of the two-body potential onto the channel states, and the kernel Ky~(e,o'), which depends upon the angular momentum and isospin variables, can be evaluated by standard techniques. 6 To solve this equation we expand the channel functions using a basis set of bicubic Hermite s p l i n e s on a r e c t a n g u l a r g r i d in the p-O coordinates. We choose an expansion o f the form: M

@a(P,O) = [ Z

m e-KP Z amn Sm(P)Sn(O)] p~/Z

N

m=l n=l

,

(g)

where ~ K 2 / M = -E, and the Sm and s n are the cubic Hermite s p l i n e s . t o have a smoother f u n c t i o n to be f i t

out the asymptotic behavior of the bound-state f u n c t i o n . t i o n s f o r the d i f f e r e n t i a l

The boundary c o n d i -

equation are i n c o r p o r a t e d i n t o the basis f u n c t i o n s .

The Hermite s p l i n e s are l o c a l f u n c t i o n s w i t h a continuous f i r s t These f u n c t i o n s are defined by d i v i d i n g the region to be f i t intervals,

In o r d e r

by the s p l i n e expansion, we have f a c t o r e d

derivative.

i n t o a number of

where the ends o f these i n t e r v a l s are c a l l e d the b r e a k p o i n t s .

At

each b r e a k p o i n t t h e r e are two s p l i n e s which are nonzero in the a d j a c e n t i n t e r vals.

In Figure I we i l l u s t r a t e

the s t r u c t u r e of these f u n c t i o n s f o r the case

w i t h f i v e i n t e r v a l s and the boundary c o n d i t i o n t h a t the f u n c t i o n be zero at each end of the r e g i o n .

Since the s p l i n e s are l o c a l f u n c t i o n s , t h e r e are at

most f o u r nonzero s p l i n e s on each i n t e r v a l

f o r the cubic s p l i n e s .

We use the

method of orthogonal c o l l o c a t i o n 7 t o determine the expansion c o e f f i c i e n t s

asmn.

G.L Payne / The Lanezos algorithm and spline expansions

142c

The collocation procedure requires that one determine the values of as mn for which the d i f f e r e n t i a l equation

c- c.B6 584

HERMITESPLINES l

-

-

is s a t i s f i e d at M d i s t i n c t values of Pi and N d i s t i n c t values of ej. I f

-

one chooses these points to be the ×

5

two point Gauss quadrature points on

6 7

-

each i n t e r v a l , the procedure is known as orthogonal collocation. The col-

-

8

location points for the example shown

9

in Figure i are marked by x's along

IO x

x

O

I

5

x

x~

x

x

IO

; x

15

x

I x

20

25

x

the x-axis. The orthogonal collocation tnethods leads to the set of simultaneous equations for the expansion coefficients a(% mn

FIGURE 1 Hermite spline functions for f i v e integral s. M

N

~ (~ +1) m

L (L +1)

m

m=IZ n=l ~ {Sm(Pi)Sn(°~)J + P-~I [¼_ c ° s ~ e J

(%

- ~

m

]

Sm(Pi)Sn(ej)

i (% + --2 Pi Sm(Pi)Sn(O~) J - 2~ s'm(Pi)Sn(Sj)}amn M

N

Z

Z

M

N

Z

Z

Nc

m=l n=l ~=I

:x

[v 6(p i cos ej)Sm(Pi)Sm(ej)]a~n 8+

Nc

m=l n=l ~=1

{~¥ v y(p i cos e j ) [ f 0 1 Ky6(ej,e')Sn(O')de']Sm(Pi)}aBmn "(i0)

In these equations ~m is the orbital angular momentum of the interacting pair 2-3, L(% is the orbital angular momentum of the spectator p a r t i c l e i , and vm6 is (M/~Y~)Vm6. This equation can be written symbolically as

M N M c Aijm,mn6 a~ ~ ~c Bijm,mn~ a6 mn Z N Z mn = ~ ~ m=l n=l 6=i m=l n=l B=I

(11)

which can be written as the generalized matrix eigenvalue problem Aa = ~Ba

,

(12)

where the eigenvector is the array of coefficients for the spline expansion. For a r e a l i s t i c calculation the size of this matrix can be quite large. A typical 34-channel (Nc = 34) calculation may require 30 splines for the p

G.L. Payne / The Lanczos algorithm and spline expansions variable (M = 30) and 30 splines for the 9 variable (N = 30).

143c

Thus, the total

number of expansion coefficients would be 34 x 30 x 30 or 30,600 unknowns. Consequently, this matrix is too large to be solved by straightforward techniques, and we use a variation of the Lanczos algorithm to find the solution. 4. LANCZOS METHOD Since we wish t o f i n d the e i g e n v a l u e s c l o s e s t t o u n i t y , e q u a t i o n (12) in the form A- i

B a = (I/~)a

we f i r s t

rewrite

or

H a = A a To s o l v e t h i s

(13)

e q u a t i o n , we use a m o d i f i e d Lanczos a l g o r i t h m t o g e n e r a t e a small

basis set which can be s o l v e d by standard t e c h n i q u e s .

Since the m a t r i c e s A and

B are not symmetric, we have t o generate a b i o r t h o g o n a l basis s e t . done using the method o f Saad. 8

Starting

w i t h an i n i t i a l

This i s

vector ai,

e r a t e the basis set a i as w e l l as the b i o r t h o g o n a l basis set bi f o r . . o ~.

These basis v e c t o r s are generated by f i r s t

we geni = 1,2,

using the r e c u r s i o n

relations N

a i + l = Hai - ~ i a i

(14)

- ~iai-1

and bi+l

= HTbi - ~ i b i

where a i = b~Hai , Bi = b T a i / Y i , basis v e c t o r s .

and ¥i =

" Yibi-1

'

(15)

I T iL I/2,

t o o b t a i n the unnormalized

The n o r m a l i z e d v e c t o r s are d e f i n e d as a i = a i / ~ i

bi = ~ i / B i . Thus, one can generate u v e c t o r s ai and t h e i r

and

b i o r t h o g o n a l v e c t o r s b i such

t h a t bT aj = a i j . Using t h i s

small basis s e t , we use t h e expansion V

a = Z ~iai

(16)

i=1

Substituting this expansion into equation (13) and taking the inner product with bj leads to the matrix equation V

v L

~ibTHai = Z Hji~i = A~j

i =1

The e i g e n v a l u e s and e i g e n v e c t o r s o f t h i s t o the l a r g e m a t r i x e q u a t i o n .

(17)

i=1 m a t r i x e q u a t i o n are a p p r o x i m a t i o n s

In p r a c t i c e one i n c r e a s e s u u n t i l

A i s l e s s than the d e s i r e d accuracy o f s i x s i g n i f i c a n t

figures.

the change in

G.L. Payne / The Lanezos algorithm and spline expansions

144c

The most d i f f i c u l t

numerical procedure of this method is solving the recur-

sion r e l a t i o n s in equations (14) and (15). sion of the matrix A.

This requires the numerical inver-

However, the choice of the equations to be solved leads

to a reduced matrix whose blocks never have more than two channels, and each of these blocks can be inverted separately.

F i n a l l y , the inversion of these

blocks is considerably s i m p l i f i e d because of the local nature of the splines. Since there are at most four nonzero splines on each i n t e r v a l , the r es u lt in g matrix is a sparse banded matrix.

By the proper ordering of the equations and

the unknowns, one obtains a block t r i d i a g o n a l matrix whose blocks are banded matrices.

This matrix can be inverted by using an algorithm which works on a

few blocks at a time.

Thus, the amount of computer memory needed does not

become excessive, and the inversion of A can be performed in an e f f i c i e n t manner. 6. CONCLUSIONS The use of spline expansions provides an e f f i c i e n t method for obtaining accurate solutions of the configuration-space Faddeev equations. i s t i c c a l c u l a t i o n one must use a modified form of these equations.

For a r e a l The Lanczos

algorithm f o r nonsymmetric matrics can be used to find the solution f o r the bound-state problem. REFERENCES i) G.L. Payne, J.L. F r i a r , B.F. Gibson, I.R. Afnan, Phys. Rev. C 22 (1980) 823. 2) C.R. Chen, G°L. Payne, JoL. F r i a r , B.F. Gibson, Phys. Rev. C 33 (1986) 1740. 3) H.P. Noyes, in Three Body Problem in Nuclear and P a r t i c l e Physics, ads. J.S.C. McKee and P.M. Rolph (North-Holland, Amsterdam, 1970), p. 2. 4) S.P. Merkuriev, C. Gignoux, and A. Laverne, Ann. Phys. (NY) 99 (1976) 30. 5) Ch. Hajduk and P.U. Sauer, Nucl. Phys. A369 (1981) 321. 6) C.R. Chen, Effects of three-nucleon forces on the trinucleon system, Ph.D. d i s s e r t a t i o n , University of lowa, 1985. 7) P.M. Prenter, Splines and V a r i a t i o n a l Methods (Wiley, New York, 1975). 8) Y. Saad, SIAM J. Num. Anal. 19 (1982) 485.