NOVEMBER/DECEMBER 2005 (Vol. 7, No. 6)

An Efficient Numerical Spectral Method for Solving the Schröger Equation: Appendix 1: Matlab Program for Scattering in a Morse Potential

This appendix contains all the Matlab m-files needed to calculate the scattering of a particle in a 1D potential for a partial wave with L = 0, using the numerical spectral integral equation (S-IEM) method we describe in the main article text.1 The potential included is an inverted Morse potential that has a negative (attractive) valley at small distances and a repulsive barrier at intermediate distances, and that decays exponentially to zero for large distances. This potential supports various resonances for different incident energies. We calculate the phase shifts for various energies contained within one such resonance and compare them with the analytic results, as we describe in the main text,1 to measure our method's accuracy. We list here the main program, IEM_Morse_res.m, and the various function-subroutines that it calls. The function-subroutine meshiem calculates a partition's size for a given input tolerance parameter; this subroutine in turn calls mapxtor and vfg. The potential and the Green's functions (sin(kr) and cos(kr)) are calculated by function-subroutine vfg. The potential is calculated by calling an appropriate subroutine (morse in the present case) that different potential subroutines can replace. The function-subroutine YZ calculates the local functions Y and Z; it calls mdiag and vfg.

## Reference

• 1. G.H. Rawitscher I. Koltracht, "An Efficient Numerical Spectral Method for Solving the Schrödinger Equation," Computing in Science & Eng., vol. 7, no. 6, 2005, pp. 58-66.

Main program IEM_Morse_res.m

% IEM_Morse_res.m is a MATLAB program that
% solves a one dimensional, one channel, positive energy (scattering)
% Lippmann-Schwinger eq. using The Chebyshev spectral expansion method,
%  for a resonance scattering in a inverted Morse potential with a
% barrier
%   Subroutines needed are C_CM1, SL_SR, meshiem, YZ, vfg, morse,
% mapxtor, and mdiag.
% meshiem calls vfg and mapxtor
% vfg calls a potential subroutine, that in this case is: morse
% YZ calls vfg and mdiag
%

clear all
tic % starts here to count the elapsed time
N=16 %number of Chebyshev mesh points in each partition, as described
% more precisely in the subroutine C_CM1 below
%T=100
nloop=10
bstart= 0.2
bstart= 0 % leftmost radial distance of the first partition
Amorse=4; REmorse=4;amorse=1/0.33; % parameters for the Morse potential
tol=1.e-12
rcut=bstart % the potential is replaced by a constant for r<rcut
%k=0.10
%K = 41.47 % K=(hbar)2/(reduced mass of nucleon)
K=1
%if K=1, this means that the potentials and energies have
% already been divided by K.
% The resulting dimensions of the energy and potential is inverse
% length squared, and the wave number k is one over a length. For
% this morse pot'l case, the unit of length is femto-meter

% anal50 is the set of L=0 phase shifts obtained analytically at a a
%distance of 50 fm. See Am. J. of Phys. vol 70, p. 935 (02).
% The first column is the value of k, the second is the phase shift
% in units of pi.
% anal100 is the same as anal50 for a matching radius of 100 fm
anal50=[
1.507100000000000e+000 -4.568483126920000e-001
1.507110000000000e+000 -4.533212648736000e-001
1.507120000000000e+000 -4.480714272385000e-001
1.507130000000000e+000 -4.394398160329000e-001
1.507140000000000e+000 -4.226587240332000e-001
1.507150000000000e+000 -3.765904018859000e-001
1.507160000000000e+000 -4.635282299884000e-002
1.507170000000000e+000 4.113560735221000e-001
1.507180000000000e+000 4.687007391308000e-001
1.507190000000000e+000 4.877779485905000e-001
1.507200000000000e+000 4.972300508370000e-001
1.507210000000000e+000 -4.971373180371000e-001
1.507220000000000e+000 -4.934015231054000e-001];
anal100=[
1.507100000000000e+000 -4.56851148617700e-001
1.507110000000000e+000 -4.53324094250100e-001
1.507120000000000e+000 -4.48074247045200e-001
1.507130000000000e+000 -4.39442620401400e-001
1.507140000000000e+000 -4.22661499186800e-001
1.507150000000000e+000 -3.76593103155600e-001
1.507160000000000e+000 -4.63555541276000e-002
1.507170000000000e+000 4.11353023527700e-001
1.507180000000000e+000 4.68697775626900e-001
1.507190000000000e+000 4.87775017641500e-001
1.507200000000000e+000 4.97227136582100e-001];
%[a,b] = size(anal100);
[a,b] = size(anal100); %this determines the number of k-values to be
% used in various do-loops later on
for n=1:a;
tandl(n)=1;
kk(n)=0;
end

