Exponentially fitted symplectic integrators of RKN type for solving oscillatory problems

A class of explicit modified Runge-Kutta-Nyström (RKN) methods for the numerical integration of second-order IVPs with oscillatory solutions is presented. The symplecticness conditions and the exponential fitting conditions for this class of methods are derived. Based on this conditions, explicit modified RKN integrators with two and three stages per step which have algebraic orders two and four, respectively, are constructed. These new integrators preserve symplecticness when they are applied to Hamiltonian problems, and they integrate exactly differential systems whose solutions can be expressed as linear combinations of the set of functions {exp(λt),exp(−λt)}, λ∈C, or equivalently {sin(ωt),cos(ωt)} when λ=iω, ω∈R. We also analyze the stability properties of the new integrators, obtaining generalized periodicity regions for the classical second-order linear test model. The numerical experiments carried out show that the new methods are more efficient than other symplectic and exponentially fitted codes proposed in the scientific literature.     

Exact analytical expressions and numerical analysis of two-center Franck-Condon factors and matrix elements over displaced harmonic oscillator wave functions

A detailed study is undertaken, using various techniques, in deriving analytical formula of Franck-Condon overlap integrals and matrix elements of various functions of power (xl), exponential (exp(−2cx)) and Gaussian (exp(−cx²)) over displaced harmonic oscillator wave functions with arbitrary frequencies. The results suggested by previous experience with various algorithms are presented in mathematically compact form and consist of generalization. The relationships obtained are valid for the arbitrary values of parameters and the computation results are in good agreement with the literature. The numerical results illustrate clearly a further reduction in calculation times.

Program summary

Program name:FRANCKCatalogue identifier:ADXX_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXX_v1_0Program obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandProgramming language:Mathematica 5.0Computer:Pentium M 1.4 GHzOperating system:Mathematica 5.0RAM:512 MBNo. of lines in distributed program, including test data, etc.:825No. of bytes in distributed program, including test data, etc.:16 344Distribution format:tar.gzNature of problem:The programs calculate the Franck-Condon factors and matrix elements over displaced harmonic oscillator wave functions with arbitrary quantum numbers (n,n₁), frequencies (a,a₁) and displacement (d) for the various functions of power (xl), exponential (exp(−2cx)) and Gaussian (exp(−cx²)).Solution method:The Franck-Condon factors and matrix elements are evaluated using binomial coefficients and basic integrals.Restrictions:The results obtained by the present programs show great numerical stability for arbitrary quantum numbers (n,n₁), frequencies (a,a₁) and displacement (d).Unusual features:NoneRunning time:As an example, for the value of Franck-Condon Overlap Integral Inn^′(d;α,α^′)=0.004405001887372332 with n=3, n₁=2, a=4, a₁=3, d=2, the compilation time in a Pentium M 1.4 GHz computer is 0.18 s. Execution time depends on the values of integral parameters n, n^′, d, α, α^′.     

Optimized implementations of rational approximations—a case study on the Voigt and complex error function

Rational functions are frequently used as efficient yet accurate numerical approximations for real and complex valued special functions. For the complex error function w(x+iy), whose real part is the Voigt function K(x,y), the rational approximation developed by Hui, Armstrong, and Wray [Rapid computation of the Voigt and complex error functions, J. Quant. Spectrosc. Radiat. Transfer 19 (1978) 509-516] is investigated. Various optimizations for the algorithm are discussed. In many applications, where these functions have to be calculated for a large x grid with constant y, an implementation using real arithmetic and factorization of invariant terms is especially efficient.     

Arbitrarily precise numerical solutions of the one-dimensional Schrödinger equation

In this paper, how to overcome the barrier for a finite difference method to obtain the numerical solutions of a one-dimensional Schrödinger equation defined on the infinite integration interval accurate than the computer precision is discussed. Five numerical examples of solutions with the error less than 10⁻⁵⁰ and 10⁻³⁰ for the bound and resonant state, respectively, obtained by the Obrechkoff one-step method implemented in the multi precision mode, which include the harmonic oscillator, the Pöschl-Teller potential, the Morse potential and the Woods-Saxon potential, demonstrate that the finite difference method can yield the eigenvalues of a complex potential with an arbitrarily desired precision within a reasonable efficiency.     

A Runge-Kutta-Nyström pair for the numerical integration of perturbed oscillators

