Spectral: Solving Schroedinger and Wheeler–DeWitt equations in the positive semi-axis by the spectral method

Computer Physics Communications 185 (2014) 380–391

Contents lists available at ScienceDirect

Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc

Spectral: Solving Schroedinger and Wheeler–DeWitt equations in the positive semi-axis by the spectral method✩ E.V. Corrêa Silva a,∗ , G.A. Monerat a , G. de Oliveira Neto b , L.G. Ferreira Filho a a

Depto. de Matemática, Física e Computação, Faculdade de Tecnologia, Universidade do Estado do Rio de Janeiro, Rod. Pres. Dutra km 298 s/n, 27537–000, Resende, RJ, Brazil b

Depto. de Física, Instituto de Ciências Exatas, Universidade Federal de Juiz de Fora, 36036–330, Juiz de Fora, MG, Brazil

article

info

Article history: Received 25 July 2012 Received in revised form 29 August 2013 Accepted 3 September 2013 Available online 18 September 2013 Keywords: Schroedinger-like equation Wheeler–DeWitt equation Galerkin spectral method

abstract The Galerkin spectral method can be used for approximate calculation of eigenvalues and eigenfunctions of unidimensional Schroedinger-like equations such as the Wheeler–DeWitt equation. The criteria most commonly employed for checking the accuracy of results is the conservation of norm of the wave function, but some other criteria might be used, such as the orthogonality of eigenfunctions and the variation of the spectrum with varying computational parameters, e.g. the number of basis functions used in the approximation. The package Spectra, which implements the spectral method in Maple language together with a number of testing tools, is presented. Alternatively, Maple may interact with the Octave numerical system without the need of Octave programming by the user. Program summary Program title: Spectral Catalogue identifier: AEQQ_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEQQ_v1_0.html Program obtainable from: CPC Program Library, Queen’s University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 20417 No. of bytes in distributed program, including test data, etc.: 2149904 Distribution format: tar.gz Programming language: Maple, GNU Octave 3.2.4 Computer: Any supporting Maple Operating system: Any supporting Maple RAM: About 4 Gbytes Classification: 1.9, 4.3, 4.6. Nature of problem: Numerical solution of Schrödinger-like eigenvalue equations (especially the Wheeler–DeWitt equation) in the positive semi-axis Solution method: The unknown wave function is approximated as a linear combination of a suitable set of functions, and the continuous eigenvalue problem is mapped into a discrete (matricial) eigenvalue problem Restrictions: Limitations are due to memory usage only

✩ This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/ science/journal/00104655). ∗ Corresponding author. Tel.: +55 2433813889. E-mail addresses: [email protected], [email protected] (E.V. Corrêa Silva).

0010-4655/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cpc.2013.09.007

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

381

Unusual features: The package may not work properly in older versions of Maple, due to a bug in that CAS; for that reason an interface with the GNU Octave system is provided, requiring no user intervention or Octave programming during calculations Running time: Seconds to hours, depending on the number of basis functions used and on the complexity of the potential used © 2013 Elsevier B.V. All rights reserved.

1. Introduction The study of differential equations is of undoubted importance in the modeling of physical systems, from quantum to cosmological phenomena. The Schroedinger equation, for instance, lies at the basis of the mathematical formulation of quantum mechanical theory, and theoretical developments regarding its solution under various scenarios are still under research. The Wheeler–DeWitt equation [1] (the form of which generalizes the Schroedinger equation) is defined on the positive semi-axis and is of key importance to contemporary attempts at a quantum description of the Early Universe (the so-called Planck Era). Several authors have applied Schutz’s variational formalism [2,3] to obtain the Wheeler–DeWitt equation for Friedmann–Robertson–Walker (FRW) spacetimes with a perfect fluid for their matter content [4–13]. The Galerkin spectral method [14] is a well-known numerical method of solution of differential equations, and has also been used in cosmology [8,9,11]. The unknown solution (rather, an approximation to it) is written as a linear combination of a suitable number of basis functions, and equations for the coefficients of such expansions are obtained. The eigendifferential problem is thus mapped into an approximate matricial eigenvalue problem, which can be solved, in its turn, by linear algebra methods. The choice of basis functions is a matter of convenience and easiness of manipulation. Generally speaking, the approximate eigenvalues and eigenfunctions obtained by any numerical method may show different degrees of accuracy. In the context of isolated quantum systems, a common touchstone is the conservation of the norm of the wave function constructed with a given set of eigenvalues/eigenfunctions. However, some other questions may be asked, for that matter. (1) How well does the set of approximate solutions (eigenvalues and eigenfunctions) satisfy the eigenvalue equation? (2) How orthonormal are the approximate eigenfunctions obtained? (In fact, the orthogonality impacts directly on norm conservation.) (3) How can we compare the solutions for two slightly different computational problems (e.g., for different values of the physical parameters of the model, or for different numbers of basis functions used in the approximation) and draw conclusions on the stability of our results? In this work we describe the package Spectral, which implements the spectral method for Wheeler–DeWitt and Schroedinger equations, and also provides tools that address the aforementioned questions. The package has been implemented in the Maple language but may also use numeric resources of the GNU Octave software [15]. The spectral method is summarized in Section 2. Section 3 offers more details on the checks outlined above. In Section 4 the usage of Spectral is described. (For more detailed explanations on command usage, the reader is referred to the user manual, distributed with the package.) In Section 5 we present an example with the Wheeler–DeWitt equation for FRW spacetimes with stiff matter perfect fluid and a negative cosmological constant [11]. 2. Outline of the spectral method We will consider Schroedinger-like time-independent equations of the form

 −