[C,CM1,xz]=C_CM1(N); % these matrices map the vector of the function into
% the vector of the expansion coefficients (all column vectors)
[SL,SR]=SL_SR(N); % these are the right and left indefinite quadrature
% matrices in the Chebyshev space, described below
% b = SL(a).
% b is the column vector of the Ch. coeffs
% of the integral from -1 to x of f(x)dx,
% a is the column vector of the Ch. coeffs of the expansion of f.
%  ditto for SR, for the integr from x to +1 SRCM1=SR*CM1;
%these are matrices that occur again and again in the calcu-
%lation of the local functions Y and Z in each partition.

SLCM1=SL*CM1;
nkmax=8;
nkmin=5
%for nk=1:13;
%for nk=nkmin:nkmax;
for nk=1:a
b1=bstart;
k=1.5071+(nk-1)*0.00001;
%determines the k-values for this particular resonance calculation
kk(nk)=k;
%K=1/0.04783258
% np = 1, 2, ...npm is the index thatlabels the partitions.
%    Partition np starts with bp(np-1) and goes to bp(np),
% also called b1 and b2.
% partition 1 starts with b1 prescribed above.
% These partition numbers are stored in the array mp,
% the expansion coeffs aY for Y(x) in partition np for n= 0, 1, ..N are %
contained in the matrix aYp, (NP1 by npm) with column
% np containing aY for partition np;
% ditto for Z, aZp; ditto for the overlap integrals, in matrix overlp.

% the method to establish the length of each partition is as follows.
% Starting from the left end b1, the length of the partition is first
%  estimated from the local wave number kbar as 3.2*pi/kbar,
%  i.e., about 10 Cheb. points per local wave length. For this
%    particular estimate of the size of the partition
%   the Cheb. Expansion coeffs. of the potential times the F-
%  component of the greens function (sin(kr)) is calculated with Cheb
% indices 0, 1, ...N. If the sum of the abs. values
% of the last three coeffs is larger than the
%  prescribed tolerance, the size of the partition is divided by two.
% The process is repeated until the accuracy criterion is satisfied,
% with a maximum of nloop repetions, and a first approximation to
% b2 is found.
%  This approximate value is then further refined by doing the full
%  calculation of the functions Y and Z in this particular partition,
% and if the tolerance criterion, based on the square root of the
% sum of the squares of the three last Cheb. coefficients
%  for both Y and Z, is met, the do-loop is exited.
% The results for the CH. Coeffs of Y and Z are stored in the matrices
% aYp and aZp, and the overlap integrals are stored in the matrix
% overlp. Each column corresponds to a particular value of np.
%
%  The last partition ends at bp(npm) = T, and usually is of a
% smaller length than what the tolerance parameter requires.

b2p=b1;
np=1;
mp(1)=1; errp(1)=0; bp(1)=0;
% j is the partition number