New Runge-Kutta-Nyström methods especially designed for the numerical integration of perturbed oscillators are presented in this paper. They are capable of exactly integrating the harmonic or unperturbed oscillator. We construct an embedded 4(3) RKN pair that is based on the FSAL property. The new method is much more efficient than previously derived RKN methods for some subclasses of problems.     

A Mathematica program for the two-step twelfth-order method with multi-derivative for the numerical solution of a one-dimensional Schrödinger equation

In this paper, we present the detailed Mathematica symbolic derivation and the program which is used to integrate a one-dimensional Schrödinger equation by a new two-step numerical method. We add the fourth- and sixth-order derivatives to raise the precision of the traditional Numerov's method from fourth order to twelfth order, and to expand the interval of periodicity from (0,6) to the one of (0,9.7954) and (9.94792,55.6062). In the program we use an efficient algorithm to calculate the first-order derivative and avoid unnecessarily repeated calculation resulting from the multi-derivatives. We use the well-known Woods-Saxon's potential to test our method. The numerical test shows that the new method is not only superior to the previous lower order ones in accuracy, but also in the efficiency. This program is specially applied to the problem where a high accuracy or a larger step size is required.

Program summary

Title of program: ShdEq.nbCatalogue number: ADTTProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADTTProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputer for which the program is designed and others on which it has been tested: The program has been designed for the microcomputer and been tested on the microcomputer.Computers: IBM PCOperating systems under which the program has been tested: Windows XPProgramming language used: Mathematica 4.2Memory required to execute with typical data: 51 712 bytesNo. of bytes in distributed program, including test data, etc.: 45 381No. of lines in distributed program, including test data, etc.: 7311Distribution format: tar gzip fileCPC Program Library subprograms used: noNature of physical problem: Numerical integration of one-dimensional or radial Schrödinger equation to find the eigenvalues for a bound states and phase shift for a continuum state.Method of solution: Using a two-step method twelfth-order method to integrate a Schrödinger equation numerically from both two ends and the connecting conditions at the matching point, an eigenvalue for a bound state or a resonant state with a given phase shift can be found.Restrictions on the complexity of the problem: The analytic form of the potential function and its high-order derivatives must be known.Typical running time: Less than one second.Unusual features of the program: Take advantage of the high-order derivatives of the potential function and efficient algorithm, the program can provide all the numerical solution of a given Schrödinger equation, either a bound or a resonant state, with a very high precision and within a very short CPU time. The program can apply to a very broad range of problems because the method has a very large interval of periodicity.References: [1] T.E. Simos, Proc. Roy. Soc. London A 441 (1993) 283.[2] Z. Wang, Y. Dai, An eighth-order two-step formula for the numerical integration of the one-dimensional Schrödinger equation, Numer. Math. J. Chinese Univ. 12 (2003) 146.[3] Z. Wang, Y. Dai, An twelfth-order four-step formula for the numerical integration of the one-dimensional Schrödinger equation, Internat. J. Modern Phys. C 14 (2003) 1087.     

Trigonometrically-fitted method with the Fourier frequency spectrum for undamped Duffing equation

With non-linearities, the frequency spectrum of an undamped Duffing oscillator should be composed of odd multiples of the driving frequency which can be interpreted as resonance driving terms. It is expected that the frequency spectrum of the corresponding numerical solution with high accurateness should contain nearly the same components. Hence, to contain these Fourier components and to calculate the amplitudes of these components in a more accurate and efficient way is the key to develop a new numerical method with high stability, accuracy and efficiency for the Duffing equation. To explore the possibility of using trigonometrically-fitting technique to build a numerical method with resonance spectrum, we design four types of Numerov methods, in which the first one is the traditional Numerov method, which contains no Fourier component, the second one contains only the first resonance term, the third one contains the first two resonance terms, and the last one contains the first three resonance terms, and apply them to the well-known undamped Duffing equation with Dooren's parameters. The numerical results demonstrate that the Numerov method fitted with the Fourier components is much more stable, accurate and efficient than the one with no Fourier component. The accuracy of the fitted method with the first three Fourier components can attain 10⁻⁹ for a remarkable range of step sizes, including nearly infinite, except individual small range of instability, which is much higher than the one of the traditional Numerov method, with eight orders for step size of π/2.011.     

