Nuclear Physics A463 (1987) 139c  144c Northltolland, 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 threebody 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 boundstate 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 threebody force as well as the usual twobody 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 threebody scattering problem.
derived by Noyes3
These equations were f i r s t
used for the
boundstate 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. (NorthHolland 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 twobody 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 threebody 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 twobody 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 pO coordinates. We choose an expansion o f the form: M
@a(P,O) = [ Z
m eKP 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 boundstate 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 xaxis. 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 23, 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 34channel (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)
 ~iai1
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 =
" Yibi1
'
(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 configurationspace 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 boundstate 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 (NorthHolland, 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 threenucleon 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.