%display(' here it calls mesh')
format short e
while b2p < T;
%  subroutine meshiem finds the first approximation b2 to the end of
% the partition.
[b2,err]=meshiem(N,nloop,b1,tol,rcut,k,K,xz,C,CM1);
delta = b2-b1;
%b2;
%display('now enter YZ loop')
% subroutine YZ calculates the functions Y and Z in the partition for
% which it is called. This routine calls the subroutine vfg, that in turn
% calls the potential sub-routine.
% In order to change the potential, it has to be done by changing the
% potential subroutine called by vfg.
for i=1:nloop;
[aY,aZ,afv,overlaps,errYZ]=YZ(N,b1,b2,rcut,k,K,xz,SLCM1,SRCM1,C,CM1);
erYZ(i)=errYZ;
%display(' YZ loop errYZ')
%errYZ,delta
delt(i)=delta;
if errYZ > tol
delta=delta/2;
b2=b1+delta;
else break
end
end
%  the final value of b2 is now found so that both Y and Z satisfy the
% tolerance criterion.
bp(np)= b2;
[v,f,g]=vfg(b2,rcut,k);
% This gives the value of the potential at b2, needed later plotting
vp(np)=v;
b2p=b2;
b1=b2;
mp(np)=np;
%  The values of err, aY, aZ, overlaps are now stored in the respective
% matrices.
errp(np)=errYZ;
aYp(:,np)=aY;
aZp(:,np)=aZ;
overlp(:,np)= overlaps;
np=np+1;
%display(' YZ end loop delta YZerr ')
%[delt',erYZ']
%display(' overlap integrals fvy gvy, fvz, gvz')
%overlaps'
%display(' converged Cheb exp coeffs for fv Y and Z')
% [afv,aY, aZ]
end

% now the last partition going to b2=T will be constructed, and the
% results of the last partition (end of the "while" loop) will be
% overridden by the calculation for the last partition with b2=T.
% This is needed since the right hand radial value of the last
% partition may not coincide with the pre-assigned value T.

b1=bp(np-2);
b2=T;
bp(np-1)=T;
[v,f,g]=vfg(T,rcut,k);
vp(np-1)=v;
mp(np-1)=np-1;
[aY,aZ,afv,overlaps,errYZ]=YZ(N,b1,b2,rcut,k,K,xz,SLCM1,SRCM1,C,CM1);
npm=np-1;
%display(' overlap integrals for final partitn: fvy gvy, fvz, gvz')
%overlaps'
%display(' converged Cheb exp coeffs for fv Y and Z')
%[afv,aY, aZ]

errp(np-1)=errYZ;
aYp(:,np-1)=aY;
aZp(:,np-1)=aZ;
overlp(:,np-1)= overlaps;
%format short g
[v,f,g]=vfg(rcut,rcut,k);
%display('v at b0')
% v;,rcut;

%display (' np b2(np) v(b2) error(np) after convergence')
%[mp',bp',vp',errp']
format short e
%display(' converged Cheb exp coeffs for Y in all partitions')
%aYp
%display(' converged Cheb exp coeffs for Z in all partitions')
%aZp
%display(' converged overlaps in all partitions')
%overlp

% *******
% now the values of the overlaps and the aY, and the aZ are stored from
% np= 1 to np = npm,
%  The values of the coeffs A, B for each partition are obtained in
% terms of the previous partition. The row [A,B]' is called alpha
% Omega(np+1)*alpha(np+1) = gamma(np)* alpha(np)= beta(np)
% alpha(np+1)=omega(np+1)\beta(np)
% beta (np+1)= gamma(np+1)*alpha(np+1), etc

%  start a do-loop from np=2 to np=npm, to carry out the above,
% and first
% obtain beta(1) by assuming A(1)=1, B(1)=0, and beta(1)= [1,-FVY]

np=1;
%below some of the overlap integrals are stored
FVY=overlp(1,np); GVY= overlp(2,np); FVZ=overlp(3,np); GVZ=overlp(4,np);
%beta=[1;-FVY]; A(1) = 1; B(1)=0;
beta=[-FVY;1]; A(1) = 1; B(1)=0;

