Numerical approaches for the analysis of physical systems can be classified into two prominent categories: the grid and basis set approaches [1] [2] [3] [4] [5] [6] [7] . The grid based methods, naturally represent the potential energy in the chosen grid range, while approximating the kinetic energy. These are commonly used in solving differential equations in a large variety of fields in science and technology. [8] [9] [10] [11] . The basis set methods are equivalent to the use of an approximate representation of the identity operator. As a result, they are accurate in a regime where the basis set represents the system states reliably.
A disadvantage of the grid based approach, in comparison to the basis set approach, is the lack of a variational principle [12] . Commonly, the calculated spectrum of the Hamiltonian under study converges from below to the exact values with increasing grid density. This typical characteristic behaviour prevents the estimation of the error based on minimization of the total energy.
Reference [13] introduced a novel Laplacian representation that provides a rigorous upper bound estimate of the true kinetic energy and its properties were illustrated by an analysis of the quantum harmonic oscillator.
The purpose of this paper is to show that upper bounds to the spectrum of any given Hamiltonian can be obtained without the modifying the Laplacian representation (as was achieved in Ref. [13] ) and by by using the same set of coupled equations as are used in the standard FDM that gives lower bounds to the spectrum. These two FDM calculations can then be combined, utilizing the upper and lower bounds to obtain a more accurate spectrum of the Hamiltonian without increasing the numerical efforts.
The strategy of the paper is as follows. First, we describe the two FDMs procedures which produce the lower and upper bounds to the spectrum of the Hamiltonian under study. We then show how the spectrum can be estimated with improved accuracy. Next, we present the numerical results, obtained for several problems, for which the discrete spectrum is analytically known. Finally, we conclude, emphasizing the generality of our approach.
The Hamiltonian operatorĤ is represented by a N ×N dimensional matrix,
in a discretized position basis. The N × N matrix T represents the kinetic energy operator, where m = 2j + 1; j = 1, 2... grid points are used to evaluate the Laplacian (second-order derivative). The maximum number of grid points is N max = N , such that,
is the size of the box which discretizes the kinetic energy spectrum. The discretized spectrum of the kinetic energy operator is given by E n = c(n/L max ) 2 where n ∈ N, and the proportionality constant c is problem dependent. In the solution of the time-independent Schrodinger equation (TISE), utilizing the FDM, the proportionality constant equals to 2 /(2µ) where is the reduced Planck's constant and µ is the particle mass. For the sake of clarity we give below a short description of the derivation of the kinetic energy matrix, calculated using m = 2j + 1 grid points. Consider a 1dimensional evenly spaced grid made up of N nodal points with total length L. The spacing between adjacent nodal points is ∆x = L N −1 . Here x i and ψ i denote the coordinate of the i th nodal point and the approximated wavefunction value at this point, e.g. ψ i = ψ (x i ). We wish to approximate the second order derivative of ψ (x) at x = x i . For this purpose we write the truncated Taylor series expansion around x i using n nodal steps, explicitly written as
where n ∈ [−j, j]. This results in a linear system of equations which relates the vector of the nodal values of the function ψ and the vector of its derivatives ψ (D)
.
(3) Notice that the first (resp. last) element of ψ is ψ i−j (resp. ψ i+j ). The nodal values and its derivatives are related through the matrix A m×m , with elements A n+j+1,k+1 = n k /k!:
By inverting equation (4), we isolate the second order derivative, which is proportionate to the third element of ψ (D) . This leads to a linear combination of the nodal values, with weights w n = A −1 3,n+j+1 . The derivative is then explicitly written as
This relation determines the matrix elements of the (2j + 1)-diagonal matrix T whose elements are given by
while 0 otherwise and i = 1, 2, . . . , N . Table I gives the coefficients for different values of m. We emphasize that the coefficients are symmetric w n = w −n ; this entails that the kinetic energy matrix is symmetric and real, e.g. positive definite. Such method for the approximation of the second order derivative shows significant effect on the eigenvalues of the kinetic energy operator. Figure 1 demonstrates that for a fixed value of N , increasing m provides a curve which approached the exact solution from below, and hence the value of m has a major influence on the derivation of the upper and lower bounds within a specific potential. (5) for different values of m. Here we present only w n for n ≥ 0 as the weights remain symmetric for every n, e.g. w n = w −n . This approach (standard FDM) bounds the spectrum of the Hamiltonian under study from below, for sufficiently large box a number of grid points. In this method the number of grid points N determines the griddifference, ∆x(N ) = L max /N . The plot of the eigenvalues of the N × N matrix T in increasing order, provides a curve which approaches the parabolic function y(n) = c(n/L max ) 2 from below, see Fig 2. As N approaches the maximal value N max the deviation of the n first lowest eigenvalues of T from an exact parabolic behavior, converges to the numerical accuracy of the computations, and are therefore considered as numerical exact results.
Alternatively, the FDM can also provide an upper bound to the spectrum of the Hamiltonian. When the grid-difference is held fixed, δx ≡ ∆x (N max ), while increasing the box-size with N : L (N ) = δxN , for a properly defined range, the numerical result converges from above towards the exact spectrum. In this case, as seen in Fig 2, by the values from the lowest value to the largest one, provide a curve which approaches the parabolic function y(n) = c(n/L max ) 2 from above. As N approaches the maximal value N max , the obtained spectrum approaches the same converged numerical result of the kinetic energy, leading to the converged Hamiltonian spectrum.
The convergence from above can be understood from simple observation that the eigenenergies of a particle confined to a restricted space, with no further external force scale as 1/L 2 . Therefore, as the box size L increases 1/L 2 approaches 1/L 2 max from above, which is taken as the numerically exact value.
To demonstrate how the two FDM schemes can be combined together to evaluate the system spectrum, we compare the FDM results to the analytical solution for two cases: the Harmonic and Rosen-Morse potentials. The Harmonic potential, V HO (x) = 1 2 µω 2 x 2 , includes an infinite number of bounded states with energies E n = ω n + 1 2 , where n = 0, 1, 2, ..., µ is the particle mass and ω is the oscillator frequency. In contrast, the Rosen Morse has finite number of bounded states n max with energies
with n ≤ n max , and the potential is of the form V RM = −V 0 /cosh 2 (ax) [14] . The calculation is performed by the following procedure: First, we evaluate the maximum box size L max uti- lizing a semi-classical approximation. The semiclassical bound state function is well described when box quantization condition is imposed on the quantum solution,
where the eigenenergies of interest lie in the range [min x (V (x)) , E max ] This evaluation is equivalent to employing the WKB method in order to recast the wavefunction in an exponential form [15] [16] [17] . This approximation is valid for large action relative to and smooth potentials, nevertheless, it produces a sufficient evaluation for L max . In our calculations we take r ≈ 10 −7 .
Next, we calculate the eigenenergies with a varying number of grid points N by the two FDM schemes. The standard procedure (L = const), typically, produces a lower bound while keeping the grid density constant with increasing grid size gives an upper bound to the spectrum. This can be observed in Figs. 3  and 4 , which present the energy error, error (E n ) = E numerical n − E exact n /|E exact n |, as a function of quantum number n for the harmonic and Rosen-Morse potentials. The two cases demonstrate the varying convergence behaviour. In the case of the harmonic potential, the standard FDM shows a faster convergence from below relative to the constant grid density method. In contrast, the later method shows a rapid convergence for the Rosen-Morse potential from above. This demonstrates the utility of applying both methods, and combining them to evaluate the exact spectrum.
This method is advantageous as it satisfies a quasivariational principle. The exact ground state eigenvalue is embedded in between the numerical result obtained by the FDM and by the minimum of the potential well. associated with the n-th exact excited state is bounded from above by the n-th numerical eigenvalue, and from below by the (n − 1)th numerical eigenvalue.
Combining the lower and upper bounds, we estimate the spectrum of the Hamiltonian with an improved accuracy. The exact n'th eigen energy is given by the weighted average:
for every n (specific values depend on the potential) where E U/L n are the upper and lower bounds of the energy by the two FDMs that respectively provide upper and lower bounds to the spectrum. The question is how to estimate the weight factor W n . In Fig. 5 we plot W n vs n for the HO and Rosen-Morse potentials. As one can see the weight factors show a smooth functional behavior that enable the predictions of W n for large n from the values obtained.
The two convergence behaviours can be understood by recalling the Hylleraas Undheim MacDonald (HUM) variational principle for excited states [18] [19] [20] . The upper bound results from the fact that the Hamiltonian matrix for N basis functions is included in the larger matrix constructed by M > N basis functions. The FDM we utilize here, for which the box size L is linearly proportional to N , is similar to applying a basis set of sinc functions, φ n (x) = sin(x − x n )/(x − x n ) where x n − x n−1 = ∆x. In contrast, the standard FDM where the box size L is held fixed and ∆x is linearly proportional to 1/N is equivalent to the use of sinc basis functions where x n = L max /N and the N × N Hamiltonian matrix is not included in the larger M × M Hamiltonian matrix and therefore HUM proof does not hold. In such a case the diagonalization of the Hamiltonian matrix will provide a spectrum that will converge from below and not from above as usual.
We emphasize that the analogy between the basis set (utilizing sinc basis functions) and the non-standard FDM method is not perfect, since the matrix representing the kinetic energy is not exact. As a result, the convergence from above is not bounded by a variational principle. Nevertheless, if an appropriate grid size and density are chosen a convergence from above is expected, leading to a quasi-variational principle. As pointed out here it is worth while to calculate the upper and lower limits to the spectrum as have been demonstrated here and find out which method converged faster to the exact spectrum.
We thank Ronnie Kosloff for fruitful discussions. This research was supported by the Adams Fellowship Program of the Israel Academy of Sciences and Humanities. The work is based on a project given in a course on Numerical Methods in Quantum Mechanics given at the Technion by ZOOM during the Coronavirus lockdown.