Finite-element DVR (FEDVR)

dense matrices. Dense matrices are very troublesome with regard to a parallel implementation of the algorithm (cf. page75). For this reason, the LobattoDVRmethod will be generalized in the following section in order to obtain sparse matrix representations for differential operators even in the1dcase.

VI.4. Finite-element DVR (FEDVR)

The DVR discretization techniques described in the previous sections, Secs. VI.2 and VI.3, provide a diagonal representation of local operators but a dense representation of differential operators in1d. One can achieve a sparse representation for both types of operators by con-structing basis functions that have a small support compared to the domain to be discretized.

This procedure is a fundamental idea behind the well-knownfinite element method(FEM) (see, e.g., [314]). In the following, this concept of theFEMwill be emulated to improve theDVR method by ensuring a sparse representation of the Hamiltonian. One proceeds as follows. After dividing the domain (interval) into a collection of small subdomains (subintervals), one applies the LobattoDVRon each of the subdomains while taking the communication of data between the subdomains into account by using suitable boundary conditions⁵. This method is employed in Refs. [4,223,232,233,255,318,342–345] for1dand3dsystems. Accordingly, this technique will henceforth be calledfinite-element discrete variable representation(FEDVR). As the concept of theFEMis considerably more widespread than the DVR method, this section will be far less extensive than the previous sections. The recipe for constructing theFEDVRbasis will be presented in conjunction with the resulting nonzero structure of the matrix which represents the Hamiltonian. This is in general a very helpful information for testing the implementation.

Let the interval of interest[a, b]be arbitrarily divided inton_e subintervals[xⁱ, xⁱ⁺¹]where x⁰ = aand x^ne−1 = b. Further, letM_i := ¹₂(xⁱ⁺¹+xⁱ)be the center ands_i = ¹₂(xⁱ⁺¹−xⁱ) be the scaling factor of thei-th interval with respect to the original interval[−1,1]. On each subinterval, one constructs an_g dimensional LobattoDVRfrom Eq. (VI.16) via scaling and translation

One can now easily define the continuous and normalizedFEDVRbasis functions in terms of theϕⁱ_m. For(i, m)∈I := ({0, . . . , ne−1} × {0, . . . , ng−2})\ {(0,0)}, one writes functions localized about the same value of adjacent intervals are merged into one continuous basis function. TheDVRfunctions localized at the boundary of the global interval are dropped, which complies with Dirichlet boundary conditions.

5Recall that the LobattoDVRwas constructed for an easy implementation of boundary conditions.

86 VI. Numerical methods

−5 −4 −3 −2 −1 0 1 2 3 4 5

x 0.0

0.5 1.0

Figure VI.4.:SmallFEDVRbasis forn_g= 3andn_e= 10in the interval[−5,5]with equally spaced boundaries of subintervalsxⁱ ∈ {±1,±3,±5}. The basis comprises 9 functions, which have zeros at each grid point with one exception. Each basis function is identified with the grid point at which it has a nonzero value:χ⁰₁(blue) is localized aroundx=−4,χ¹₀(green) aroundx=−3, etc. This small basis is not used in practice; it is shown for illustrative purposes only.

An example of a very small FEDVR basis is depicted in Fig. VI.4. The interval size is constant, and thus the scaling factor is triviallys_i= 1for all intervals. Therefore, theFEDVR basis essentially consists of the shiftedDVRbasis functions from Fig.VI.2a. Adopting the term ofBalzer et al. (2010)[232], the mergedFEDVR basis functions localized at the boundaries of the intervals are called bridge functions. In this special case, the bridge functions are, according to Eq. (VI.22), scaled in value by ^√1₂ with respect to the DVRbasis Fig. VI.2ato ensure normalization of the resulting bridge function. Moreover, one shall keep in mind that the bridge functions are obviously not differentiable at the associated grid point. This has further implications: firstly, regarding the representation of functions within theFEDVR, differentiable functions in the Hilbert spaceF := span(B)spanned by theFEDVRare not arbitrary linear combinations of basis functions. This means that one might be restricted to subspaces such as F ∩C¹([a, b])for certain purposes. Secondly, regarding the representation of operators within the FEDVR, differential operators are only well-defined on the subspace of “well-behaved”

