A new Fortran 90 program to compute regular and irregular associated Legendre functions

https://doi.org/10.1016/j.cpc.2010.08.038Get rights and content

Abstract

We present a modern Fortran 90 code to compute the regular Plm(x) and irregular Qlm(x) associated Legendre functions for all x(1,+1) (on the cut) and |x|>1 and integer degree (l) and order (m). The code applies either forward or backward recursion in (l) and (m) in the stable direction, starting with analytically known values for forward recursion and considering both a Wronskian based and a modified Miller's method for backward recursion. While some Fortran 77 codes existed for computing the functions off the cut, no Fortran 90 code was available for accurately computing the functions for all real values of x different from x=±1 where the irregular functions are not defined.

Program summary

Program title: Associated Legendre Functions

Catalogue identifier: AEHE_v1_0

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEHE_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.: 6722

No. of bytes in distributed program, including test data, etc.: 310 210

Distribution format: tar.gz

Programming language: Fortran 90

Computer: Linux systems

Operating system: Linux

RAM: bytes

Classification: 4.7

Nature of problem: Compute the regular and irregular associated Legendre functions for integer values of the degree and order and for all real arguments. The computation of the interaction of two electrons, 1/|r1r2|, in prolate spheroidal coordinates is used as one example where these functions are required for all values of the argument and we are able to easily compare the series expansion in associated Legendre functions and the exact value.

Solution method: The code evaluates the regular and irregular associated Legendre functions using forward recursion when |x|<1 starting the recursion with the analytically known values of the first two members of the sequence. For values of the argument |x|<1, the upward recursion over the degree for the regular functions is numerically stable. For the irregular functions, backward recursion must be applied and a suitable method of starting the recursion is required. The program has two options; a modified version of Miller's algorithm and the use of the Wronskian relation between the regular and irregular functions, which was the method considered in [1]. Both approaches require the computation of a continued fraction to begin the recursion. The Wronskian method (which can also be described as a modified Miller's method) is a convenient method of computations when both the regular and irregular functions are needed.

Running time: The example tests provided take a few seconds to run.

References:

  • [1]

    A. Gil, J. Segura, A code to evaluate prolate and oblate spheroidal harmonics, Comput. Phys. Commun. 108 (1998) 267–278.

Introduction

Regular Plm(x) and irregular Qlm(x) associated Legendre functions of integer l and m are needed for many applications in applied mathematics and computational physics. For example, the regular Legendre functions on the interval [1,+1] (cut) (x=cos(θ)) are used as the fundamental angular functions in spherical harmonics which are ubiquitous in atomic and molecular physics and are also used in numerous applications in applied mathematics as spectral expansion functions to solutions of differential and integral equations. There are also applications where both the regular and irregular functions appear outside the cut. One such example is the Neumann expansion of the interaction of two electrons, 1/|r1r2|, in prolate spheroidal coordinates. The prolate spheroidal coordinate system is a near optimal choice to treat quantum mechanical problems involving diatomic molecules not far from their equilibrium geometries. In this application numerous two-electron integrals involving the Neumann expansion are required to form a matrix representation of the molecular Hamiltonian and the basic two-electron integrals involve the associated Legendre functions off the cut.