d2 dx2

+ f (x)



ψ(x) = E g (x)ψ(x)

(1)

for x ∈ [0, ∞), in which f (x) and g (x) are known functions of the spatial variable x, ψ(x) is the spatial part of the wave-function of a quantum system, and E is the energy eigenvalue. The Schroedinger equation is obtained as a special case, by letting g (x) = 1. Restricting the domain to x ∈ [0, L] is equivalent to placing two infinite impenetrable barriers at x = 0 and x = L, thus imposing the boundary conditions ψ(0) = ψ(L) = 0. The inner product

(Ψ , Φ ) =

∞



dx Ψ ∗ (x) w(x)Φ (x)

(2)

0

involves a suitable weight function w(x), which must be chosen so that the Hamiltonian operator

Hˆ ≡ −

d2 dx2

+ f ( x)

(3)

be self-adjoint [5], that is,

(φ, Hˆ ψ) = (Hˆ φ, ψ).

(4)

Moreover, the weight function plays a fundamental role in the notion of orthogonality of eigenfunctions corresponding   to different values of the eigenvalue E in (1). From that equation and boundary conditions of the type limx→0 C1 ψ(x) + C2 ψ ′ (x) =

382

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

limx→∞ C1 ψ(x) + C2 ψ ′ (x)



(Ei − Ej )

= 0, in which (C1 , C2 ) are constants, one can easily obtain [16] that



∞



dx ψi∗ (x) g (x) ψj (x) = 0,

(5)

0

in which {Ei , Ej } are eigenvalues corresponding to the eigenfunctions {ψi (x), ψj (x)}, respectively. For Ei ̸= Ej , we must have ∞



dx ψi∗ (x) g (x) ψj (x) = 0.

(6)

0

Hence we may identify the l.h.s. of the last equation as the inner product (ψi , ψj ), with a weight function w(x) equal, up to a constant factor, to the function g (x) in (1), so that two physical states for which a measurement of the eigenvalue E is certain to give different results are orthogonal [17]. It is required that the wave-function ψ(x) fall sufficiently fast for large values of the spatial variable x; hence for the sake of numerical computations we may reduce the initially unbounded domain for x to a bounded domain, say, a ∈ (0, L), in which L ∈ ℜ is suitably chosen. Let ψp (x) be the eigenfunction related to the p-th eigenvalue Ep . The Galerkin spectral method (SM) [14] for the solution of the eigenvalue equation (1) starts with the choice of an orthonormal basis of functions for the expansion of ψp (x). A possible choice is that of sine functions1 ; we may write ψp (x) as

ψp (x) ≈

N 

A(np)



2

sin

L

n=1

 nπ x  L

,

(7)

(p)

in which the coefficients An are yet to be determined, and a finite number N of basis functions has been chosen. The same basis is used for the expansion of the terms in Eq. (1), N 

f (x)ψp (x) ≈

Fn(p)



2 L

n=1 N 

g (x)ψp (x) ≈

(p)



Gn

L

n=1

(p)

The coefficients Fn Fn(p) =

N 

2

sin

sin

 nπ x  L

 nπ x  L

,

(8)

.

(9)

(p)

(p)

and Gn are related to the coefficients An of Eq. (7) by

Φn,m A(mp) ,

(10)

Γn,m A(mp) ,

(11)

m=1

G(np) =

N  m=1

in which

Φn,m =

2

L



sin

L

 nπ x  L

0

f (x) sin

 mπ x  L

dx

(12)

dx.

(13)

and

Γ n ,m =

2 L



L

sin

 nπ x 

0

L

g (x) sin

 mπ x  L

Introducing the approximate expansions (7), (8) and (9) into Eq. (1) we obtain

 nπ 2 L

A(np) +

N  m=1

(p) Φn,m Am =E

N 

Γn,m A(mp) .

(14)

m=1

In matricial form, this generalized eigenvalue equation2 may be written as

∆ A(p) = Ep Γ A(p) ,

1 It is convenient, but not mandatory, that the basis functions satisfy the boundary conditions at x = 0 and x = L. 2 In Spectral, this equation is handled either by the Maple package LinearAlgebra or by the Octave function qz.

(15)

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

383

(p)

in which Γ is the N × N square matrix with elements given by (13), A(p) is a 1 × N column vector with elements Am and ∆ is an N × N square matrix with elements

∆n,m =

 nπ  2 L

δn,m + Φn,m .

(16)

Hence the eigenvalue problem (1) has been mapped into the approximate discrete problem (15), the solution of which provides the approximate eigenvalues Ep and eigenvectors A(p) , and hence the corresponding approximate eigenfunctions ψp (x) of (7). The larger the value of N the closer our results will be to the exact ones (apart from rounding errors and memory limits, naturally). Once the eigenfunctions and eigenvalues are obtained, the wave package Ψ (x, t ) can be built as a linear combination of the eigenfunctions, and the expected position as a function of time can be calculated as

