Very-high-precision normalized eigenfunctions for a

class of Schrodinger type equations

arXiv:1105.1460v1 [math-ph] 7 May 2011

Amna Noreen , Kare Olaussen

Institutt for fysikk, NTNU, Trondheim Norway

AbstractWe demonstrate that it is possible to compute wave

function normalization constants for a class of Schrodinger type
equations by an algorithm which scales linearly (in the number of
eigenfunction evaluations) with the desired precision P in decimals.
KeywordsEigenvalue problems; Bound states; Trapezoidal rule;
Poisson resummation

I. I NTRODUCTION
N a recent paper[1] it was demonstrated that it is possible to
solve some eigenvalue problems of the Schrodinger type to
almost arbitrary high precision. The cases presented explicitly
were (i) the ground state eigenenergy of the anharmonic
oscillator,
00 + x4 = ,
(1)

which was found to an accuracy of more than 1 000 000

decimal digits, (ii) the eigenstate number 50 000 of the same
equation where the eigenenergy was found to an accuracy of
more than 50 000 decimal digits, and (iii) the lowest even and
odd parity states of the double-well potential,
00
s2
+ (x2 1)2 = ,

(2)

where the two eigenenergies were found to an accuracy of

more than 30 000 decimal digits for the case s = 1/50 000.
In the latter case the two eigenenergies are degenerate to
almost 29 000 decimals, with the difference directly computable by the WKB method. The 10th order WKB expansion
given in [2] provide an accuracy of about 48 decimals for the
difference, all of which agrees with the difference between the
two numerical calculated eigenenergies.
Also the eigenenergies of the highly excited states of
equation (1) can be calculated by the WKB method. The 12th
order WKB expansion given in [3] provide an accuracy of
about 67 decimals, all of which agrees with the numerical
result.
It is certainly difficult to find physical systems where one
need to know eigenenergies to tens of thousands of decimals
or more. However, if one need to compute the wavefunction
very accurately to very high values of x (to f.i. evaluate
matrix elements) the value depends extremely sensitively on
the eigenvalue parameter. Thus, in our opinion, algorithms for
very-high-precision evaluation of eigenvalues may be useful in
combination with routines for very-high-precision evaluation
of matrix elements and normalization integrals.
In this paper we demonstrate that the latter can be achieved
very simply, with the number of eigenfunction evaluations
growing only linearly with the desired precision P in decimal digits. Each eigenfunction evaluation will in turn require a number of high-precision multiplications which grows
linearly with P , and each such multiplication require a
Amna.Noreen@ntnu.no
Kare.Olaussen@ntnu.no

CPU time which scales asymptotically between P 1.6 and

P log P log log P (depending on which high-precision multiplication algorithm is used). Thus, the total time to evaluate
one normalization integral to P decimals precision can be
expected to increase somewhat faster than P 3 . However, since
the eigenfunction evaluations are independent they can be run
in parallel.
The rest of this paper is organized as follows. In section II
we make some general remarks on numerical integration
rules. We base our analysis on the Euler-Maclaurin summation
formula and the Poisson resummation formula. Our conclusion
is that the trapezoidal rule is not only the simplest one but also
the best one. For integrals over finite intervals some there are
endpoint corrections which should be considered separately;
these corrections vanish for our normalization integrals. In
section III we consider some simple example cases similar
to our wavefunction normalization integrals. For an infinite
(in principle) integration range and a fixed number M of
integration steps one must strive for a balance between the
error (h) due to using a finite stepsize h, and the error (xmax )
due to covering only a finite integration range. These errors can
be estimated by respectively analysing the Fourier transform
f(p) of the integrand f (x) as |p| (by the method of
steepest descent), and the asymptotic behaviour of f (x) as
|x| . It appears that both analyses can be extended
to wavefunction normalization integrals by use of the WKB
approximation. This extension is done in section IV.
II. R EMARKS ON INTEGRATION FORMULAE
There is no scarcity of numerical integration formulae in
the
R b literature[4]. Standard choices are rules for I {f }
dx f (x) which reproduces integrals of polynomials below
a
a certain order exactly,
1
(f0 + f1 ) h T (h),
2
1
I {f } (f0 + 4f1 + f2 ) h S(h),
3
(3)
3
I {f } (f0 + 3f1 + 3f2 + f3 ) h S3/8 (h),
8
2
I {f }
(7f0 + 32f1 + 12f2 + 32f3 + 7f4 ) h B(h),
45
I {f }

