The integrations occur frequently in the inner product of two wave functions or an expan-sion in orthogonal basis. Its specific and non-uniform grids are also used in the generalized-pseudo spectral method throughout the numerical calculations in the present research.
To numerical calculate these integrations, the quadrature and specific grids are necessary.
Gauss-Legendre quadrature and Gauss-Legendre-Lobatto quadrature are introduced in this section.
2.1.1 Gaussian Quadrature
For the numerical integration, Gaussian quadrature has the highest accuracy in theory, but its grids and coefficients are hard to calculate. The general formula of Gaussian quadrature is
∫b
a
w(x)f (x)dx =∑
i
Wif (xi) (2.1.1)
where w(x) is the weighting function, xi are specific grids called Gauss nodes or Gaus-sian grids, and Wi are the weights on the Gaussian grids. When the interval [a, b] and weighting function w(x) are given, the Wi and xi are determined.
There are many typical forms for the interval [a, b] and weighting function w(x) derived from different orthogonal polynomials in (2.1.1). For a example, the interval [−1, 1] and weighting function 1/√
1− x2 in common use are derived from Chebyshev polynomials of the first kind.
2.1.2 Gaussian Grids and Weights
We only chose the [a, b] and w(x) based on the Legendre polynomials throughout the present research. For Legendre polynomials, the valid interval is [−1, 1]. The weighting function w(x) is a constant as follow
w(x) = 1 (2.1.2)
In practice, the unity lets us avoid calculating the weighting function for the integral function. With the given interval [−1, 1] and weighting function w(x), there are three unique sets of Gaussian grids xi and weights Wi. We used two of them in calculation.
The first set is the roots of Legendre polynomial PN(x), x = cos θ, with degree N . The roots can be numerically calculated by the recurrence relations, (A.0.3) and (A.0.4), and Newton’s method. The roots of PN(x) form following collocation points,
xi ∈ {x1, x2, . . . , xN−1, xN}. (2.1.3) In this thesis, the set (2.1.3) of points is called Gauss-Legendre grids, abbreviated to GL grids. They are non-uniformly distributed over the abscissa and symmetric (or antisymmetry) to zero point. The relation between the polynomial degree N and the number of GL grids NGL is
NGL = N (2.1.4)
The weights Wi for GL grids are (Abramowitz and Stegun, 1972, p.887)
Wi = 2
(1− xi)[PN′ (x)]2 (2.1.5)
The domain of GL grids is [−1, 1]. It is suitable for expressing the polar angle in spherical coordinate in discreteness without mapping.
The other set is the roots of the first derivative of Legendre polynomial PN′ (x), x = cos θ, with degree N . The roots can be numerically calculated by recurrence relations, (A.0.4) and (A.0.5) and Newton’s method. In addition, it includes the end points of interval [−1, 1]. The roots and end points form the collocation points,
xi ∈ {x0 =−1, x1, x2, . . . , xN−1, xN = 1} (2.1.6) where x0 and xN are end points, and x1,· · · , xN−1 are the roots of PN′ (x). In this thesis, the set (2.1.6) is called Gauss-Legendre-Lobatto grids, abbreviated to GLL grids.
Their distribution is the same as GL grid. The relation between the polynomial degree N and the number of GLL grids NGLL′ is
NGLL′ = N + 1. (2.1.7)
The weights Wi for GLL grids are (Abramowitz and Stegun, 1972, p.888)
Wi = 2
N (N + 1)[PN(xi)]2. (2.1.8) Because of the inaccuracy close to the end points of interval, the subset of GLL grids was used more often than the GLL grids in actual calculation. The subset excluding the end points, x0 and xN, only has the roots of PN′ (x),
xi ∈ {x1, x2, ..., xN−1}. (2.1.9) The relation between the polynomial degree N , the number of all GLL grids NGLL′ , and the number of this subset NGLL is
NGLL = N− 1 = NGLL′ − 2. (2.1.10)
2.1.3 Mapping Funciton
The GLL grids xi of which valid domain is [−1, 1] can not express the radius r of which valid domain is [0,∞] in spherical coordinates. The connection of different domain is called mapping. To achieve mapping, it requires a function which is also called mapping function.
Between the two intervals, x∈ [−1, +1] and r ∈ [0, ∞], there are two nonlinear map-ping, algebraic mapping and exponential mapmap-ping, to use. We only used algebraic mapping in the present research because the algebraic mapping is more accurate than exponential mapping in practice. The algebraic mapping is
r(x) = L1 + x + β
1− x + α (2.1.11)
where α and β are constants to determine the boundaries, rmin and rmax. With x =−1, r(−1) is the boundary rmin,
rmin = L β
2 + α. (2.1.12)
With x = 1, r(1) is the boundary rmax,
rmax= L2 + β
α . (2.1.13)
When α = β = 0, the interval [−1, 1] can be map to the interval [0, ∞].
The values of two constants α and β in (2.1.11) should be noted. The constant β may not be used and let itself be zero.
β = 0 (2.1.14)
The other constant α should be added and let itself be non-zero to avoid divergence of r(x) at the last GLL point, r(1). r(x), rmin, and rmax with the given β are redefined as
r(x) = L 1 + x
1− x + α (2.1.15)
rmin = 0 (2.1.16)
rmax = L2
α (2.1.17)
The other constant L in mapping functions, (2.1.11) and (2.1.15), is called length scale. L affects the distribution of all collocation points, and it is important to accurately calculate. Figure 2.1.1 shows how to change the distribution of 22 GLL grids by mapping function with various L in a fixed range [0, rmax] = [0, 100]. There are nine horizontal lines for nine distributions by mapping function with nine length scales.
0 20 40 60 80 100
0 10 20 30 40 50 60 70 80 90 100
L=lengthscale
r
L=1 L=5 L=10 L=20 L=25 L=40 L=50 L=80 L=100
Figure 2.1.1: Influence of various length scale L by mapping function on grid-point dis-tribution
L = 1 does not change the distribution. The distribution shows that the most grids are close to the origin, r = 0, and very dense near the origin, r = 0. It can not correctly describe the shape of a function on the middle or far grids, and it may lead to inaccuracy calculation when L = 1. On increase of L, the influence of L on the distribution decreases gradually, and the distribution is changed a little. When L is over a value, the grids start to approach to the rmax.
In practice, the length scale L depends on the variation of function. With the adapted length scales for different problems, the grids can describe the functions well. In the present research, the most functions are the wave functions in stationary state. The wave functions oscillate near the nucleus and decay quickly, so we multiplied rmax/5 or rmax/10 to the mapping function r(x) as a length scale.
Map the interval [−1, 1] to interval [0, rmax] by the mapping function (2.1.15), and the Gaussian quadrature requires a variable change,
r′(x) = d
dxr(x) = L 2 + α
(1− x + α)2, (2.1.18)
in term of x to work. The Gaussian quadrature with interval [−1, 1] and weighting function w(x) = 1 becomes
r∫max
0
f (r)dr =∑
i
Wif (r(xi))dr(xi)
dx . (2.1.19)
where xi are GLL or GL grids, Wi are weights on GLL or GL grids.
More mapping functions that map finite interval to semi-infinite interval or infinite interval can be found in Canuto et al. (2007,§ 8.8 )