⟨x⟩ (t ) =

∞



dx Ψ ∗ (x, t ) w(x) x Ψ (x, t ).

(17)

0

The uncertainty of the position, in its turn, is

σ (t ) =

 

x2 − ⟨x⟩2



∞

dx Ψ (x, t ) w(x) x Ψ (x, t ) − ∗

=

2

∞



0

2 1/2 dx Ψ (x, t ) w(x) x Ψ (x, t ) . ∗

(18)

0

3. Eigenvalue and eigenvector checking In addition to the norm conservation, the following checks can be performed.

• Equation checking. Generally, if the eigenvalues and eigenvectors are substituted, under the same number of digits used in spectrum calculations, back into the l.h.s. of the matricial equation (15), the latter is not satisfied exactly, i.e., one does not obtain the vector in its r.h.s. Rather, the difference between the l.h.s. and r.h.s. is expected to be a vector with small non-vanishing components; we may use the largest absolute value of such components as a measure of error. • Eigenfunction orthonormality checking. The eigenfunctions of the Hamiltonian (3) must be orthonormal under the inner product (2), ∞



dx ψn∗ (x) w(x) ψm (x) = δn,m .

(19)

0

Given the expansion of eigenvectors in Eq. (7), the dot product of two eigenvectors A(p) and A(q) of Eq. (15) related to the eigenvalues Ep and Eq can be defined as A(p) · A(q) ≡

N  N 

(q) An(p) Wnm Am

(20)

n=1 m=1

in which the weight matrix W is obtained from the weight function w(x) as Wnm ≡

2 L

L



dx sin 0

 nπ x  L

w(x) sin

 mπ x  L

.

(21)

Hence, the orthonormality condition of eigenvectors reads A(p) · A(q) = δp,q .

(22)

For each eigenvector (or pair of eigenvectors), deviation from normality (or orthogonality) can be used as a measure of error, computed by evaluating the difference between δp,q and the actual value of the l.h.s. of (22). The eigenvectors A(p) are defined up to a scalar factor. In Spectral, we have established the convention that, given an eigenvector (p)

A(p) with components Am , m ∈ {1, 2, . . . , N }, its first non-vanishing component should be positive. (Octave and Maple may choose it otherwise; hence the need for establishing some convention.) The conservation of the norm of the wave function Ψ (x, t ), commonly used as a test for the accuracy of numerical calculations, stems from the orthogonality of the (approximate) eigenfunctions ψn (x). The more accurate the orthogonality, the better the norm conservation. 4. General features of the Spectral package In what follows we show how to set up the package (Section 4.1) and how to use its commands (Section 4.2). A short discussion on limitations of a purely symbolic approach to the implementation of the spectral method is presented in Section 4.3. The package Spectral has been written in Maple language and it can rely on Maple resources alone, but it may also interface with the GNU Octave system [15]; Section 4.4 explains why and how. A few examples of the accuracy of eigenvalues and eigenfunctions are presented in Section 4.5.

384

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

For the installation of Maple and Octave, the reader is referred to their respective documentation [15,18].

4.1. Loading the package The package source code is contained in the file Spectral and should be loaded by the standard maple command read. Supposing that this file is stored in the local working directory,3 this would amount to the following.

> read Spectral: The code has been written and tested for Maple 13, but it should work fine in any compatible version. Please refer to Maple documentation for further information [19]. For conjoint use with the Octave system, the path of its binary file must be indicated to Maple by the command octavedir. For more details, the reader is referred to the user manual, distributed with Spectral.

4.2. Usage This section will give an general description of the commands in the package (for a detailed description, the reader is referred to the user manual, distributed with it), which include: the definition of eigenvalue problems of type (1); the solution (eigenvalues and eigenfunctions) of such problems, with the optional help of Octave resources; the assembly of a normalized wave function and corresponding probability density function, with arbitrarily superposition of eigenfunctions; the calculation of the expected position and uncertainty; the tests described in Section 3; and the comparison of the eigenvalue spectra of different problems (with slightly different values of parameters, for instance). The concept of a module is central to the package implementation. A module is a Maple language structure which works as a container for other variables, akin to the concept of a class in object-oriented (OO) languages. Even though the concepts of classes and objects are not fully realized by the module structure (e.g., encapsulation is implemented and procedure overloading is allowed, but no inheritance mechanism is available), programmers in Maple language can benefit, with limitations, from the OO paradigm when projecting their software. In Spectral, every eigenvalue problem is stored in a module, created by a ‘‘constructor’’ procedure, and that module in its turn contains all data and procedures for its solution and the manipulation of that solution. For more details on the module structure, the reader is referred to Maple documentation [19]. In what follows we give brief examples of the commands provided by the package.4 The column ‘‘CPU time’’ shows the CPU time used by each command of the package.5 The time values have been returned by the standard Maple command time; a null value actually indicates a CPU time of less than 10−3 s. More details on Spectral commands can be found in the user manual, distributed with the package. Desired action The  quantum  harmonic oscillator problem 2

− dxd 2 + x2 Ψ (x) = E Ψ (x) is defined for x ∈ [0, 10].