A Mathematica program for the approximate analytical solution to a nonlinear undamped Duffing equation by a new approximate approach

According to Mickens [R.E. Mickens, Comments on a Generalized Galerkin's method for non-linear oscillators, J. Sound Vib. 118 (1987) 563], the general HB (harmonic balance) method is an approximation to the convergent Fourier series representation of the periodic solution of a nonlinear oscillator and not an approximation to an expansion in terms of a small parameter. Consequently, for a nonlinear undamped Duffing equation with a driving force Bcos(ωx), to find a periodic solution when the fundamental frequency is identical to ω, the corresponding Fourier series can be written as

     

A trigonometrically-fitted one-step method with multi-derivative for the numerical solution to the one-dimensional Schrödinger equation

In this paper we present a new multi-derivative or Obrechkoff one-step method for the numerical solution to an one-dimensional Schrödinger equation. By using trigonometrically-fitting method (TFM), we overcome the traditional Obrechkoff one-step method (or called as the non-TFM) for its poor-accuracy in the resonant state. In order to demonstrate the excellent performance for the resonant state, we consider only the simplest TFM, of which the local truncation error (LTE) is of O(h⁷), a little higher than the one of the traditional Numerov method of O(h⁶), and only the first- and second-order derivatives of the potential function are needed. In the new method, in order to solve two unknowns, wave function and its first-order derivative, we use a pair of two symmetrically linear-independent one-step difference equations. By applying it to the well-known Woods-Saxon's potential problem, we find that the TFM can surpass the non-TFM by five orders for the highest resonant state, and surpass Numerov method by eight orders. On the other hand, because of the small error constant, the accuracy improvement to the ground state is also remarkable, and the numerical result obtained by TFM can be four to five orders higher than the one by Numerov method.     

Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Pöschl-Teller-Ginocchio potential wave functions

The fast computation of the Gauss hypergeometric function with all its parameters complex is a difficult task. Although the function verifies numerous analytical properties involving power series expansions whose implementation is apparently immediate, their use is thwarted by instabilities induced by cancellations between very large terms. Furthermore, small areas of the complex plane, in the vicinity of , are inaccessible using power series linear transformations. In order to solve these problems, a generalization of R.C. Forrey's transformation theory has been developed. The latter has been successful in treating the function with real parameters. As in real case transformation theory, the large canceling terms occurring in analytical formulas are rigorously dealt with, but by way of a new method, directly applicable to the complex plane. Taylor series expansions are employed to enter complex areas outside the domain of validity of power series analytical formulas. The proposed algorithm, however, becomes unstable in general when |a|, |b|, |c| are moderate or large. As a physical application, the calculation of the wave functions of the analytical Pöschl-Teller-Ginocchio potential involving evaluations is considered.

Program summary

Program title: hyp_2F1, PTG_wfCatalogue identifier: AEAE_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEAE_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 6839No. of bytes in distributed program, including test data, etc.: 63 334Distribution format: tar.gzProgramming language: C++, Fortran 90Computer: Intel i686Operating system: Linux, WindowsWord size: 64 bitsClassification: 4.7Nature of problem: The Gauss hypergeometric function , with all its parameters complex, is uniquely calculated in the frame of transformation theory with power series summations, thus providing a very fast algorithm. The evaluation of the wave functions of the analytical Pöschl-Teller-Ginocchio potential is treated as a physical application.Solution method: The Gauss hypergeometric function verifies linear transformation formulas allowing consideration of arguments of a small modulus which then can be handled by a power series. They, however, give rise to indeterminate or numerically unstable cases, when b−a and c−a−b are equal or close to integers. They are properly dealt with through analytical manipulations of the Lanczos expression providing the Gamma function. The remaining zones of the complex plane uncovered by transformation formulas are dealt with Taylor expansions of the function around complex points where linear transformations can be employed. The Pöschl-Teller-Ginocchio potential wave functions are calculated directly with evaluations.Restrictions: The algorithm provides full numerical precision in almost all cases for |a|, |b|, and |c| of the order of one or smaller, but starts to be less precise or unstable when they increase, especially through a, b, and c imaginary parts. While it is possible to run the code for moderate or large |a|, |b|, and |c| and obtain satisfactory results for some specified values, the code is very likely to be unstable in this regime.Unusual features: Two different codes, one for the hypergeometric function and one for the Pöschl-Teller-Ginocchio potential wave functions, are provided in C++ and Fortran 90 versions.Running time: 20,000 function evaluations take an average of one second.     

A numerical differentiation library exploiting parallel architectures

We present a software library for numerically estimating first and second order partial derivatives of a function by finite differencing. Various truncation schemes are offered resulting in corresponding formulas that are accurate to order O(h), O(h²), and O(h⁴), h being the differencing step. The derivatives are calculated via forward, backward and central differences. Care has been taken that only feasible points are used in the case where bound constraints are imposed on the variables. The Hessian may be approximated either from function or from gradient values. There are three versions of the software: a sequential version, an OpenMP version for shared memory architectures and an MPI version for distributed systems (clusters). The parallel versions exploit the multiprocessing capability offered by computer clusters, as well as modern multi-core systems and due to the independent character of the derivative computation, the speedup scales almost linearly with the number of available processors/cores.

Program summary

Program title: NDL (Numerical Differentiation Library)Catalogue identifier: AEDG_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEDG_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 73 030No. of bytes in distributed program, including test data, etc.: 630 876Distribution format: tar.gzProgramming language: ANSI FORTRAN-77, ANSI C, MPI, OPENMPComputer: Distributed systems (clusters), shared memory systemsOperating system: Linux, SolarisHas the code been vectorised or parallelized?: YesRAM: The library uses O(N) internal storage, N being the dimension of the problemClassification: 4.9, 4.14, 6.5Nature of problem: The numerical estimation of derivatives at several accuracy levels is a common requirement in many computational tasks, such as optimization, solution of nonlinear systems, etc. The parallel implementation that exploits systems with multiple CPUs is very important for large scale and computationally expensive problems.Solution method: Finite differencing is used with carefully chosen step that minimizes the sum of the truncation and round-off errors. The parallel versions employ both OpenMP and MPI libraries.Restrictions: The library uses only double precision arithmetic.Unusual features: The software takes into account bound constraints, in the sense that only feasible points are used to evaluate the derivatives, and given the level of the desired accuracy, the proper formula is automatically employed.Running time: Running time depends on the function's complexity. The test run took 15 ms for the serial distribution, 0.6 s for the OpenMP and 4.2 s for the MPI parallel distribution on 2 processors.     

Evolutionary fitting methods for the extraction of mass spectra in lattice field theory

We present an application of evolutionary algorithms to the curve-fitting problems commonly encountered when trying to extract particle masses from correlators in lattice QCD. Harnessing the flexibility of evolutionary methods in global optimization allows us to dynamically adapt the number of states to be fitted along with their energies so as to minimize overall χ²/(d.o.f.), leading to a promising new way of extracting the mass spectrum from measured correlation functions.     

A new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation

In this paper we present a new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation. The scheme uses a symmetrical multi-point difference formula to represent the partial differentials of the two-dimensional variables, which can improve the accuracy of the numerical solutions to the order of Δx2Nq+2 when a (2Nq+1)-point formula is used for any positive integer Nq with Δx=Δy, while Nq=1 equivalent to the traditional scheme. On the other hand, the new scheme keeps the same form of the traditional matrix equation so that the standard algebraic eigenvalue algorithm with a real, symmetric, large sparse matrix is still applicable. Therefore, for the same dimension, only a little more CPU time than the traditional one should be used for diagonalizing the matrix. The numerical examples of the two-dimensional harmonic oscillator and the two-dimensional Henon-Heiles potential demonstrate that by using the new method, the error in the numerical solutions can be reduced steadily and extensively through the increase of Nq, which is more efficient than the traditional methods through the decrease of the step size.     

Numerical solving of the vibrational time-independent Schrödinger equation in one and two dimensions using the variational method

A program package for variational solving of the time-independent Schrödinger equation (SE) in one and two dimensions is described. The first part of the the program package includes the fitting program (FIT) with which the ab initio or DFT calculated points are fitted to a computationally inexpensive functional form. Proper fitting of the potential energy surface is crucial for the quality of the results. The second part of the package consists of a program for variational solving of the SE (2DSCHRODINGER) using either a shifted Gaussian basis set or the rectangular basis set proposed by Balint-Kurti and coworkers [J. Chem. Phys. 91 (1989) 3571]. The third part of the program package consists of the calculation of the expectation values, IR and Raman spectra XPECT), and the visualization of results (PLOT). The program package is applied to study a quantum harmonic oscillator and an intramolecular, strong hydrogen bond in picolinic acid N-oxide. For the former system analytical solutions exist, while for the latter system a comparison with the experimental data is made. The advantages and disadvantages of the applied methods are discussed.     