for np=2:npm;
FVY=overlp(1,np); GVY= overlp(2,np); FVZ=overlp(3,np); GVZ=overlp(4,np);
%omega =[1-GVY,-GVZ;0,1];
omega =[0,1;1-GVY,-GVZ];
alpha =omega\beta;
% this value of beta is carried over from the previous cycle
A(np)=alpha(1); B(np)=alpha(2);
%BoA(np)=B(np)/A(np);
%gamma = [1,0;-FVY,1-FVZ];
gamma = [-FVY,1-FVZ;1,0];
beta = gamma*alpha;
tand(np)=beta(1)/beta(2);
%ctand(nk)=beta(2)/beta(1);
%delta(nk)=(atan2(beta(1),beta(2)))/pi;
%delta(nk)=(atan(tand(nk)))/pi;

end
%this is the end of the partition loop
tandl(nk)=tand(npm);
end %This is the end of the k-loop
display (' np b2(np) v(b2) error(np) after convergence')
[mp',bp',vp',errp']

% Delta is the phase shift
delta=(atan(tandl))/pi;

format long e
% now
% print out the errors, compared to the values stored in the analytic
% values stored in anal50 or anal100 The corresponding matching point
% T was inserted manually at the start of the program , line 13. Also
% the value of the tolerance that determines the size of a partition
% so as to achieve a desired accuracy for the local functions Y and Z
% is read in at the start of the program