Approximate eigenfunctions with 100 levels (i.e., corresponding to the linear combination of 100 sine functions of the form shown in Eq. (7)) will be employed. Output files should be named after the label HalfOscillator. The weight function is equal to unity. The problem representation is stored in the module p. Problem data are printed. The path of the Octave executable file is defined, so that Maple can find it. (This is mandatory only if Octave is employed, of course.) Eigenvalues and eigenvectors are calculated by Maple. Eigenvalues and eigenvectors are calculated by Maple, but using Octave resources. The eigenvalue spectrum (or part of it) is shown. In the examples: all; just the lowest one; just the three lowest ones. The (approximate) eigenfunction associated to the first lowest eigenvalue is shown, as a function of x.

Examples of calling sequences

p := Problem(Ffunction = x∧2, Gfunction=1, Variable = x, Length=10, NumberOfLevels = 100, Weigth=1, Label=‘‘HalfOscillator’’)

CPU time (s) 1.154

p:-IsDescribed() octavedir(’’C:/Programs/Octave/bin’’)

0 0

p:-IsSolved() p:-IsSolvedByOctave()

8.299 3.775

p:-ItsEigenvalues() p:-ItsEigenvalues(1) p:-ItsEigenvalues(1..3) p:-ItsEigenfunction(x,1)

0 0.015 (continued on next page)

3 The working directory is defined/recalled by Maple native command currentdir. 4 The character / indicates the continuation of a command line in Maple language. 5 On a VAIO VPCF111FB notebook with 6GB of RAM and Intel i5 CPU with 2.40 GHz, running Windows 7 (64 bits). However, execution times may vary not only with CPU, RAM and operational system, but also with Maple version and may even depend on the current status of the running Maple system.

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

A normalized wave function, in terms of the spatial variable x and the time variable t, is built out of the eigenfunctions. Each eigenfunction is represented by a coefficient as a list. Here, we use 100 coefficients, in which those related to the three lowest energy levels are equal to unity, and the remaining ones are null. The normalization factor of the wave function is calculated automatically, taking into account the weight function of the inner product (2). The probability density for the same situation described above is shown. WaveFunction does not need to be used at all, if what one really wants is the probability density. The probability density function for the same situation described above is plotted as a function of the spatial variable x, for a single instant of time t = 0. ProbabilityDensity does not need to be used before this command, if what one wants is the plot of the probability density, rather than its expression. The probability density function for the same situation described above is animated as a function of the spatial variable x, for the time interval [0, 1]. ProbabilityDensity does not need to be used before this command, if what one wants is the plot of the probability density, rather than its expression. The expected position, Eq. (17), is returned as a function of time, in the form of a procedure. The same wave-function coefficients are used. The position uncertainty, Eq. (18), is returned as a function of time, in the form of a procedure. The same wave-function coefficients are used. The expected position is calculated for the time interval [0, 10], for the same wave package used above, and having the time interval divided into 100 sub-intervals of the same length. Data are saved to a text file named HalfOscillator_Position. Similarly to the command ExpectedPositionIsCalculated, both the expected position and its uncertainty are calculated. Data are saved to the text file HalfOscillator_PositionUncertainty. The data file generated by ExpectedPositionIsCalculated is read and the expected position is plotted as a function of time. The data file generated by ExpectedPositionAndUncertainty is read and the expected position, together with the uncertainty, is plotted as a function of time. The norm of the wave function (represented by its coefficients) as a function of time is returned as a procedure. The spectrum of eigenvalues and eigenvectors is tested (cf. Section 3) and results are saved in three human-readable report files: HalfOscillator_Check_Equations, HalfOscillator_Check_Normality and HalfOscillator_Check_Orthogonality. Data in the aforementioned report files are shown in graphical form. The spectra of two problems p and q may be compared, and the difference in eigenvalues is saved into a file. Data from the files thus created are displayed graphically.

385

p:-WaveFunction(x,t,[1$3, 0$97])

0

p:-ProbabilityDensity(x,t,[1$3, 0$97])

81.573

p:-ProbabilityDensityIsPlotted(0, [1$3, 0$97])

0.686

p:-ProbabilityDensityIsPlotted(0..1, [1$3, 0$97])

37.721

p:-ExpectedPosition([1$3,0$97])

0.015

p:-PositionUncertainty([1$3,0$97])

0.078

p:-ExpectedPositionIsCalculated(/ 0..10,[1$3,0$97], 100)

594.083

p:-ExpectedPositionAndUncertainty/ AreCalculated(0..10,[1$3,0$97], 100)

1535.674

p:-ExpectedPositionIsPlotted()

0.031

p:-ExpectedPositionAndUncertainty/ ArePlotted()

0.047

p:-NormOfWaveFunction([1$3,0$97])

0.016

p:-CheckSolutions:-Numerically()

105.862

p:-CheckSolutions:-Graphically()

3.697

p:-IsComparedTo:-Numerically(q)

0.047

p:-IsComparedTo:-Graphically(q)

0.109