with h = (b a)/M and fm = f (a + mh) when there

is M + 1 terms in the integration rule. These are known
respectively as the trapezoidal, Simpsons, Simpsons 83 , and
Booles rule. They are automatically exact for polynomials
which are antisymmetric about the midpoint x
= (b a)/2.
The M independent coefficients are chosen to give exact
results for all symmetric (about x
) polynomials of order below
2M . However, as M increases the weight coefficients develop
in a suspicious way. To integrate functions over large intervals
to very high precision it seems dubious to extend the procedure
above.

rules

A. Extended Simpsons rule

An alternative way to handle integration over large intervals
is to divide them into many smaller ones, and apply one of the
rules above to each of the subintervals. We denote this by a
bar over the rule, T (h) T(h) etc. A rather common choice
is to extend Simpsons rule, leading to the formula
Z

dx f (x) S(h)
(f0 + 4f1 + 2f2 +
3
(4)
+2fM 2 + 4fM 1 + fM ) h,

which looks curious to any member of an equal society. What

is wrong with half of the points? Are they of the wrong
gender? Which mysterious force of numerical error analysis
leads to this spontaneous breakdown of translation invariance?

#
Z b
M
1

X
X
1
1
B2k (2k1)
f0 +
fm + fM h =

+
dx f (x) +
2
2
(2k)!
a
m=1
k=1
Z b
(1)
(3)
(5)
(7)
=
dx f (x) +

+ ,
(7)
12
720
30240
1209600
a

"



where (k) = f (k) (b) f (k) (a) hk+1 , and the ellipsis in
the first line denote non-perturbative terms depending on

how f (x) varies in the full integration range. With S(h)

=
4
1
T
(h)

T
(2h)
one
finds
the
corresponding
quantity
for
the
3
3
extended Simpsons rule:
Z b

B2k (2k1)
1X
4 22k
S {f } =
dxf (x) +

+ .
3
(2k)!
a
k=1
(8)

105
1010
1015
1020
1025

(h) I4

I.e., the (1) -coefficient has been eliminated at the cost of

increasing all the remaining coefficients1 .
A natural way to eliminate the (1) -coefficient would be
to calculate f 0 explicitly at the endpoints, or use a numerical
approximations like

1015

(h) I12

1
f (x + ) f (x )
f (3) (x) 3 + . . . ,
2
6
3f (x) 4f (x + ) + f (x + 2)
1
0
f (x) =
+ f (3) (x) 3 + . . . ,
2
3
3f (x) 4f (x ) + f (x 2)
1
f 0 (x) =
+ f (3) (x) 3 + . . . ,
2
3
f 0 (x) =

1030
1045
1060
1030

(h) J

1060
1090
10120
104

103

Simpsons rule
Simpsons (fit)
Trapezoidal rule
Trapezoidal (fit)
102

101

Fig. 1 Comparison of accuracy of the extended trapezoidal and

Simpsons rules for some integrands with very smooth boundary
behaviour.

In our opinion only T(h) is a logically sensible procedure.

We claim that this is also the best procedure. To support and
exemplify this claim we have evaluated the integrals below:
Z

n
dx 1 x2 ,
0

Z 1 
=
dx 1 tanh

In =
J

(5)
2x 1
1 (2x 1)2


,

(6)

by the T(h) and S(h)

= 43 T(h) 13 T(2h) rules. The relative
errors of the two methods as function of discretization length h
are shown in figure 1. The result is as expected that S behave
as well as T with twice the discretization length h.

for sufficiently small . One should take somewhat smaller