The Gröbner basis of the ideal of vanishing polynomials

We construct an explicit minimal strong Gröbner basis of the ideal of vanishing polynomials in the polynomial ring over Z/m for m≥2. The proof is done in a purely combinatorial way. It is a remarkable fact that the constructed Gröbner basis is independent of the monomial order and that the set of leading terms of the constructed Gröbner basis is unique, up to multiplication by units. We also present a fast algorithm to compute reduced normal forms, and furthermore, we give a recursive algorithm for building a Gröbner basis in Z/m[x₁,x₂,…,x_n] along the prime factorization of m. The obtained results are not only of mathematical interest but have immediate applications in formal verification of data paths for microelectronic systems-on-chip.     

On the evaluation of H and X functions

We indicate a Double Exponential Formula based numerical integration method for the evaluation of the Ambarzumian-Chandrasekhar H function and the X function of neutron transport for the single speed and isotropic case. This method is significantly more economical than our earlier scheme, which was based on IMT quadrature. For c<5, the present method converges faster than our earlier IMT scheme. This will be adequate for all radiative transport and transport theory applications. These findings are supported by appropriate error analysis. Unlike the IMT method, the DE quadrature nodes are generated by a simple algebraic expression which is a great advantage.     

f1: a code to compute Appell's F1 hypergeometric function