Some differences in precision may be found in results for different versions of Maple, due to differences in the way Maple native procedures deal with the numerical integrals and/or handle linear systems of equations. Occasionally, the command CheckSolutions:-Numerically may issue a warning message like ‘‘warning: scalar product of levels [...] yielded a complex number. I will record just the real part’’. This is due to the way different versions of Maple deal with the numeric integrals involved in the calculation of scalar products. These integrals are calculated by Maple native procedures; sometimes they issue a small spurious imaginary part or even a null float imaginary part (0. I) which is not automatically discarded. By ‘‘small’’, we mean a number which is of order 10−D , in which D is about the number of significant digits used. It gets smaller as D increases; so, we know it is a numerical artifact that can be safely discarded. The same holds for the result 0. I issued by Maple, indicating that such a quantity is null up

386

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

to the working precision, so it should not simply vanish from sight. By the very definition of the scalar product, addressed in the warning message, we know that it is a real number; so, any small spurious imaginary part that eventually shows up can be safely ignored. 4.3. Symbolic × Numeric computation Regarding the implementation of the spectral method, there are advantages but also limitations regarding the use of symbolic computation, more specifically to the use of exact arithmetic. For instance, let p be a problem-module, i.e., the module output by the Problem procedure. The integrals corresponding to the matrix elements Φn,m and Γm,n in Eqs. (10) and (11) are calculated symbolically, whenever possible, by the procedures p:-C and p:-C2, exported by p. It is important that, whenever possible, an exact evaluation of the integrals of the aforementioned equations be carried out. The evaluation of arguments of the sine and cosine functions within the integrands is a delicate issue; transforming those arguments into floats may hinder the evaluation of those integrals. Here, a symbolic calculation has an advantage over a purely numerical approach. With that in mind, the procedure Problem automatically converts all floats within its parameters into rational numbers, so that no spurious approximations are introduced. Let us imagine an example such as

> p:-C(1,1) 50 −3 + 2 π 2 3

π2

which involves the literal constant π . However, the procedure CalculateEigenVectors attempts to solve the (high-order polynomial) eigenvalue equation. Maple can only solve that equation if its coefficients are numeric: they must not contain any symbolic constants such as π — from that point on, one must work with floats. The package allows, in principle, calculations in Maple (i.e., by the command IsSolved) with an arbitrary number of digits, which is another useful characteristic of symbolic systems. Its default value is 15, but it can be set otherwise by the option NumberOfDigits of the Problem command. (More details can be found in the user manual, distributed with the package.) Sadly, our attempts to use more than 15 digits have caused the Maple system to crash and issue the ‘‘kernel connection has been lost’’ error message, due to undocumented bugs in that system (at least for version 13, in which our package has been tested) that are out of our control. Nevertheless, the package will issue a warning message to the user, if they try to use more than 15 digits. Calculations with Octave (by the command IsSolvedByOctave), on the other hand, are always performed with 15 digits: the option NumberOfDigits is ignored. 4.4. Maple and Octave By the time Spectral was first developed and run in Maple 12, we only had the command IsSolved, which relies solely on Maple resources. A number of problems occurred, though. For instance, defining the problem

> p1 := Problem(NumberOfLevels=100, Length=10, Ffunction=-36*x^2 - 12 * Lambda * x^4, Gfunction=12, Variable=x, Weight=1); and then letting Maple try to solve it with

> p1:-IsSolved(); resulted, after a long while, in a dialog box with the error message ‘‘kernel connection has been lost’’ and the puzzling suggestion that the firewall might have something to do with it. This is not a bug of Spectral, but rather a known Maple bug described as ‘‘a logic error in processing’’ [20]. After a series of unsuccessful attempts to circumvent this bug, we decided to write the command IsSolvedByOctave, delegating the calculation of eigenvalues and eigenvectors to Octave. Despite the 15-digit limitation on precision imposed by Octave, it turned out to be a fast and reliable calculating tool. The same problem is solved by

> p1:-IsSolvedByOctave(); in about 5 s. More recent tests of Spectral in Maple 13 indicate that the kernel connection issue seems to have disappeared (at least for our particular situation), and both IsSolved and IsSolvedByOctave are completely operational. We can only hope that newer versions of Maple will be free from such odd behavior. In order to have Maple and Octave in cooperation, the path of the Octave executable file must be indicated to Maple by the command octavedir, as described at the beginning of Section 4.1. For the sake of calculations using Octave, the environment variable Digits, which controls the number of significant digits used, is set to 15. Maple and Octave interaction occurs in five steps, without user intervention (e.g., Octave programming), and it is represented in Fig. 1. Supposing that the problem to be solved has been labeled as foo, in step 1 the procedure IsSolvedByOctave (within the package Spectral, interpreted by the Maple system) generates two data files, foo_Input_FirsMatrix.tmp and foo_Input_SecondMatrix.tmp, corresponding to the matrices ∆ and Γ in Eq. (15). In step 2 Spectral automatically generates the Octave script file foo_OctaveScript.m to be interpreted by Octave in step 3, with the aforementioned data files as input files, for the sake of solving the eigenproblem. In step 4 the output files foo_Output_Eigenvalues.tmp and foo_Output_Eigenvectors.tmp, containing data corresponding to the eigenvalues Ep and eigenvectors A(p) in Eq. (15), are generated by Octave and, in their turn, are read by Maple in step 5. All files are ASCII text files, and they are stored in the current working directory. Essentially, the native Maple commands upon which this interaction scheme has been built are the writeline command (which allows the generation of the Octave script file), the ImportMatrix and ExportMatrix commands (which allow the writing/reading of matrices into/from files, respectively), and the ssystem command (which allows invoking the Octave system from within Maple). The generalized eigenvalue equation (15) is actually solved by the Octave function qz.

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