functions inF. Consequently, differential operators do not have a unique continuation onF. This ambiguity is, however, rather commonplace than irritating if one is accustomed to thefinite difference method(FDM).

Analogous to the previous sections, Secs.VI.3andVI.2, one can define a quadrature rule and likewise an inner product onF by making use of the following discreteFEDVRmeasure

D = X

(i,m)∈I

wⁱ_mδ x−xⁱ_m

, (VI.23)

wherexⁱ_m :=s_ix_m+M_iare theFEDVRgrid points andw_mⁱ areFEDVRgrid weights defined bywⁱ_m := (χⁱ_m(xⁱ_m))⁻².

The corresponding quadrature formula can be obtained from Eq. (VI.17) by splitting the inter-val into subinterinter-vals and transforming all arising integrals into integrals of the form Eq. (VI.17).

As usual, the basis functions from Eq. (VI.22) are orthonormal with respect to the inner prod-ucth,i^D. Before presenting the utilized standard representation, which determines the nonzero

87

VI.4. Finite-element DVR (FEDVR)

structure of matrices that represent operators, the canonical mapping of the two indices(i, m)is given as

M:I → {0, . . . , d−1}, (i, m)7→i(ng−1) +m−1, (VI.24) whered:= dim(F) = #I =n_e(ng−1)−1is the dimension ofF. Note thatMis obviously a bijection. The standard representationφB can be now written as

φB :F →C^d, φB(g) = X

(i,m)∈I

g(xⁱ_m)

χⁱ_m(xⁱ_m)eˆM((i,m)), (VI.25) whereeˆ_k is thek-th element of the standard basis ofC^d.

Again, analogous to Eq. (VI.19) this mapping may be continued onHas given in Eq. (VI.25) if required. Note that for any functiong ∈ H, even ifg /∈ F, the values on the grid points are reconstructed properly. This condition was also fulfilled by theDVR.

To appreciate the sparse matrices, the matrix representation of the operator for the first derivativeDis analyzed in a little more detail. For further information and representation of the Hamiltonian, see App.B. The derivation matrixD∈C^d×dis characterized byD=φB◦∂x^∂ ◦φ⁻¹_B on the subspaceF ∩C¹([a, b])and adapted to the Dirichlet boundary conditions. The explicit calculation is given in App.B.3.1. It turns out thatDcan be written in terms of another matrix De ∈ C^d×d so thatD =BDBe whereB = diag not depend on the interval sizes. To illustrate the structure of the sparse matricesDandDe (the nonzero structure is identical in both cases),De is given below for a very small system ofn_g = 4 andn_e = 6:

88 VI. Numerical methods

where D^loc_mk =√ wkwm

Z1

−1

ϕk(x)ϕ⁰_m(x)dx. (VI.27) In general, the nonzero diagonals ofD_mk^loc (cf. the derivation operator in theDVR) cancel each other out, thus generating a skew-symmetric matrix due to Dirichlet boundary conditions. For test purposes, there is an instructive example for differential operators in App.B.3.3. It is also shown that in this special case, theFEDVRcan be reduced to a nontrivialFDM, which may be more intuitive for the reader.

The matrix representing the kinetic-energy operator can be chosen to have the same nonzero structure as the derivation operator in Eq. (VI.26) except for the diagonal. This is shown in App.B.3.2, expressed in simpler terms as in the original publications [232,233] without loss of generality. It is worth mentioning that the matrix for the kinetic energy isnotutilized in the form

1

2D^TDor−¹₂D²because, first, there is a sparser representation and, second, problems due to finite floating-point precision can be mitigated (see also Sec.VI.5).

The FEDVR provides sparse representations of the relevant operators in a very flexible way, allowing for an arbitrary enhancement of the local accuracy (increasen_g) at the expense of sparsity. On a side note, the dense matrices appearing in the orbital form of the TDHF equation Eq. (II.2) due to the nonlocal exchange operator do not necessarily mean that the DVR is preferable over theFEDVR because a dense matrix is involved anyway⁶. Common time-propagation algorithms may distinguish between terms with different properties, so one can still benefit from the sparsity (cf. the ensuing section).