In this work we present the FORTRAN code to compute the hypergeometric function F₁(α,β₁,β₂,γ,x,y) of Appell. The program can compute the F₁ function for real values of the variables {x,y}, and complex values of the parameters {α,β₁,β₂,γ}. The code uses different strategies to calculate the function according to the ideas outlined in [F.D. Colavecchia et al., Comput. Phys. Comm. 138 (1) (2001) 29].

Program summary

Title of the program: f1Catalogue identifier: ADSJProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADSJProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputers: PC compatibles, SGI Origin2∗Operating system under which the program has been tested: Linux, IRIXProgramming language used: Fortran 90Memory required to execute with typical data: 4 kbytesNo. of bits in a word: 32No. of bytes in distributed program, including test data, etc.: 52 325Distribution format: tar gzip fileExternal subprograms used: Numerical Recipes hypgeo [W.H. Press et al., Numerical Recipes in Fortran 77, Cambridge Univ. Press, 1996] or chyp routine of R.C. Forrey [J. Comput. Phys. 137 (1997) 79], rkf45 [L.F. Shampine and H.H. Watts, Rep. SAND76-0585, 1976].Keywords: Numerical methods, special functions, hypergeometric functions, Appell functions, Gauss functionNature of the physical problem: Computing the Appell F₁ function is relevant in atomic collisions and elementary particle physics. It is usually the result of multidimensional integrals involving Coulomb continuum states.Method of solution: The F₁ function has a convergent-series definition for |x|<1 and |y|<1, and several analytic continuations for other regions of the variable space. The code tests the values of the variables and selects one of the precedent cases. In the convergence region the program uses the series definition near the origin of coordinates, and a numerical integration of the third-order differential parametric equation for the F₁ function. Also detects several special cases according to the values of the parameters.Restrictions on the complexity of the problem: The code is restricted to real values of the variables {x,y}. Also, there are some parameter domains that are not covered. These usually imply differences between integer parameters that lead to negative integer arguments of Gamma functions.Typical running time: Depends basically on the variables. The computation of Table 4 of [F.D. Colavecchia et al., Comput. Phys. Comm. 138 (1) (2001) 29] (64 functions) requires approximately 0.33 s in a Athlon 900 MHz processor.     

JADAMILU: a software code for computing selected eigenvalues of large sparse symmetric matrices

A new software code for computing selected eigenvalues and associated eigenvectors of a real symmetric matrix is described. The eigenvalues are either the smallest or those closest to some specified target, which may be in the interior of the spectrum. The underlying algorithm combines the Jacobi-Davidson method with efficient multilevel incomplete LU (ILU) preconditioning. Key features are modest memory requirements and robust convergence to accurate solutions. Parameters needed for incomplete LU preconditioning are automatically computed and may be updated at run time depending on the convergence pattern. The software is easy to use by non-experts and its top level routines are written in FORTRAN 77. Its potentialities are demonstrated on a few applications taken from computational physics.