than h. The modified trapezoidal rule described in [5] is
generally worse than the extended Simpsons rule because of
an inaccurate approximation of f 0 . However, there are few
reasons to use the same discretization length for computing
derivatives at the endpoints as for computing the bulk contribution to the integral. To avoid introduction of significant new
errors it is sufficient
to choose < h/2 in the first formula.
By choosing = h/ 20 in the last two one also eliminates
the (3) -term in equation (7).
C. Poisson resummation formula
There exist expressions for the non-perturbative terms in
equation (7), but we find the Poisson resummation formula
more clarifying: Assume f (x) is an entire, integrable function,
with Fourier transform
Z
f(p) =
dx f (x) eipx .
(9)

Then one version of the Poisson resummation formula reads

X
X
2k
).
(10)
h f (mh) =
f(
h
m=
k=

This can be used to estimate the accuracy of the numerical

integration formula
Z

X
X 2k
dx f (x) =
h f (mh)
f(
).
(11)
h

m=
k6=0

B. Euler-Maclaurin summation formula

The examples in figure 1 are untypical in the sense that the
integrands are very smooth at the endpoints (which is typical
for the normalization integrals we are primarily interested in).
One may use the Euler-Maclaurin summation formula to see
which quantity is actually computed by the various integration

When f (x) is an entire function its Fourier transform f(p)

will vanish faster that any inverse power of p as p . This
means that the numerical approximation to the (infinite range)
1 Note that the error terms given for the extended rules in reference [4]
are in disagreement with our results.

integral of an entire function will converge very fast towards

the exact value as h 0.
Integrals of f (x) over a finite range [a, b] can be viewed
as integrals of g(x) (b x) (x a) f (x) over an
infinite range. But g(x) will usually be discontinuous, with
a fourier transform g(p) which vanishes only algebraically as
p . The perturbative endpoint corrections of the EulerMaclaurin formula can be used to account for these algebraic
terms in a systematic way. But there is no point in making
endpoint corrections beyond the error in the bulk contribution,
2
the latter being of magnitude f( 2
h ) + f ( h ).
III. E XAMPLE INTEGRALS
In this section we will analyze the behaviour of some simple
cases which are similar to typical ground state normalization
integrals.
2

A. ex
Consider first the ground state of the harmonic oscillator,
2
f (x) = ex , in which case

2
f(p) = ep /4 .
(12)
Hence we have
Z

2
2
2
2
dx ex =
he(mh) 4 e /h . (13)

m=

In practise we must approximate the infinite sum by a finite

one,
Z
M
X
2
2
dx ex = h + 2h
e(mh)

m=1
2

+ 2he(M +1)

h2

4 e

/h2

If we make h too small the first correction term (on the second
line above) becomes too large. If we make h too large the second correction term becomes too large. I.e., the optimal choice
2 2
2
2
of h occurs approximately when e(M +1) h = e /h ,
r

.
(14)
h=
M +1
This leads to an expected error of order
(M ) = e(M +1) .

(15)

I.e., to evaluate this integral numerically to P decimals accuracy one must choose
M +1

log 10
P 0.73 P.

(16)

x2n

B. e
2n
In more generality consider f (x) = ex with n a positive
integer. In this case the integral is
 
Z
1
1
x2n
dx e
=
Kn
.
(17)
n
2n

We use the saddle point method to estimate its Fourier

transform
Z
Z
x2n +ipx

f (p) =
dx e

dx e(x,p) .
(18)

The saddle point equation becomes

0 (xs , p) = 2nx2n1
+ ip = 0,
s

with 2n 1 solutions

xs = ei(4k+1)/(4n2)

(19)

 p 1/(2n1)
,
2n

k = 0, 1, . . . 2n 2.
(20)
The two most relevant saddle points2 occur for k = 0 and
n 1, at which the real part of the exponent (x, p) is (when
p = 2/h)

 

2n/(2n1)
Re (xs , p) = (2n 1) sin
4n 2
nh
an h2n/(2n1) .

(21)

I.e., the leading error due to a finite stepsize h is of magnitude

2n/(2n1)
eRe (xs ,p) = ean h
(apart from a prefactor of less
importance). On the other hand, including only the M first
2n
terms of the sum leads to an error of magnitude e[(M +1)h] .
For a given M we should choose h to balance these errors,
i.e.
(11/2n)
h = bn (M + 1)
,
(22)
2
h

i
(2n1)/4n

1/2n

with bn = n
(2n 1) sin 4n2
. This
leads to an error of magnitude
(M ) = ecn (M +1) ,
(23)
h

i11/2n

with cn = n (2n 1) sin 4n2

. I.e., to evaluate
the integral numerically to P decimals accuracy one must
choose
log 10
M +1
P (0.12 + 0.467 n) P.
(24)
cn
where the numerical approximation is very good for n 2.
2

C. e(x

a2 )2

2
2
Finally consider f (x) = e(x a ) . In this case
Z
2
2 2
I(a) =
dx e(x a )

i
h
1 4
= a e 2 a 12 K 14 ( 12 a4 ) + I 41 ( 12 a4 )

21 (1/4) as a 0,

/a
as a .

(25)

The saddle point approximation to the Fourier transform

Z
Z
(x2 a2 )2 +ipx

f (p) =
dx e

dx e(x,p) ,
(26)

leads to the saddle point equation

0 (xs , p) = x3s + a2 xs + i
We introduce
xs = y +

p
= 0.
4

a2
,
3y

(27)

(28)

which leads to a quadratic equation for y 3 ,

p
a6
y3 i +
.
4 27y 3

(29)

2 The deformation of the integration path so that it passes through these

saddle points is an interesting exercise[6]


3
It is convenient to rewrite p so that (p/8) = a/ 3 sinh 3.
Then the solutions for equation (29) becomes
! 
r
3
p2
a6
a
p
3

+
= e(+i/6) . (30)
y =i
8
64 27
3

a2
a
= e(+i/6) ,
3y
3

which in both cases of leads to the solution

2a
xs = cosh ( + i/6) .
3

1/2
,
xmax = a 1 + s

0
1/2
xmin =
a (1 s)

(36)
if s 1,
otherwise.

(37)

I.e., we must use M = (xmax xmin )/h evaluation steps in

the numerical integration.

I.e.,
a
y = e(+i/6) ,
3

(4/3) a4 sinh2 cosh 2,

(31)

(M )

IV. WAVEFUNCTION NORMALIZATION INTEGRALS

Our investigation of the example integrals gives us confidence that we can obtain a fairly good a priori estimate
of the obtainable precision (M ) at a given number M of
discretization steps, at least asymptotically for large M . For
this we need to (i) estimate the behaviour of the Fourier
f2 (p), of the integrand at large p = 2/h (to find
transform,
the obtainable accuracy (h) at a given stepsize h), and (ii)
estimate how the integrand decays away from its maxima (to
find the required x-range of summation for the same accuracy).
Both quantities can be obtained to reasonable accuracy by
use of the WKB approximation.

10100

10200

A. WKB estimates of wavefunctions for x2n potentials

For large x an estimate of solutions to the eigenvalue
problems


00 + x2n E = 0,
(38)

(x2 a2 )2

Obtained (M )

can be written in the form3

 R


x
(x) = exp x0 dt t2n E



1
x x2n E ,
C exp n+1

x4
Obtained (M )

x2
Obtained (M )

10300
0
50
100
150 M
200
Fig. 2 Predicted (lines) and obtained (points) precision as function of
M , for the threeR example integralsR discussed in
the text. In order from
2
4 R
2
2 2
top to bottom: dx e(x a ) , dx ex , dx ex . The number
of function evaluations in each case is M + 1. As can be seen, the
obtained precision agrees well with the theoretical estimate.

where x2n
0 = E, and

Z
n
E
dt
n+1 x
t2n E
 1 n10

B( 2 , 2n ) (n+1)/2n
= exp
E
2(n + 1)


 

1
tan
,
= exp
N+
2
2n
2


C = exp

There are two more ways to take the cube root. One of them
leads to an equally relevant saddle point,
2a
xs = cosh ( + i5/6) ,
3

(32)

while the last one is irrelevant. The real part of the exponent
at the relevant saddle points is
4
Re (xs , ) = a4 sinh2 cosh 2.
3

(33)

This provides an error estimate in parametric form: If we

choose a finite stepsize

27
h=
,
(34)
3
4 a sinh 3
the corresponding error will be of magnitude
4

(h) e 3 a

sinh2 cosh 2

(35)

The error caused by summing over only a finite range

of x-values, 0 xmin x xmax , should be chosen to be of the same magnitude as (h). I.e., with s

(39)

with B(x, y) = (x) (y)/(x + y) the Beta function. We

have evaluated the WKB integral using partial integration, and
let x in the remainder. Further, in the last equality
we have assumed E EN to be eigenvalue number N , and
evaluated it by the WKB approximation.
We use this approximation to estimate the Fourier transform,
Z
f
2
(p) dx (x)2 eipx ,
(40)
by the saddle point method. The saddle point equation can be
written
2
x2n
(41)
s = E (p/2) ,
giving




 f2 
n
Re (ipxs )
 (p)  C 2 exp n+1
n


1/2n o
n

= C 2 exp n+1
sin 2n
p (p/2)2 E
.
(42)
3 We

ignore the algebraic prefactor (x2n E)1/4 .

Thus, to compute the normalization integral to P decimals

precision one must choose p = 2/h so that the right hand
side of (42) becomes equal to 10P (or smaller). The value
of h is best found numerically. But note that one at least must
have (p/2)2 > E, or
h1 > 1 E 1/2 .
For large x the normalization integrand behaves like

 Z x p
2
2n
dt t E
(x) exp 2
 0

p
2
C 2 exp n+1
x x2n E .
(43)
Thus, to compute the normalization integral to P decimals
precision one should integrate to a value x = xmax for which
the right hand side of (43) becomes less than 10P . The value
of xmax is best found numerically. But note that one at least
must have x2n
max > E, or
xmax > E 1/2n .
This provides a lower limit on the number of integration steps,
M +1=

N + 12
xmax
1
> E (n+1)/2n R 1
. (44)

2n
h

u
du
1

Here we have assumed E EN to be eigenvalue number N ,

and used the WKB approximation to estimate its value. Since
N is the number of nodes in the eigenfunction, the inequality
(44) quite reasonably says that the number of evaluation points
must be larger than the number of oscillations.
B. Ground state wavefunction of x4 potential
Specializing the results of the previous section to the ground
state of the x4 -potential, and approximating E0 by zero, one
finds from equation (42) that


3/2
 f2 
(45)
 (p) ep /3 ,

(x)2 e2x

/3

(46)

From this one deduces that the optimal choices for h = 2/p
and xmax = M h are so that
e

23 (M +1)3 h3

I.e.,
h = 21/9 1/3 (M + 1)2/3 1.58 (M + 1)2/3 .

(47)

The corresponding estimated error becomes

13 24/3 (M +1)

(M ) = e

2.64 (M +1)

(48)

I.e., if one wants to compute the normalization integral to

P decimals precision one has to choose
M +1=

105

E0
E100
EWW

104

103

102
102

Pest

3 ln 10
P 0.87 P.
24/3

(49)

103

Fig. 3 This figure displays the computer time (measured in milliseconds) required to compute some normalization integrals to P
decimals estimated precision. The cases considered are (i) the ground
state wavefunction of the x4 -potential (E0 ), (ii) the 100th excited state
wavefunction of the x4 -potential (E100 ), and (iii) the ground state
wavefunction of the (x2 1)2 -potential (EWW ) with s = 1/100.

C. Excited state wavefunctions of x4 potential

The main effect of having highly excited states EN lies
in the extra factors C 2 in the error estimates, with C found
in equation (39). Clearly C becomes large for large N . To
account for this factor we must replace 10P by C 2 10P
in the error analysis. I.e., make the replacement

 
1
P log 10 P log 10 + tan
(50)
N+
2n
2
in the expressions like (49). There are also E-dependent
corrections to the remainding factors, but they are best treated
numerically.
D. Ground state wavefunction of (x2 1)2 potential
We next consider the lowest even parity eigenstate of
s2 00 + (x2 1)2 = ,

and from equation (43) that

13 (2/h)3/2

t [ms]

(51)

with s is small and positive. By WKB analysis one finds that

the solution behaves like
(x) e(x1)

(x+2)/3s

(52)

(up to an algebraic prefactor) in the region of interest (Rex >

0, |(x 1)2 (x + 2)/3s|  1). We again estimate Fourier
transform of (x)2 by the saddle point method, assuming the
relevant saddle point to be in a region where the approximation
(52) is valid. The saddle point equation,


d 2
0 (x, p) =
(x 1)2 (x + 2) ipx = 0,
(53)
dx 3s
has a solution consistent with this assumption,
ps
xs = (1 + i )1/2 .
(54)
2
The corresponding exponent is
i
1
4 h
(xs , p) =
(1 + ips/2)3/2 1 (1 i)(ps)3/2 ,
3s
3s
(55)

where the last equality is valid for ps  1. To evaluate the

integral to P decimals precision one must choose p = 2/h
so that (xs , p) = P log 10 (or smaller). For very large P
one may use the last approximation in equation (55) to find
h

2
(3 log 10)

2/3

s1/3 P 2/3 1.73 s1/3 P 2/3 ,

(56)

otherwise it is simplest to determine h numerically. Further,

to evaluate the integral to P decimals precision one must add
contributions from the (positive) x-range where
2

(x)2 e2(x1)

(x+2)/3s

> 10P .

(57)

I.e., we must take xmax to be the largest solution to the equation

(x 1)2 (x + 2) = (3/2) log 10 sP . For very large P the
solution becomes

1/3 1/3 1/3
xmax 32 log 10
s P
1.51 s1/3 P 1/3 .
(58)
Hence the number of required function evaluations again grow
asymptotically for large P like
xmax
0.87 P,
(59)
M +1=
h
cf. equation (49).

E0

E100

P = Pobt Pest

EWW

4
2

8
6
4
2
3
2
1
0
40

200

500

Pest

1000

Fig. 4 These figures display the difference between the actually

obtained precision Pobt and the estimated one Pest , for some
normalization integrals, with precisions measured in decimal digits.
The cases considered are (i) the ground state wavefunction of the
x4 -potential (E0 ), (ii) the 100th excited state wavefunction of the
x4 -potential (E100 ), and (iii) the ground state wavefunction of the
(x2 1)2 -potential (EWW ) with s = 1/100.

ACKNOWLEDGEMENT
This work was supported in part by the Higher Education

Commision of Pakistan (HEC) and the SMAFORSK

program of the Research Council of Norway. The calculations
were performed using the CLN[7] high-precision numerical
library.
R EFERENCES
[1] A. Mushtaq, A. Noreen, K. Olaussen, and I. verb, Very-high-precision
solutions of a class of Schrdinger type equations, Computer Physics
Communications, (in press) doi:10.1016/j.cpc.2010.12.046 (2011)
[2] J. Zinn-Justin, Expansion around instantons in quantum mechanics,
J. Math Phys. 22, 511 (1981)
[3] C.M. Bender, K. Olaussen and P.S. Wang, Numerological analysis of the
WKB approximation in large order, Physical Review D16, 1740 (1977)
[4] A collection is given in M. Abramowitz and I.A. Stegun, Handbook of
mathematical functions, sect. 25.4, Dover Publications, (1970).
[5] Ibid, equation. 25.4.4
[6] . Tafjord, Master thesis, Institutt for fysikk, NTH 1994.
[7] B. Haible and R.B. Kreckel, CLN Class Library for Numbers,
http://www.ginac.de/CLN/