387

Fig. 1. Schematic representation of the several steps in Maple and Octave interaction (which takes place without user intervention). Data files are produced in step 1, whereas an Octave script is produced in step 2. In step 3 that script is interpreted by Octave and uses the aforementioned files as sources of input data; eigenvalues and eigenvectors are calculated. In step 4 the solution is written to output files, which in their turn will be read by Maple in step 5.

4.5. Accuracy of eigenvalues and eigenfunctions: examples In this section we evaluate the accuracy of eigenvalues and eigenvectors for two systems with well-known exact solutions: the halfharmonic oscillator and the linear potential. Consider the time-dependent Schroedinger equation for the half-harmonic oscillator

−

∂ 2ψ + x2 ψ = E ψ, ∂ x2

(23)

which can be defined by the command line

> p := Problem(NumberOfLevels = 100, Length=14.96, Ffunction = x^2, Gfunction=1, Variable = x): in which 100 levels will be used in calculations with 15 significant digits in the working interval [0, 14.96]. (The upper limit of the working interval has been obtained empirically, in order to minimize errors.) The eigenvalue/eigenfunction problem is solved by either p:-IsCalculated() or p:-IsCalculatedByOctave(). In either case, we may compute the absolute error

  δ En ≡ Enexact − Ennum 

(24)

between the n-th exact eigenvalue Enexact = 4n − 1

(25)

and the approximate eigenvalue

> ExactEigenvalue :=

Ennum .

The procedure

n -> 4 * n - 1;

may represent Enexact ; the errors δ En may be calculated by

> AbsErrorEigenvalue := proc(n) global p; abs(p:-ItsEigenvalues(n) - ExactEigenvalue(n)); end; and shown in graphical form by the procedure

> PlotAbsErrorEigenvalue := proc(xmin,xmax) local i; plot([seq([i,’AbsErrorEigenvalue’(i)],i=xmin..xmax)]); end: For instance, the procedure call

> PlotAbsErrorEigenvalue(1,100); will plot the error δ En for n ∈ {1, 2, . . . , 100}. In our example, we obtain that δ En ≤ 4 × 10−13 for n ≤ 38. We may also define the absolute error on the n-th eigenfunction as

δψn ≡ max ∆n (x),

(26)

x∈[0,L]

in which

  ∆n (x) ≡ ψnexact (x) − ψnnum (x) ,

(27)

ψnnum (x) is the n-th approximate eigenfunction and ψnexact (x) is the corresponding exact eigenfunction ψnexact (x) = 

H2n−1 (x)e−x

2/2

√ ,

22n (2n − 1)! π

(28)

388

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

so that Hp is the p-th Hermite polynomial. A procedure for computing ψnexact (x) would look like

> ExactEigenfunction := unapply(subs(n=2*n-1,sqrt(2)*(1/Pi)^(1/4)* 1/sqrt(2^n * n!) * HermiteH(n, x) * exp(-x^2/2)),[x,n]); The function ∆n (x) can be calculated by the procedure

> AbsErrorEigenfunction := proc(x,n) global p; abs(evalf(p:-ItsEigenfunction(x,n) - ExactEigenfunction(x,n))); end: and shown in graphical form by

> PlotAbsErrorEigenfunction := proc(n,xmin,xmax) local x; global L; plot(’AbsErrorEigenfunction’(x,n),x=xmin..xmax); end: For instance, for the first eigenfunction and x ∈ [0, 14.96] the procedure call

> PlotAbsErrorEigenfunction(1,0,14.96); will plot ∆n (x) in that interval. Its largest value is δψ1 . In our example, we have δψn < 5 × 10−13 for n ≤ 30. Consider yet another example with the time-dependent Schroedinger equation for the linear potential

−

∂ 2ψ + x ψ = E ψ, ∂ x2

(29)

which can be defined by the command line

> q := Problem(NumberOfLevels = 100, Length=18.22, Ffunction = x, Gfunction=1, Variable = x): in which 100 levels will be used in calculations with 15 significant digits in the working interval x ∈ [0, 18.22], the upper limit of which has been determined empirically. The eigenvalues En of (29) are the solutions of the equation Ai(−E ) = 0

(30)

in which Ai is the Airy function [21] of the first kind. For a given root En , the corresponding normalized eigenfunction is Ai(x − En )

ψn (x) =  ∞ 0

Ai(x − En )2 dx

.

(31)

The Airy function Ai(x) at a point x can be evaluated in Maple by the procedure call AiryAi(x), and its i-th real root can be calculated by the procedure call AiryAiZeros(i). Solving the eigenvalue/eigenfunction problem with either q:-IsCalculated() or q:-IsCalculatedByOctave() and calculating the errors as in the previous example, we obtain for the error in eigenvalues δ En < 1.2 × 10−9 for n ≤ 9, and for the eigenfunctions δψn < 5 × 10−6 for n ≤ 9. 5. An example in quantum cosmology: the Wheeler–DeWitt equation In Section 4.2 we have presented the commands of the package and their usage, taking as an example the Schroedinger equation for the half-harmonic oscillator. In what follows we will present a further example: the solution of the Wheeler–DeWitt equation for the Friedmann–Robertson–Walker (FRW) cosmological model with stiff matter and negative cosmological constant [11]. The quantization of a FRW cosmological model for a stiff matter fluid and negative cosmological constant can be carried out from the classical Hamiltonian [11] H =−