diff=abs(delta'-anal100(:,2));
%diff=abs(delta'-anal100(:,2));
logdiff=log10(diff);
%display('coeffs A and B for each partition, and tan(delta)')
display(' k, delta/pi')
format short e
%[tol,N,T,bstart]
[Amorse,REmorse, amorse]
%format long e
kk=(kk-1.5071)*1.e5;
[kk',diff,logdiff]
% tic and toc determine the time elapsed in between
toc
% This is the end of main program, IEM_Morse_res.m,

file C_CM1.m
function [C,CM1,xz]=C_CM1(N)

%f(x)=sum(1:NP1)a(n) T(j=N-1,x), j is the max power of x in T(j,x)
% here there is no a_0, but a(1)= a_0/2
%calculate the coefficients a(j);
% j=0 to N from the Clenshaw Curtis alg.
% column vector a = CM1 * (column vector of the f's)
% column vector f at Ch. zeros of T(j=N+1) = C * (column vector a)
% xz is the column vector of the zeros of T(j=N+1)
% format short g
% format short e

% N is the highest order Chebyshev included in the expansion of f(x)

NP1=N+1;
format long e
%calc the zeros of Chebyshev of order NP1. Here m=j+1
for m=1:NP1
theta(m)=(-1+2*m)*pi/(2*NP1);
y(m)=cos(theta(m));
%f(m)=exp(y(m));
end
xz=y'
%disp( 'the values of x and f(x)')
%[y',f'];
%Now calculate the matrix CM1(NP1,NP1) (i.e, the inverse of C)
% this requires the matrix
% T(i,j)=T(0, x0,x1,..xN);T(1,x0,x1,..xN);

for i=1:NP1
for j=1:NP1
T(i,j)=(cos((i-1)*theta(j)));
if i==1
A(i,j)=1/2;
else A(i,j)=T(i,j);
end
end
end
C=T';
CM1=A*(2/NP1);
%O=CM1*C;
%O;
%test if CM1 is the inverse of C.
% yes, O= unity matrix
return

file SL_SR.m
function [SL,SR]=SL_SR(N)
% N is the highest order Chebyshev included in the expansion of f(x)
%b=SL*a, b is the column vector of the coefficients of the expansion
%into Chb's of integral from -1 to x of f(x) dx
%c=SR*a, c is the column vector of the coefficients of the expansion
%into Chb's of integral from x to 1 of f(x) dx

%first construct the integr routines SL=M11*U11
NP1=N+1
for i=1:NP1
for j=1:NP1
M11(i,j)=0;
U11(i,j)=0;
B(i,j)=0;
if i==j
M11(i,j)=1;
M22(i,j)=-1;
end
end

end

% now fill the top row
i=1; M11(i,2)=1; M22(1,1)=1; M22(1,2)=1;
for j= 3:NP1
M11(i,j)= M11(i,j-1)*(-1);
M22(i,j)= 1;
end
%disp (' M')
%M11;
%M22;
%now construct U11
aux=0;
for i=2:NP1-1
aux=aux+2;
U11(i,i+1)=-1/aux;
U11(i,i-1)=1/aux;
end
U11(2,1)=1;
U11(NP1,NP1-1)=1/(aux+2);
%disp('U; SL=M*U ')
%U11;
SL=M11*U11;
SR=M22*U11;
% These are the NP1*NP1 integration matrices
return
% file C_CM1.m ends here

file maprtox.m
function x=maprtox(b1,b2,r)

% x is the column vector of NP1 numbers in (-1,+1)
% r is the corresponding column vector stretching from (b1, b2)

auxm = (b2-b1)/2;
auxp = (b2+b1)/2;
auxr = auxp/auxm;
x = (1/auxm)*r-auxr;
return

file mapxtor.m
function r=mapxtor(b1,b2,x)

% x is the column vector of NP1 numbers in (-1,+1)
% r is the corresponding column vector stretching from (b1, b2)

auxm = (b2-b1)/2;
auxp = (b2+b1)/2;
auxr = auxp/auxm;
r = auxm*(x+auxr);
return

file mdiag.m
function D=mdiag(x)

% x is a column vector of MP1 elements to be placed on the diagonal of
% matrix D

[NP1,p] = size(x);
A = zeros(NP1,NP1);
for i= 1:NP1
A(i,i)=x(i);
end
D=A;
return

file meshiem.m

function [b2,err]=meshiem(N,nloop,b1,tol,rcut,k,K,xz,C,CM1)

% calculates the value of the b2, the right limit of the partition that
% starts at the left with b1,
% based on the magnitudes of the last few coeffs of the Chebyshev
% expansion of v.*f. If they are larger than tol then the
% partition is cut in half
[v,f,g] =vfg(b1,rcut,k);
%display('v(b1)')
%v
v=v/K;

kloc=sqrt(abs(k*k-v));
if kloc < k
kloc=k;
end

delta=10/kloc;

for i=1:nloop
del(i)=delta;
b2t = b1+delta;

%b2t is the trial partition right hand limit

r    = mapxtor(b1,b2t,xz);
% maps the vector xz into the radial distance vector r
[v,f,g]=vfg(r,rcut,k);
v=v/K; %now v stands for VBAR
av=CM1*(v);
af=CM1*(f);
afv=CM1*(f.*v);
[afv];
aux=afv(N+1)^2 + afv(N)^2 + afv(N-1)^2;
err(i)=sqrt(aux);
b2=b2t;
if err(i) > tol
delta=delta/2;
else
break
end
end
[del',err'];
return
% end of meshiem

file morse.m
function v=morse(rr);
VBAR = 4;
alpha = 1/0.3 ;
x = (rr-4)/alpha;
y=exp(-x);
z= (y-2);
v=-VBAR*(y.*z);
return;

file vfg.m
function [v,f,g]=vfg(r,RCUT,k)

% given the range of points in r from b1 to b2,
% then calculate the corresponding potential values
[NP1,p]=size(r);
for i=1:NP1

if(r(i)>=RCUT)
rr(i)=RCUT;
else
rr(i)=r(i);
end
end
vv=morse(rr);
% For other potentials replace the above subroutine
%   by one that corresponds to a different potential
% for example vv=MalfTj(rr);

v=vv';
%Dv=mdiag(v) this diag matrix is not needed here;
f=sin(k*r);%Df=mdiag(f);
g=cos(k*r);%Dg=mdiag(g);
return
% end of file vfg.m

file YZ.m
function [aY,aZ,afv,OVERLAPS,errYZ]=YZ(N,b1,b2,rcut,k,K,xz,SLCM1,SRCM1,C,CM1)
% solves the integral eq. for the Chebyshev expansion coefficients of
% the local functions Y and Z in the interval b1 <r < b2,
% N is the order of the last Chebyshev used in the expansion of Y and Z
% k = the wave number for the Greens functions
% K = h^2/(2M) in MeV fm^2.
% The potentials called by the subroutines are in MeV,
% hence VBAR=V/K and EBAR = k^2 = E/K are in fm^(-2)
%
NP1=N+1;
FACTOR=(b2-b1)/(2*k);
% (b2-b1)/2 converts the integration over dr into integr over dx
%  dividing by K converts the potential calculated in MeV to the
% units of fm^(-2) and k is the wave number in fm^(-1). Hence, after
%  the multiplication by the FACTOR, the integrals over products of
%   functions (f, g, or X, or Y) times potential times dx
% have no dimension.

I=eye(NP1);
% I is the unit matrix NP1xNP1
% xz is the column vector of the zeros of the Ch pol of order NP1=N+1
%  af = CM1(f(xz)) is the column vector of the Ch. Coeffs of the
% expansion of f.
% Here af(1) is the coeff of T(0,x) , i.e., =af0/2; f(xz)=C(af)
% is the column vector of the values of f(x) at the zeros xz
format short e
%[SL,SR]=SL_SR(N);
% b = SL(a).
% b is the column vector of the Ch. coeffs
% of the integral from -1 to x of f(x) dx,
% and a b is the column vector of the Ch. coeffs of the expansion of f.
%   ditto for SR, for the integr from x to +1

r   = mapxtor(b1,b2,xz);
% maps the vector xz into the radial distance vector r
[v,f,g]=vfg(r,rcut,k);
v=v/K; %now v stands for VBAR
Dv=mdiag(v); %Dv is the diagonal matrix that has the vector v at its diagonal
Df=mdiag(f);
Dg=mdiag(g);
%display(' x r v f g');
[xz,r,v,f,g];
av=CM1*(v);
af=CM1*(f);
afv=CM1*(f.*v);
ag=CM1*(g);
%SLCM1=SL*CM1;
%SRCM1=SR*CM1;

MG = CM1*(Df*C*SRCM1*Dg + Dg*C*SLCM1*Df)*Dv*C;
M=(I+FACTOR*MG);
aY = M\af;
aZ = M\ag;
aux=aY(N+1)^2 + aY(N)^2 + aY(N-1)^2;
aux=aux + aZ(N+1)^2 + aZ(N)^2 + aZ(N-1)^2;
errYZ=sqrt(aux);

%display(' Cheb exp coeffs for v v*f Y Z f g');
[av, afv,aY,aZ,af, ag];

Y=C*aY;
Z=C*aZ;
%Now calculate the overlap integrals
% first calculate the integrals from -1 to x of g*v*z
bgvz=SLCM1*(g.*v.*Z);
bfvz=SLCM1*(f.*v.*Z);
bgvy=SLCM1*(g.*v.*Y);
bfvy=SLCM1*(f.*v.*Y);
% next sum the coefficients of the Cheb expansion of the integral,
% in order to get the integral from -1 to +1 or from b1 to b2.
% This is so because when x = 1, all chebyshevs are equal to 1.
GVZ=FACTOR*sum(bgvz);
FVZ=FACTOR*sum(bfvz);
GVY=FACTOR*sum(bgvy);
FVY=FACTOR*sum(bfvy);
OVERLAPS=[FVY;GVY;FVZ;GVZ];
Return
% this is the end of subroutine YZ,
% All programs are now listed above


Article Contents