In this paper we present a Fortran 90 code to evaluate Plm(x) and Qlm(x) for all real values of the argument x different from ±1 and integer degree and order. The code is based on well-known recursion relations for these functions in degree and order. Computation of both the regular and irregular functions on the cut is easily and stably implemented by forward recursion starting with the first two members of the sequence. As we demonstrate below, the two starting values may be easily obtained analytically. For values of x outside the cut, Plm(x) may also be accurately and stably computed for all degree and order using upward recurrence over the degree l. Unfortunately, this is not the case for Qlm(x). For Qlm(x), the stable direction for recursion in l is downward and the stable direction for recursion in m is upward. In addition, for values close to the cut, one must be very careful about losing accuracy as the irregular functions blow up as |x| approaches 1. The approach taken by us is to find an accurate downward recursion in l for m=(0,1) and to then recur upward in m for all the remaining m values. We have provided two alternatives to perform the downward recursion in l. The first employs a modified version of Miller's algorithm. In the standard version of Miller's algorithm one begins the downward recursion at some high value of Lmax, with the first two starting values of the function being chosen as zero and one. It has been demonstrated that this process will produce values proportional to the correct values of the function provided a high enough value of Lmax is chosen. By repeating the process for a higher value of Lmax, one is able to check the accuracy and stop when that is sufficient for the intended purpose. If the functions are always different from zero for any x, which is the case of Legendre functions off the cut, the proportionality constant required to “renormalize” the functions may be easily found by examining the ratio of the computed to exact value of the lowest member in the sequence, which may be found analytically. This is then used to scale all of the functions to their correct values. The fundamental problem with Millers algorithm is that it is often unclear what starting value is needed to obtain the required accuracy and it is becomes necessary to repeat the process a number of times which can be inefficient.

The modified Miller approach considered in this paper does not require prior knowledge of Lmax. All that is required is to be able to compute the ratio Tl=QlmQl1m for the last l needed. This ratio may be expressed as a continued fraction by simply rewriting the recursion relation in terms of ratios. The continued fraction is computed using the Lentz–Thompson algorithm [6] (see also [1, Section 6.6.2]). By setting Ql1m=e with e a small number, say, 1.e–300, and Qlm=eTl, the downward recursion is initiated with correct values of the functions to within an overall scaling. The renormalization is done precisely as in the standard Miller method. The second approach, which was also considered in [2], is to use the continued fraction for the ratio of the two highest members of the sequence and by combining that with the known Wronskian for the regular and irregular functions, compute two accurate values for the starting values of the sequence. To compute the Wronskian, it is necessary to calculate values of the regular function but these are easily obtained by stable upward recursion or are already known from previous computation. This can also be seen as a modified Miller's method, but with the advantage that no renormalization is needed. When both solutions are needed (regular and irregular), this is a convenient method of computation.

Section snippets

Theoretical background

The appropriate definition of the regular and irregular associated Legendre functions depend on whether the argument x is on or off the cut. To treat both cases on the same footing, we define,sf=sign(1x2), and write the recurrences over l and m as follows,(ml1)Pl+1m(x)+(2l+1)xPlm(x)(l+m)xPl1m(x)=0,Plm+1(x)+2mxsf(1x2)Plm(x)+sf(l+m)(lm+1)Plm1(x)=0. The same recurrences apply to the irregular solutions Qlm(x).

Overview of the software structure

The program is written in Fortran 90 and uses many of the new features of that language such as optional arguments to subroutines, dynamic memory allocation, and derived types. The code consists of two short driver programs Legendre_Function_Driver.f90 and Prolate_Driver.f90 and a number of modules containing shared variables, allocated arrays, often used data and the needed subroutines. The driver Legendre_Function_Driver.f90 is used to compute the associated Legendre functions for the l, m

Description of the individual software components

Here we provide descriptions of the individual subroutines.

Installation instructions

The code may be downloaded as the tar file associated_legendre_functions.tar. After the code is untarred, there will be a makefile, driver routines and some subroutines and Modules. If you run the makefile, the executables will be produced. Running them will produce the output files which can then be compared with those in this article.

Test run description

To produce an executable file, type make F90=your_fortran_compiler. We have compiled the code with a number of compilers and have not encountered any difficulties. The makefile produces both executables discussed above, Legendre_Function_Driver_x and Prolate_Driver_x. There are two input files for testing the drivers and their names are obvious if one looks. There are two options. The first is to simply compute the regular, irregular or both associated Legendre functions for a range of

Acknowledgements

We thank Nico Temme for bringing a few of us (BIS, AG, JS) together on this project. This work was supported, in part by the United States National Science Foundation under grant PHY-0757755 (XG and KB). AG and JS acknowledge support from Ministerio de Ciencia e Innovación, project MTM2009-11686.

References (7)

There are more references available in the full text version of this article.

Cited by (0)

This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).

1

Corresponding authors.

View full text