p2a 12a

+ 3a + Λa3 +

pT a3

,

(32)

in which Λ is the cosmological constant, a is the scale factor of the Universe, T is the canonical variable associated to the stiff matter fluid, and (pa , pT ) are the canonically conjugate momenta to (a, T ), respectively. The Schutz variational formalism [2,3] yields the Wheeler– DeWitt equation [11]

 −

1 ∂2

  2  1 ∂ 4 Ψ (a, τ ) − 3a + Λ a Ψ (a, τ ) = i 2 2 12 ∂ a a ∂τ

(33)

in which Ψ (a, τ ) is the wave-function of the Universe and τ = −T plays the role of ‘‘time’’ in this formalism. Separation of the variables yields a trivial solution for the temporal part, leading to Ψ (a, τ ) = η(a)e−iE τ . The time-independent equation then reads



  2  ∂2 12 4 − 2 − 36a + 12Λa η(a) = 2 E η(a), ∂a a

(34)

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

which corresponds to letting f (a) = − 36a2 + 12Λa



Hˆ ≡ −

 4

and g (a) =

12 a2

389

in the general equation (1). Moreover, the Hamiltonian operator

  ∂2 − 36a2 + 12Λa4 2 ∂a

(35)

is self-adjoint [5] with respect to the inner product (2) for the weight function m(a) =

1

. (36) a2 In Ref. [11] we have used Λ = −0.1, N = 100 levels of energy and L = 6 as the interval amplitude for our calculations. Let us define the problem and assign it to the variable p, associating the label ‘‘StiffMatter’’ to it. (In what follows, some output has been removed.)

> p := Problem(NumberOfLevels = 100, Length=6, Ffunction = - 36*a^2 + 12/10 *a^4, Gfunction=12/a^2, Variable = a, Weight=1/a^2, Label="StiffMatter"): The solution can be calculated (with the aid of Octave) by the exported procedure IsSolvedByOctave. > p:-IsSolvedByOctave(): Let us show the 20 lowest eigenvalues.6

> Lowest20 := p:-ItsEigenvalues(1..20); # 20 lowest eigenvalues Lowest20 := -380.2201269503364, -342.1147708535757, -305.9147159189577, -271.6318932398069, -239.2791773149366, -208.8705127744706, -180.4210661160512, -153.9474092349063, -129.4677439579943, -107.0021803515411, -86.5730868917408, -68.20553885594896, -51.9279046312318, -37.77263223488593, -25.77733893053901, -15.98638552857081, -8.453286388628996, -3.24473177627632, -.4483871548924153, 0.4083506212601217e-1

The expected value (17) of the scale factor ⟨a(t )⟩ can be obtained by the procedure ExpectedPosition; its argument is the list of coefficients that represents the linear combination of eigenfunctions. For instance, let us consider an equiprobable superposition (i.e., a linear combination with all of its coefficients equal; normalization is automatic) of the 3 lowest eigenfunctions.

> av := p:-ExpectedPosition([1$3, 0$97]):#Issues a procedure. Now we have a procedure that calculates the expected value of the scale factor for any given instant of time.

> av(0); #Expected scale factor at t=0. 4.03457860034719

The uncertainty (18) on the scale factor as a function of time is obtained by the procedure PositionUncertainty as follows.

> un =: p:-PositionUncertainty([1$3, 0$97]):#Issues a procedure. > un(0); #Uncertainty at t=0. 0.244415887109042

The procedure ExpectedPositionAndUncertaintyAreCalculated generates simultaneously the expected value and uncertainty for several instants of time. In what follows, an equiprobable superposition of the 3 lowest eigenvalues is considered; the expected scale factor and uncertainty are calculated for t ∈ [0, 100], for 100 instants of time, and the results are stored in the file StiffMatter_PositionUncertainty.

> p:-ExpectedPositionAndUncertaintyAreCalculated(0..10,[1$3,0$97],100); "Expected position and uncertainty data stored in file:", "StiffMatter_PositionUncertainty"

The procedure ExpectedPositionAndUncertaintyArePlotted reads the data file and plots the expected value ⟨a(t )⟩ curve (the central curve, shown in red in the online version) as well as the curves ⟨a(t )⟩ ± σ (t ) (the upper and lower curves, shown in green in the online version).

> p:-ExpectedPositionAndUncertaintyArePlotted();

0

6 The spectrum shown here is more accurate than that presented in [11]. In the latter work, the default 10-digit precision of Maple has been inadvertently used in combination with the 15-digit precision of Octave. In the present work, we have used 15-digit precision throughout.

390

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

The solutions can be checked according to the criteria established in Section 3 by the command CheckSolutions:-Numerically, as in

> p:-CheckSolutions:-Numerically(); "Reporting errors in equations in file: ",StiffMatter_Check_Equations "Reporting errors in eigenvector normality in file: ", StiffMatter_Check_Normality "Reporting errors in eigenvectors orthogonality in file: ", StiffMatter_Check_Orthogonality "Done."