Program summary

Program title: JADAMILUCatalogue identifier: ADZT_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADZT_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 101 359No. of bytes in distributed program, including test data, etc.: 7 493 144Distribution format: tar.gzProgramming language: Fortran 77Computer: Intel or AMD with g77 and pgf; Intel EM64T or Itanium with ifort; AMD Opteron with g77, pgf and ifort; Power (IBM) with xlf90.Operating system: Linux, AIXRAM: problem dependentWord size: real:8; integer: 4 or 8, according to user's choiceClassification: 4.8Nature of problem: Any physical problem requiring the computation of a few eigenvalues of a symmetric matrix.Solution method: Jacobi-Davidson combined with multilevel ILU preconditioning.Additional comments: We supply binaries rather than source code because JADAMILU uses the following external packages:

•

MC64. This software is copyrighted software and not freely available. COPYRIGHT (c) 1999 Council for the Central Laboratory of the Research Councils.

•

AMD. Copyright (c) 2004-2006 by Timothy A. Davis, Patrick R. Amestoy, and Iain S. Duff. All Rights Reserved. Source code is distributed by the authors under the GNU LGPL licence.

•

BLAS. The reference BLAS is a freely-available software package. It is available from netlib via anonymous ftp and the World Wide Web.

•

LAPACK. The complete LAPACK package or individual routines from LAPACK are freely available on netlib and can be obtained via the World Wide Web or anonymous ftp.

•

For maximal benefit to the community, we added the sources we are proprietary of to the tar.gz file submitted for inclusion in the CPC library. However, as explained in the README file, users willing to compile the code instead of using binaries should first obtain the sources for the external packages mentioned above (email and/or web addresses are provided).

Running time: Problem dependent; the test examples provided with the code only take a few seconds to run; timing results for large scale problems are given in Section 5.     

Quasilinearization approach to computations with singular potentials

We pioneered the application of the quasilinearization method (QLM) to the numerical solution of the Schrödinger equation with singular potentials. The spiked harmonic oscillator r²+λr−α is chosen as the simplest example of such potential. The QLM has been suggested recently for solving the Schrödinger equation after conversion into the nonlinear Riccati form. In the quasilinearization approach the nonlinear differential equation is treated by approximating the nonlinear terms by a sequence of linear expressions. The QLM is iterative but not perturbative and gives stable solutions to nonlinear problems without depending on the existence of a smallness parameter. The choice of zero iteration is based on general features of solutions near the boundaries.We show that the energies of bound state levels in the spiked harmonic oscillator potential which are notoriously difficult to compute for small couplings λ, are easily calculated with the help of QLM for any λ and α with accuracy of twenty significant figures.     

Gauss-Legendre and Chebyshev quadratures for singular integrals

Exact expressions are presented for efficient computation of the weights in Gauss-Legendre and Chebyshev quadratures for selected singular integrands. The singularities may be of Cauchy type, logarithmic type or algebraic-logarithmic end-point branching points. We provide Fortran 90 routines for computing the weights for both the Gauss-Legendre and the Chebyshev (Fejér-1) meshes whose size can be set by the user.

New program summary

Program title: SINGQUADCatalogue identifier: AEBR_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEBR_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 4128No. of bytes in distributed program, including test data, etc.: 25 815Distribution format: tar.gzProgramming language: Fortran 90Computer: Any with a Fortran 90 compilerOperating system: Linux, Windows, MacRAM: Depending on the complexity of the problemClassification: 4.11Nature of problem: Program provides Gauss-Legendre and Chebyshev (Fejér-1) weights for various singular integrands.Solution method: The weights are obtained from the condition that the quadrature of order N must be exact for a polynomial of degree?(N−1). The weights are expressed as moments of the singular kernels associated with Legendre or Chebyshev polynomials. These moments are obtained in analytic form amenable for computation.Additional comments: If the NAGWare f95 compiler is used, the option, “-kind = byte”, must be included in the compile command lines of the Makefile.Running time: The test run supplied with the distribution takes a couple of seconds to execute.