The Laguerre pseudospectral method for the radial Schrodinger equation

13 Pages·2014·0.42 MB·English

by  H. Alıcı

AppliedNumericalMathematics87(2015)87–99
ContentslistsavailableatScienceDirect
Applied Numerical Mathematics
www.elsevier.com/locate/apnum
The Laguerre pseudospectral method for the radial
Schrödinger equation
H. Alıcıa,∗, H. Tas¸elib
aDepartment of Mathematics, Harran University, 63290, S¸anlıurfa, Turkey
bDepartment of Mathematics, Middle East Technical University, 06800, Ankara, Turkey
a r t i c l e i n f o a b s t r a c t
Article history: By transforming dependent and independent variables, radial Schrödinger equation is
Received 21 April 2014 converted into a form resembling the Laguerre differential equation. Therefore, energy
Received in revised form 5 August 2014 eigenvalues and wavefunctions of M-dimensional radial Schrödinger equation with a wide
AAvccaeilpatbelde o3n Slienpet e1m6b Seerp 2te0m14ber 2014 rmanetgheo dosf . iCsootmroppairci spoont ewntitiahl st haer er eosbutlatisn efrdo mnu mliteerriactaulrlye bshyo wussi ntgh aLta gtuheer rme epthseoudd oiss pheicgthrlayl
competitive.
Keywords:
Laguerre pseudospectral methods ©2014 IMACS. Published by Elsevier B.V. Allrightsreserved.
Radial Schrödinger equation
Quantum mechanical potentials
1. Introduction
In this article, we consider the radial part of the Schrödinger equation in M dimensions
(cid:2) (cid:3)
− d2 − M−1 d + (cid:2)((cid:2)+M−2) +V(r) R(M)(r)=E(M)R(M)(r), R(M)(r)∈L2(0,∞) (1)
dr2 r dr r2 n,(cid:2) n,(cid:2) n,(cid:2) n,(cid:2)
(cid:4)
with a variety of quantum mechanical potentials of the form V( x2+x2+...+x2 ) = V(r), where r denotes the
(cid:5) 1 2 M
M-dimensional spherical coordinates such that r2 = iM=1x2i. In (1), L2(0, ∞) is the Hilbert space of square integrable
functions on the half line, n, (cid:2) =0, 1, ... stand, respectively, for the radial and angular quantum numbers of the energy
eigenvalues E and the corresponding wavefunctions R(r). M-dimensional radial Schrödinger equation has been the subject
of many computational methods. The most commonly studied form of (1) is the three-dimensional case with M=3
(cid:2) (cid:3)
d2 (cid:2)((cid:2)+1)
− + +V(r) Ψ(r)=EΨ(r) (2)
dr2 r2
in which the first derivative term is removed by making use of the transformation Ψ(r) =rR(r) on the dependent variable.
It is interesting to notice that the last equation can still be derived from (1) without any transformation. To be specific,
in (1), regarding M as a parameter with the value M=1 and replacing (cid:2) with (cid:2) +1 we obtain (2). Therefore the link
between the eigenpairs {En,(cid:2), Ψn,(cid:2)(r)} and {En(M,(cid:2)), Rn(M,(cid:2))(r)} of (2)and (1), respectively, is given by the relations
En,(cid:2)=En(1,(cid:2))+1, Ψn,(cid:2)(r)=Rn(1,(cid:2))+1(r). (3)
* Corresponding author.
E-mail addresses:[email protected](H.Alıcı), [email protected](H.Tas¸eli).
http://dx.doi.org/10.1016/j.apnum.2014.09.001
0168-9274/©2014 IMACS. Published by Elsevier B.V. Allrightsreserved.
88 H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99
Various approximation methods have been proposed for computing the eigenvalues of this problem by many authors.
Among them we may recall variational methods [32–34,36,37,39], Hill determinant method[7], pseudospectral methods[2],
constant perturbation methods [17,20–22], Prüfer transformation followed by a shooting procedure [4,3], Pekeris-type ap-
proximations [12,15,16,18,25,41], asymptotic iteration technique [5] and Hamiltonian hierarchy picture [10].
In particular, in [2], we have studied the Schrödinger equation by an equation of the form
σ(ξ)y(cid:5)(cid:5)+τ(ξ)y(cid:5)+Q(ξ)y=−λy, (4)
called the equation of hypergeometric type with a perturbation (EHTP), by means of transformations on both independent
and dependent variables. Here σ(ξ) and τ(ξ) are polynomials of degrees at most two and one, respectively, and λ is a
parameter. In fact, in the special case of Q(ξ) ≡0, (4)reduces to the well known equation of the hypergeometric type whose
suitable solutions are the classical orthogonal polynomials. Accordingly, in [2], we have converted the radial Schrödinger
equation in (1) into an EHTP
(cid:6) (cid:7) (cid:6) (cid:7)
(cid:8) (cid:9) (cid:10) (cid:11)(cid:12)
ξy(cid:5)(cid:5)+ (cid:2)+ 1M−ξ y(cid:5)+ 1 ξ−c−2V c−1 ξ y= 2(cid:2)+M− 1c−2E y (5)
2 4 4
first by introducing the scaled quadratic variable ξ=(cr)2 where c is a positive constant and then proposing a solution of
the type R(ξ) =ξ(cid:2)/2e−ξ/2y(ξ) satisfying the asymptotic boundary condition at infinity and the regularity condition at the
origin. However, the quadratic transformation, which was necessary to convert (1) into an EHTP, does not seem suitable
for potentials that contain both odd and even powers of the variable r. For the sake of a form that is well suited for all
potentials, we have to sacrifice the EHTP. Therefore, in this study, instead of the quadratic transformation we start with a
more flexible one
ξ=(cr)α, α,c>0, ξ∈(0,∞) (6)
where the parameters α and c may be exploited to accelerate the convergence rate of the method. Then we suggest a
solution of type
R(ξ)=ξ(cid:2)/αe−ξ/2y(ξ) (7)
which transforms the radial Schrödinger equation (1) into
ξy(cid:5)(cid:5)+(γ +1−ξ)y(cid:5)+Q(ξ)y=λξα2−1y (8)
where
(cid:9) (cid:11)
Q =Q(ξ;α,c,γ)= 1ξ− 1(γ +1)− 1 ξα2−1V ξ1/α/c (9)
4 2 (αc)2
represents the modified potential, γ is the parameter
1
γ =γ(α,(cid:2),M)= (2(cid:2)+M−2), (10)
α
and
λ=λn(α,c,γ)=−(α1c)2En(M,(cid:2)) (11)
the rescaled energy eigenvalues.
The positive nonconstant term ξα2−1 on the right hand side of (8)can be seen as a weight function, and hence, the new
form may be considered as a weighted and perturbed Laguerre (WPL) equation. Now the EHTP in (5) is a special case of (8)
with α=2. Both (5) and (8) resemble the Laguerre differential equation
ξy(cid:5)(cid:5)+(γ +1−ξ)y(cid:5)=−ny, ξ∈(0,∞) (12)
where γ is a real parameter and n a non-negative integer. In the next section, we construct the pseudospectral formulation
of the WPL equation based on the associated Laguerre polynomials Lγ(ξ) which are the polynomial solutions of (12). Sec-
n
tion3introduces the numerical examples. Section4concerns with the implementation notes and the last section concludes
the paper with some remarks.
H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99 89
2. The Laguerre pseudospectral method (LPM) for the WPL equation
A pseudospectral method, also known as spectral collocation method, is based on the N-th degree polynomial interpola-
tion of a function y(ξ) denoted by PN(ξ),
(cid:13)N
PN(ξ)= (cid:2)n(ξ)yn, (13)
n=0
where the yn=y(ξn) are the actual values of y(ξ) at the specified nodes ξ =ξn for n =0, 1, ..., N [6,9,38]. Such a pseu-
dospectral scheme in which the N-th degree Lagrange polynomials
γ
(cid:2)n(ξ)= (ξ−ψξnN)+ψ1N(cid:5)(ξ+)1(ξn) = (ξ−ξn)[LdNdξ+L1γN(+ξ1)(ξ)]ξ=ξn, n=0,1,...,N (14)
are defined by the normalized
(cid:14)
ψn(ξ)=h1nLnγ(ξ), hn= Γ(n+nγ! +1), γ >−1 (15)
oorf LsγNta+n1d(aξr)d. ALpapgruoexrirme aptionlygn tohme isaolsl uLtinγo(nξs) oifs dcifafleleredn tai aLl PeMqu, awtiohnerse b yth Lea gnuoedrerse ξpnolyarneo mthiael sr eaarle, udsisutainllcyt naontd s tpaobsleit ifvoer rlaorogtes
N due to their wild behaviors at infinity, and hence, one usually works with the Laguerre functions ψn(ξ) =e−ξ/2Lnγ(ξ)/hn
instead. This situation is theoretically investigated, for example in [11,24,28,29], and it is shown that the Laguerre functions
have better stability properties than Laguerre polynomials. On the other hand, it is stated in [29] on p. 214 that the gener-
alized Laguerre polynomials are useful for the approximation of functions which decay at infinity. Therefore, in this study,
we continue with the normalized Laguerre polynomials in (15) since we search for the square integrable solutions of (1)
which behave suitably at the origin and vanish exponentially at infinity.
Notice that, y(ξn) =PN(ξn)at least at the nodes since the Lagrange polynomials have the well-known property (cid:2)n(ξm) =
δmn where δmn is Kronecker’s delta. In this article, the discretization procedure of (1) by the LPM for any arbitrary γ
parameter is presented keeping in mind that we will eventually take γ as the parameter in (10).
It is also possible to approximate the derivatives of the function y(ξ) by differentiating the interpolant PN(ξ). Further-
more, the derivative values at the nodes ξn may be determined in terms of function values yn= PN(ξn) by means of a
differentiation matrixdefined by
D(k):=(cid:8)dm(kn)(cid:12)= ddξkk(cid:8)(cid:2)n(ξ)(cid:12)(cid:15)(cid:15)ξ=ξm, k=1,2,...,N (16)
for m, n =0, 1, ..., N. The approximate derivative values y(k)=[P(Nk)(ξ0), P(Nk)(ξ1), ..., P(Nk)(ξN)]T may therefore be written
in matrix–vector form
y(k)=D(k)y (17)
where y=[y0, y1, ..., yN]T is the vector of function values at the nodes. In particular, the entries of the first and the
second order differentiation matrices can be obtained as
⎧
dm(1n)= 1⎨⎩ξm−2ξn ψψN(cid:5)N(cid:5)++11((ξξmn)) ifm(cid:7)=n (18)
2 ξ1n(ξn−γ −1) ifm=n
and
⎧
dm(2n)= 1⎨⎩ξm−3ξn [ξ1m(ξm−γ −1)− ξm−2ξn]ψψN(cid:5)N(cid:5)++11((ξξmn)) ifm(cid:7)=n (19)
3 ξ1n[ξ1n(ξn−γ −1)(ξn−γ −2)−N] ifm=n
by making use of (14) and (16) [35].
On the other hand, the three-term recursion [35]
(cid:10) (cid:10)
n(n+γ)ψn−1(ξ)−(2n+γ +1−ξ)ψn(ξ)− (n+1)(n+γ +1)ψn+1(ξ)=0 (20)
for the normalized Laguerre polynomials may be used to determine the zeros of ψn(ξ) and therefore those of Lnγ(ξ). Ac-
tually, running the above recursion over the range n =0, 1, ..., N we obtain an inhomogeneous linear algebraic system
(W −ξI)t=b, or in matrix-vector form
90 H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99
⎡ (cid:10) ⎤⎡ ⎤ ⎡ ⎤
γ +1−ξ − γ +1 0 ··· 0 ψ0 0
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢⎢⎢−(cid:10)γ +1 γ +3−ξ −(cid:10)2(γ +2) ... ... ⎥⎥⎥⎢⎢⎢ ψ1 ⎥⎥⎥ ⎢⎢⎢ 0 ⎥⎥⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ ⎥⎢ . ⎥ ⎢ . ⎥
⎢⎢⎢ 0 −(cid:10)2(γ +2) ... ... 0 ⎥⎥⎥⎢⎢⎢ ... ⎥⎥⎥=⎢⎢⎢ ... ⎥⎥⎥ (21)
⎢ ⎥⎢ . ⎥ ⎢ . ⎥
⎢ ⎥⎢ . ⎥ ⎢ . ⎥
⎢⎢⎢⎣ ... ... ... γ +2N−1−ξ −(cid:10)N(γ +N) ⎥⎥⎥⎦⎢⎢⎢⎣ψN−1⎥⎥⎥⎦ ⎢⎢⎢⎣ 0 ⎥⎥⎥⎦
(cid:10)
0 ··· 0 − N(γ +N) γ +2N+1−ξ ψN bN+1
(cid:10)where ψi = ψi(ξ) and the right-hand side is a vector with only one nonzero component bN+1 =
(N+1)(N+γ +1)ψN+1(ξ). Therefore, if we require ψN+1(ξ) =0 or, equivalently, LγN+1(ξ) =0 then the system reduces
to a standard eigenvalue problem Wt=ξt with the eigenvalue parameter ξ, which provides us the roots ξi, i =0, 1, ..., N
of LγN+1(ξ) as required [2,35].
Since the eigenvector associated to each eigenvalue of W is unique up to a constant factor, the m-th computed eigen-
vector
vm=[v0,m,v1,m,...,vN−1,m,vN,m]T (22)
of the matrix W associated to the eigenvalue ξm is a constant multiple of tm=[ψ0(ξm), ψ1(ξm), ..., ψN−1(ξm), ψN(ξm)]T,
that is vm=atm. The value of a can be determ(cid:10)ined by considering the first entries v0,m and ψ0(ξ(cid:10)m) of the eigenvectors vm
and tm, respectively, since ψ0(ξ) =1/h0=1/ Γ(γ +1) is a constant polynomial we obtain a = Γ(γ +1)v0,m.
Therefore, for n =0, 1, ..., N, the values ψn(ξm) of the orthonormal Laguerre polynomials at the zeros of ψN+1(ξ) may
be computed as
⎡ ⎤ ⎡ ⎤
ψ0(ξm) v0,m
⎢ ⎥ ⎢ ⎥
⎢⎢ ψ1(ξm) ⎥⎥ ⎢⎢ v1,m ⎥⎥
⎢⎢⎢⎢ ... ⎥⎥⎥⎥= (cid:10)Γ(γ +11)v0,m ⎢⎢⎢⎢ ... ⎥⎥⎥⎥ (23)
⎣ψN−1(ξm)⎦ ⎣vN−1,m⎦
ψN(ξm) vN,m
in terms of the computed eigenvector vm of tridiagonal symmetric matrix W.
Then we put the interpolant (13) into (8)
(cid:13)N (cid:8) (cid:12) (cid:13)N
ξ(cid:2)n(cid:5)(cid:5)(ξ)+(γ +1−ξ)(cid:2)n(cid:5)(ξ)+Q(ξ)(cid:2)n(ξ) yn=λ ξα2−1(cid:2)n(ξ)yn (24)
n=0 n=0
and demand its satisfaction at the grid points ξm for m =0, 1, ..., N, to get the discrete representation
Ay=λBy (25)
of the WPL. The general entries Amn and Bmn of the resulting matrices A and B are given by
Amn=ξmdm(2n)+(γ +1−ξm)dm(1n)+Q(ξm;α,c,γ)δmn (26)
and
Bmn=ξmα2−1δmn (27)
respectively. Since B is a diagonal matrix, the generalized matrix–eigenvalue problem in (25) can be replaced with the
standard one
T(cid:25)y=λy (28)
where T(cid:25)=B−1A and its general entry reads as
(cid:8) (cid:12)
T(cid:25)mn=ξm1−α2 ξmdm(2n)+(γ +1−ξm)dm(1n)+Q(ξm;α,c,γ)δmn (29)
where m, n =0, 1, ..., N. Here, the vector yn=[y0,n, y1,n, ..., yN,n]T involves the values of the n-th eigenfunction of (8)
associated with the eigenvalue λn at the nodal points.
H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99 91
By using (18) and (19) the first two terms in (29) can be incorporated to define
⎧
K(cid:25)mn=−16⎪⎨⎪⎩ξ(1nξ12m−ξ−m2α2−ξn{α2)22NψψN(cid:5)N+(cid:5)++11ξ(1(nξξmn[)()γ −ξn)2−1]} iiffmm(cid:7)==nn (30)
which represents the effect of kinetic energy terms independent of a specified potential.
It seems that the evaluation of K(cid:25)mn requires the computation of the derivatives ψN(cid:5)+1(ξn) of the normalized Laguerre
polynomials at the nodes. Fortunately, a nice similarity transformation T =S−1T(cid:25)S in which
S=smδmn=ξm1−α1ψN(cid:5)+1(ξm)δmn (31)
is a diagonal matrix, makes it possible to get rid of such a cumbersome labor. Furthermore, the matrix in (28) reduces to a
symmetric one, say T =S−1(K(cid:25)+Q)S, whose entries are given by
Tmn=Kmn+Qmn (32)
where
⎧
Kmn=−16⎨⎩ξ1n21((−ξξmmα2−ξ{nξ2)n1N)−2α1+ ξ1n[(γ −ξn)2−1]} iiffmm(cid:7)==nn (33)
and
(cid:2) (cid:3)
(cid:9) (cid:11)
Qmn=ξm1−α2Q(ξm;α,c,γ)δmn=ξm1−α2 14ξm− 12(γ +1)− (α1c)2ξmα2−1V ξm1/α/c δmn. (34)
Thus, the eigenvalues of (28), and, hence, the approximate eigenvalues of the WPL can be determined by the symmetric
matrix eigenvalue problem
Tu=λu (35)
since the similar matrices share the same spectrum. Note that, the Laguerre pseudospectral formulation of the WPL leads
to the symmetric matrix eigenvalue problem whose construction requires only the knowledge of the roots of the Laguerre
polynomial LγN+1(ξ).
On the other hand, the n-th eigenvector y of (28) is given by the formula
n
yn=Sun (36)
iimn ptelirems st hoaft tT(cid:25)he[ Snu-t]h = eλig[eSnuv]e. cTtohru su, nth=e [mu0-,tnh, ue1n,ntr, y. ..y,m u,nN,=n]Tyno(fξ mth)eo fs ytmhem ne-ttrhic emigaetnrivxe cTtor= ySn−m1T(cid:25)ayS bsei nwcer itSte−n1 T(cid:25)asS uym=,nλ=u
smum,n, or in nodal notation we have
yn(ξm)=ξm1−α1ψN(cid:5)+1(ξm)um,n (37)
upon using (31). Then (7) determines the original wavefunction
Rn(M,(cid:2))(ξm)=ψN(cid:5)+1(ξm)ξm1+(cid:2)−α1e−ξm/2um,n (38)
at a collocation point ξm. Turning back to the original variable r via (6), the last equation reads as
(cid:9) (cid:11)
Rn(M,(cid:2))(rm)=ψN(cid:5)+1 (crm)α (crm)α+(cid:2)−1e−(crm)α/2um,n. (39)
Thus, we may state the following theorem.
Theorem 2.1. The approximate eigenvalues E(M)of the radial Schrödinger equation in (1)are connected with the eigenvalues of the
n,(cid:2)
linear system in (35)by the formula
En(M,(cid:2))=−(αc)2λn(α,c,γ), n=0,1,... (40)
and the values of the corresponding normalized eigenfunctions Rn(M,(cid:2))(rm)(in L2ρ sense) at the points rm=ξm1/α/care given by
(cid:27)
(cid:10)
Rn(M,(cid:2))(rm)= αcMRn(M,(cid:2))(rm)=− αcM(NΓ+(1γ)(+N1+)γ +1)vvN0,,mm(crm)α+(cid:2)−1e−(crm)α/2um,n (41)
whenever unis the normalized (in Euclidean 2-norm) eigenvector of (35). The v0,mand vN,mare the first and the last entries of (22),
respectively.
92 H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99
Proof. The first part follows from (11). When written in the Sturm–Liouville form, it is not difficult to see that the weight
functions of (1)and (8)are given by ρ(r) =rM−1 and ρ(ξ) =ξγ+α2−1e−ξ, respectively. First using the scaling transformation
in (6)we write
(cid:28) (cid:28) (cid:29)∞(cid:8) (cid:12) (cid:29)∞(cid:8)(cid:10) (cid:12) (cid:29)∞(cid:8) (cid:12)
(cid:28)R(M)(r)(cid:28)2 = R(M)(r) 2rM−1dr= αcMR(M)(r) 2rM−1dr= R(M)(ξ) 2ξMα−1dξ (42)
n,(cid:2) L2ρ(r) n,(cid:2) n,(cid:2) n,(cid:2)
0 0 0
which, upon the use of (7), reads as
(cid:29)∞
(cid:28) (cid:28)
(cid:28)R(M)(r)(cid:28)2 = y2(ξ)ξγ+α2−1e−ξdξ. (43)
n,(cid:2) L2ρ(r) n
0
Applying the N+1 point Laguerre–Gauss quadrature to the function ξα2−1y2(ξ), we obtain, in the limiting case when
n
N→∞
(cid:29)∞
(cid:28) (cid:28) (cid:13)N
(cid:28)Rn(M,(cid:2))(r)(cid:28)2L2ρ(r)=0 ξα2−1yn2(ξ)ξγe−ξdξ=Nl→im∞m=0ξmα2−1yn2(ξm)ωm (44)
where
ωm= (N+1)(N1+γ +1)ψ2ξ(mξ ), m=0,1,...,N (45)
N m
are known as the Christoffel numbers or weights of the quadrature in terms of the normalized Laguerre polynomials [13].
Now, the differential-difference relation [35]
(cid:10)
ξψn(cid:5)(ξ)=nψn(ξ)− n(n+γ)ψn−1(ξ) (46)
of the normalized Laguerre polynomials with n =N+1 and ξ=ξm leads to
(cid:10)
ψN(cid:5)+1(ξm)=−ξ1 (N+1)(N+γ +1)ψN(ξm) (47)
m
since ψN+1(ξm) =0 for m =0, 1, ..., N and hence, (37) rewritten as
(cid:10)
−1
yn(ξm)=− (N+1)(N+γ +1)ψN(ξm)ξmαum,n. (48)
Inserting (45) and (48) into the right hand side of (44) we obtain
(cid:28) (cid:28) (cid:13)N
(cid:28)Rn(M,(cid:2))(r)(cid:28)2L2ρ(r)=Nl→im∞m=0um2,n=Nl→im∞(cid:9)un(cid:9)22=Nl→im∞1=1. (49)
On the other hand, (39) together with (47) and (23) gives (41) which completes the proof. It is clear from (43) and (49)
that
(cid:28) (cid:28)
(cid:28)yn(ξ)(cid:28)L2ρ(ξ)=1 (50)
is also true under the same assumption where ρ(ξ) =ξγ+α2−1e−ξ. (cid:2)
3. Numerical examples
It is clear from (8) and (10) that the spectrum of the eigenvalue problem remains unchanged for a prescribed value of
the sum 2(cid:2) +M. Hence the eigenvalues in M dimensions denoted by E(M) are degenerate in such a way that
n,(cid:2)
E(2) ≡E(4)
n,1 n,0
E(2) ≡E(4) ≡E(6)
n,2 n,1 n,0
.
.
.
En(2,(cid:2))≡En(4,(cid:2))−1≡En(6,(cid:2))−2≡...≡En(2,2(cid:2)−2)≡En(2,1(cid:2))≡En(2,0(cid:2)+2) (51)
H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99 93
Table1
The energy eigenvalues En(3,(cid:2)) of the potential, V(r) =r2+v4r4, as a function of v4where αopt=2.
v4 copt N n l En(3,(cid:2))
10−4 1 8 0 0 3.0003748969361210983378468299
1 40 25 1 105.4103438524395595536210145910
1 66 50 5 214.6749649908048220255705112162
1 124 100 10 429.514482011916008399592238938
1 3 30 0 10 54.1849846104544399241234801756a
3 65 25 5 483.0222074133947094286082707291
2.5 100 50 1 1062.889853853655671834975735691
2 160 100 0 2604.43248571463930745940568155
104 9 30 0 0 81.903316953284467567471308555
9 80 25 1 9253.92349941549971482158637398
12 100 50 5 23756.5339836909761084585149553
13 162 100 10 59302.0603134555154912941546049
a Reference [36]: E(3) =54.184 984 610 454 439 924 123 480 175 7.
0,10
Table2
The effect of second optimization parameter αon the accuracy of ground state energies of isotropic quartic oscillator when v4=10−4 and
the Airy equation.
α (EN0(3,0)=o8f, icsoopttr=op1ic) quartic oscillator (EN0(1,1)==66E,0 c,0opotf= A6ir)y equation
3 3.066750 9.111969
2.5 3.026398 4.842756
2 3.0003748969361210983378468299 2.632534
1.5 2.993411 2.338107893
1 2.946460 2.3381074104597670384891972524
0.5 3.027127 2.338111
for even values of space dimension M where E(2) is single in the system. Similarly, if M is odd, then
n,0
En(3,0) ≡En(1,1) ≡En,0
E(3) ≡E(5)
n,1 n,0
E(3) ≡E(5) ≡E(7)
n,2 n,1 n,0
.
.
.
En(3,(cid:2))≡En(5,(cid:2))−1≡En(7,(cid:2))−2≡...≡En(2,2(cid:2)−1)≡En(2,1(cid:2)+1)≡En(2,0(cid:2)+3) (52)
where En,0 are the eigenvalues of the Schrödinger equation in (2). The degenerate structure of the spectrum of (1) now
suggests evidently that we may consider only two- and three-dimensional cases, without any loss of generality.
As a first example we consider the M-dimensional isotropic quartic oscillator
V(r)=r2+v4r4 (53)
where v4>0. Tas¸eli and Zafer [36] expanded the wave function into a Fourier–Bessel series to solve the radial Schrödinger
equation with isotropic polynomial potentials and Tas¸eli [32]proposed an alternative series solution to the isotropic quartic
oscillator in M dimensions. It seems that the transformation (6) with α=2 is the most suitable one since the potential
contains even powers of r only. In Table 1, we present eigenvalues of isotropic quartic oscillator in 3-dimensions for some
pairs of (n, l). In all tables, n and (cid:2) stand for the quantum numbers of the state, N the truncation order for which the
desired (machine) accuracy of the corresponding eigenvalue is obtained, and c and α denote some scaling or optimization
parameters which may be used to accelerate the convergence rate of the method. The accuracy of the results in all tables
reported here has been checked by inspecting the number of stable digits between two consecutive truncation orders. For
comparison, at the bottom right of Table 1we included a result from [36] when v4=1.
On the other hand, the effect of the second optimization parameter α is reported in Table 2. In order to regain the
machine accuracy for E(3) when, for example, α=1 we should choose N=50 and c=15 at the same time. Therefore,
0,0
we may say that our prediction of α is experimentally justified by Table 2. If further c=1, that is, in the absence of
scaling factors, the truncation order N increases dramatically. For unbounded domains scaling factors are crucial to obtain
the desired accuracy with the possible smallest truncation order.
The second example is an odd polynomial potential
94 H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99
Table3
Several eigenvalues of the Airy equation when αopt=1. The last column includes the negatives of first ten zeros of Airy function Ai(x).
copt N n En,0=En(1,1) Reference[1]
6 66 0 2.3381074104597670384891972524 2.33810741
1 4.0879494441309706166369887014 4.08794944
2 5.5205598280955510591298555129 5.52055983
3 6.7867080900717589987802463845 6.78670809
4 7.9441335871208531231382805558 7.94413359
5 9.0226508533409803801581908399 9.02265085
6 10.0401743415580859305945567373 10.04017434
7 11.0085243037332628932354396496 11.00852430
8 11.9360155632362625170063649029 11.93601556
9 12.8287767528657572004067294072 12.82877675
6 66 10 13.6914890352107179282956967794
6 70 20 21.2248299436420953684599203593
8 110 30 27.5883878098824448119503644141
9 135 40 33.2848846819014018796197398960
10 155 50 38.5288083050942488226298967447
15 280 100 60.8589317646089237955214557538
V(r)=r, r∈(0,∞) (54)
for which the Schrödinger equation reduces to the Airy equation −y(cid:5)(cid:5)+ry =λy when we set M=1 and (cid:2) =0 or (cid:2) =1
in (1). In this case, eigenvalues are given implicitly by Ai(−λ) =0 where Ai(x) is the Airy function [8]. Several states of the
Airy equation in one dimension is reported in Table 3. The last column includes the negatives of zeros of the Airy function
taken from [1]. It is clear from Table2that for the Airy equation the optimum value of the second scaling factor is αopt=1.
As a third example, we consider the exponential cosine partially screened Coulomb potential (ECPSC)
(cid:2) (cid:3)
1
V(r)=−2ZVec(r,ν,μ)−2Zas −Vec(r,ν,μ) , Z>0, Zas>0 (55)
r
where
Vec(r,ν,μ)= 1e−νrcos(μr) (56)
r
with the two screening parameters ν and μ[17]. In particular, when Zas=0the potential reduces to the exponential cosine
screened Coulomb potential (ECSC) V(r) =−2Zexp(−νr) cos(μr)/r. If further, μ =0 at the same time, it is known as the
Yukawa potential V(r) =−2Zexp(−νr)/r. On the other hand, Zas=Z corresponds to the pure attractive Coulomb potential
V(r) =−2Z/r which has countably many discrete states given by
E(M)= −4Z2 , n=0,1,... (57)
n,(cid:2) (2n+2(cid:2)+M−1)2
together with the continuous spectrum over the entire positive real axis.
These potentials have been subject of several studies. For example, Lai [19] determined several states of ECSC within the
framework of the hypervirial Padé scheme. Tas¸eli [33] used modified Laguerre basis for the ECSC and Yukawa potentials.
Ixaru and co-workers [17] developed accurate, robust and safe approach for ECPSC.
Several eigenvalues of the ECPSC potential in three dimensions are reported in Table 4 for the parameter values Z=50,
Zas=1 and ν=μ =0.025 when (cid:2) =0, 10. It is clear from Table 4 that if we keep α=1, higher levels become expensive
or even impossible to obtain. In this stage, transformations (6) and (7) together with the decreasing character of c suggests
that we should take smaller α values. Indeed, when we take α=0.7 while keeping c=1, it is possible to obtain higher
states with relatively small truncation orders.
Table5demonstrates the bound states of the ECSC potential. In all computations we used quadruple-precision arithmetic
on a main frame computer with machine accuracy of 32 digits, by truncating the results to 27–28 significant figures. For
comparison, at the bottom right of the table, we attached results from [33] and [21] which are in good agreement with
those of the present study to the accuracy quoted.
Fig. 1 illustrates the two normalized eigenfunctions Ψ30,0(r) and Ψ4,0(r) of the ECPSC and ECSC potentials, respectively,
for the specified parameter values. Our results are in good agreement with those of MATSLISE, which is typical for all
eigenfunctions considered in this study.
The present method works well also for the pure attractive Coulomb potential. For instance, in three dimensions when
Zas=Z =1, (cid:2) =0, α=1 and N=500, the eigenvalues En(3,0) with 2 ≤n ≤125 are obtained within the machine accuracy
by taking c=0.05 while c=0.005 leads to the same accuracy for E(3) with 130 ≤n ≤430. That is, the optimization
n,0
parameters may be seen as the position identifier of a flash lamp which illuminates the specific part of the eigenvalue
H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99 95
Table4
Several states of the ECPSC potential in three dimensions when Z=50, Zas=1and ν=μ =0.025as (cid:2)varies.
αopt copt N n (cid:2) En(3,(cid:2))=En(1,(cid:2))+1=En,(cid:2) En,(cid:2)(Reference[17])
1 90 20 0 0 −2497.5500006121173026119994770 −2497.5500006120
60 20 1 −622.55000855817107243365113285 −622.550008557
30 20 2 −275.327819864885534764663309750 −275.327819864
8 30 10 −18.2182548645294488912560639439 −18.218254864
5 40 20 −3.3012939237449879468250075666 Notreported
3 62 30 −0.4779793951088033625235770389 Notreported
0.5 325 50 −0.001.5318333743193636649758064 −0.001531833374
0.7 1 350 100 −0.00017450705849151 −0.00017450705849
1 8 25 0 10 −18.2144512404084595291668336118 −18.21445124040
8 25 1 −14.9165994843486359598726352512 −14.9165994843
8 25 2 −12.3512992295081145310345160781 −12.3512992294
5 40 10 −3.2899432840178999972361084807 −3.28994328401
4 80 20 −0.4601174206377746472193161698 Notreported
0.7 200 30 −0.0033680864235131844680688099 Notreported
.3 355 50 −0.000699631092664645508 −0.000699631
0.7 1 350 100 −0.000129394044816757730 −0.0001293940
Table5
Bound energy eigenvalues of the ECSC potential in three dimensions when Z=1and ν=0.05, as μvaries. The case, μ =0corresponds to the Yukawa
potential.
μ copt αopt N n (cid:2) En(3,(cid:2))
0 0.6 1 65 0 0 −0.9036328570490110877124341515a
1 −0.1635423915905062483469788275
2 −0.0387051096295046845907959936
3 −0.0061833198003226429693179006
1 0.5 250 4 −0.000003138989336707735885268b
0.4 1 56 0 1 −0.1614807740755692194242054872
1 −0.0371155037668119932097879878
2 −0.0051961177051437079305223825
0.4 1 51 0 2 −0.0338311411396316857722295165
1 −0.0031617432537420099057678961
0.05 1 1 55 0 0 −0.9002349328413753360905526418
1 −0.1528991925004954887681337617
2 −0.0231511284141215916951113223
55 0 1 −0.1521180248834620947098731497
1 −0.0218586596451123224971294624
55 0 2 −0.0191097586454757388818428140
a Appropriately scaled result from Reference [33]: E(3)=−0.903 632 857 049 011 087 712 434.
0,0
b MATSLISE [21]: E4(3,0)=E4(1,1)=E4,0=−0.000 003 138 989 337.
Fig.1.Normalized eigenfunction Ψ30,0(r) =R(310),1(r) of the ECPSC potential corresponding to E3(10),1=E30,0 when Z=50, Zas=1, ν=μ =0.025 (left).
Normalized eigenfunction Ψ4,0(r) =R(41,1)(r)of the ECSC potential corresponding to E4(1,1)=E4,0when Z=1, ν=0.05, μ =0(right).
96 H.Alıcı,H.Tas¸eli/AppliedNumericalMathematics87(2015)87–99
Table6
Bound states of the Hulthén screening potential in three dimensions when Z=50and ν=0.025, as (cid:2)varies (αopt=1).
copt N n (cid:2) En(3,(cid:2))=En(1,(cid:2))+1 En,(cid:2)(Reference[17])
60 60 0 0 −2498.7501562500000000000000000198
1 −623.7506250000000000000000000007
2 −276.5291840277777777777777777795
3 −155.0025000000000000000000000016
0.5 350 59 −0.0069444444444444444444444444
60 −0.0032686525799516259070142435
61 −0.0009892039542143600416233091
62 −0.0000378322625346
8 30 0 10 −19.4243353045265954430260492336 −19.42433530452
1 −16.1278839619036748840561546591 −16.1278839619
2 −13.5635796169491111377056995371 −13.5635796169
3 −11.5300024482500496959324878452 −11.5300024482
1 250 47 −0.0156756090848458797234145628 Notreported
48 −0.0092609697111221593780144693 −0.00926096971
49 −0.044487838376324670737199982 −0.0044878383
50 −0.0011881662588910964623754334 −0.001188166259
column. Therefore, adjusting the position of the flashlight, it is possible to lightening the whole set of eigenvalues for a
prescribed truncation order N.
The next potential we consider is the partially screening Hulthén potential
(cid:2) (cid:3)
1
V(r)=−2ZVH(r,ν)−2Zas −VH(r,ν) , Z>0, Zas>0 (58)
r
where
νe−νr
VH(r,ν)= 1−e−νr (59)
in which ν is the screening parameter. It reduces to the Hulthén screening potential [14] when Zas=0 which is exactly
solvable when M=(cid:2) =1 (or (cid:2)¯=0). In this case, bound states are given by
(cid:2) (cid:3)
En(1,1) =En,0=− n+Z 1 − (n+21)ν 2, n=0,1,...,k (60)
√
where k =(cid:11) 2Z/ν(cid:13) −1 [8]. Here (cid:11)a(cid:13) denotes the integer part of a real number a.
The partially screening Hulthén potential is considered by Ixaru, De Meyer and Vanden Berghe [17]. Hulthén potential is
studied by many authors. For instance, Roy [26]applied the generalized pseudospectral approach to approximate the bound
states, Stubbins [30]used the generalized variational method to compute the eigenvalues for n ≤6, Bayrak and Boztosun [5]
used asymptotic iteration method for any (cid:2)state and Gönül and co-workers [10]considered the potential in the Hamiltonian
hierarchy picture to approximate the eigenvalues when (cid:2) (cid:7)=0.
Table 6 illustrates the bound states of the Hulthén screening potential in three dimensions when Z=50 and ν=0.025.
For (cid:2) =0 the results are in good agreement with the exact eigenvalues. The last column includes results from [17] when
(cid:2) =10, which are also in good agreement to the accuracy quoted.
Finally, we take into account the Woods–Saxon potential defined by
(cid:2) (cid:3)
50 5t
V(r)=− 1− (61)
1+t 3(1+t)
where t=exp[3/5(r−7)]. The potential has been considered by several authors; Zakrzewski [40] used a power series
method, Lo and Shizgal [23]applied quadrature discretization method, Shao and Wang [27] considered Obrechkoff one-step
method to approximate the eigenvalues of the problem.
Bound states of the Woods–Saxon potential in two dimensions when (cid:2) =2 and in three dimension when (cid:2) =0 are
reported in Table 7. In both cases, there exist 14 discrete states before the start of continuous spectrum over the entire
positive real axis.
Fig.2shows the first eigenfunction Ψ1,0(r) =R(11,1)(r)of the Woods–Saxon potential by using N=120points for different
values of the optimization parameter c while keeping α=1. At this truncation order, the optimum c value is copt=45.
Notice that, copt collects the grid points to the region where the eigenfunction is nonzero. Neither the points are wasted in
the region where the wavefunction is too close to zero, nor they are insufficient to recover the shape of the eigenfunction.
In this way, it reduces the number N of collocation points used to get the desired accuracy. Here, for c=15 and c=45,
the energy E1,0=E1(1,1) is correct to 20 and 27 digits respectively, but no convergence occurs when c=100 for the same
truncation order of N=120.