The procedure CheckSolutions:-Graphically, in its turn, is able to read data from the aforementioned report files, the root name of which should be passed on to the procedure as an argument, and convey them in graphical form.

> p:-CheckSolutions:-Graphically(); We highlight some traits of the output of such a graphical form in what follows. The first group of data relates to the error obtained when substituting eigenvalues and eigenvectors back into the matricial eigenvalue equation, in the same precision employed in the calculations. The errors are displayed in logarithmic scale. (‘‘Zero’’ errors obtained within the number of digits used in calculations are represented by green dots.) Checking the energy spectrum and eigenvectors using , 15, Digits. Maximum absolute error obtained after substitution of (v,E) in the eigenvector equation M(E)v=0 (Green dots indicate zero deviation.) Reading data from file , StiffMatter_Check_Equations

0

In the example shown, the error ranges from 10−14 to 10−7 , approximately. The next check is related to the normality of eigenvectors, the weight function having been taken into account. Checking absolute deviation from normality of eigenvectors (dot product). (Green dots indicate zero deviation.) Reading data from file , StiffMatter_Check_Normality

0

In this example, normality deviation ranges from 10−15 to 10−11 , approximately. For two specific energy levels, the deviation is inferior to 10−15 , which is rounded to zero (green dot). The last check is that of the orthogonality of eigenvectors. The i-th eigenvector (i ∈ {1, 2, 3, . . .}) is compared to the j-th eigenvectors (j ∈ {i + 1, i + 2, i + 3, . . .}). Here, we show the comparison of the first lowest energy eigenvector to the remaining ones. Checking absolute deviation from orthogonality of eigenvectors (dot product). (Green dots indicate zero deviation.) Reading data from file , StiffMatter_Check_Orthogonality Comparing eigenstate at level , 1, to upper levels.

E.V. Corrêa Silva et al. / Computer Physics Communications 185 (2014) 380–391

391

The deviation from orthogonality ranges from 10−15 to 10−11 , approximately. Acknowledgments E.V. Corrêa Silva (Researcher of CNPq, Brazil) thanks CNPq for partial financial support. G. Oliveira-Neto thanks FAPEMIG for partial financial support. E.V. Corrêa Silva and G.A. Monerat thank UERJ for the Prociência grant, via FAPERJ. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

B.S. DeWitt, Phys. Rev. 160 (1967) 1113. B.F. Schutz, Phys. Rev. D 2 (1970) 2762. B.F. Schutz, Phys. Rev. D 4 (1971) 3559. V.G. Lapchinskii, V.A. Rubakov, Theoret. Math. Phys. 33 (1977) 1076. N.A. Lemos, J. Math. Phys. 37 (1996) 1449. M.J. Gotay, J. Demaret, Phys. Rev. D 28 (1983) 2402. F.G. Alvarenga, J.C. Fabris, N.A. Lemos, G.A. Monerat, Gen. Relativity Gravitation 34 (2002) 651. P. Pedram, et al., Gen. Relativity Gravitation 40 (2008) 1663. P. Pedram, M. Mirzaei, S.S. Gousheh, Comput. Phys. Commun. 176 (2007) 581. E.M. Barboza Jr., N.A. Lemos, Gen. Relativity Gravitation 38 (2006) 1609. G. Oliveira Neto, G.A. Monerat, E.V. Corrêa Silva, C. Neves, L.G. Ferreira Filho, Int. J. Mod. Phys.: Conference Series 3 (2011) 254. G.A. Monerat, G. Oliveira-Neto, E.V. Corrêa Silva, L.G. Ferreira-Filho, P. Romildo Jr., J.C. Fabris, R. Fracalossi, S.V.B. Gonçalves, F.G. Alvarenga, Phys. Rev. D 76 (2007) 024017. G.A. Monerat, L.G. Ferreira-Filho, G. Oliveira-Neto, E.V. Corrêa Silva, C. Neves, Phys. Lett. A 374 (2010) 4741. J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover, New York, 2001. GNU Octave, Available at http://www.gnu.org/software/octave/, (accessed 28.08.2013). GNU Octave, as part of the GNU Project, it is free software under the terms of the GNU General Public License (GPL). H. Sagan, Boundary and Eigenvalue Problems in Mathematical Physics, Dover, New York, 1989. P.A.M. Dirac, The Principles of Quantum Mechanics, fourth ed., Clarendon Press, Oxford, 1967. Maplesoft. Maple 13 Installation and Licensing Guide: Release 13.02. Available at http://www.maplesoft.com/support/install/maple13_install.html (accessed 28.08.2013) Maple is a registered trademark of Waterloo Maple Inc. M.B. Monagan, et al. Maple 13 Advanced Programming Guide. Waterloo, Canada: Maplesoft, 2009. Available at http://www.maplesoft.com/products/maple/history/ documentation.aspx (accessed 28.08.2013). Maple online help. http://www.maplesoft.com/support/help/Maple/view.aspx?path=kernelconnectlost, (accessed 28.08.2013). E. Merzbacher, Quantum Mechanics, third ed., John Willey & Sons, New York, 1998.