SciELO - Scientific Electronic Library Online

Home Pagealphabetic serial listing  

Services on Demand




Related links

  • On index processCited by Google
  • Have no similar articlesSimilars in SciELO
  • On index processSimilars in Google


Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales

Print version ISSN 0370-3908

Rev. acad. colomb. cienc. exact. fis. nat. vol.37 no.145 Bogotá Oct./Dec. 2013









I. A. Selezneva, Yu. L. Ratis*, E. Hernández**, J. Pérez-Quiles and P. Fernández de Córdoba***

* Samara State Aerospace University named after Academician S.P. Korolev Samara, Russia
** Department of Mathematics, Faculty of Informatics, Universidad de Pinar del Río, Pinar del Río, Cuba
*** Instituto Universitario de Matemática Pura y Aplicada, Universitat Politècnica de València, València, Spain


In this paper we present, with a pedagogical aim, a method to calculate the associated Legendre functions and polynomials. The method uses stable recurrence relations involving these functions.

Keywords: legendre associated functions, Legendre polynomials, recurrence relations, stability.


En este trabajo se presenta, con una finalidad pedagógica, un método para evaluar funciones y polinomios de Legendre. El método utiliza relaciones de recurrencia estables que involucran a estas funciones.

Palabras clave: funciones asociadas de Legendre, polinomios de Legendre, relaciones de recurrencia y estabilidad.



Spherical harmonic functions and associated functions of Legendre of the first and second kind [1,2] are related to most of the special functions of mathematical physics. For example, they appear in the computation of atomic electron configurations [4], in the representation of electromagnetic fields in materials and at surfaces and intefaces [5], and in the characterization of some forms of electromagnetic radiation [6]. Apart from these works, we refer to the reader to [11,12,13] and references therein for other applications and the recent discoveries about Legendre's functions. They are solutions of Legendre's equat.

Solution of this equation are analytical in z, n and m. Usually n is called the degree and m the order. Solutions of this equation are polynomials for n integer and even, and m integer such that 0 ≤ m n. In the case of n odd they are not polynomials, although the whole set of solutions are known as Legendre Polynomials. If m = 0, then the solution are the well-known Legendre polynomials. Usually, when one uses these functions in applied problems, it is necessary to compute them for all n, m in a wide range, which is computationally expensive.

In the case of Legendre's Polynomials, there exists an explicit form of the polynomials Pn(z), given by the formula of Rodrigues,

Expanding this formula, one obtains,

where [n/2] is the integer part of the number n/2.

Whereas for small values of n, the computation of Legendre's polynomials using (3) is trivial, when N » 1 several accuracy and stability problems arise. One should be tempted to use the recurrence relation,

since P0(z) = 1 and P1(z) = z. However, the ascending recursion formula (4) is also unstable. This is related to the fact that in the interval (-1,1) the norm of the polynomial Pn(z) diminishes as 1/√(n) [6]. In other words, the use of a recursion that produces a monotonically descending sequence in norm, leads to a rapid accumulation of round-off errors.

In this paper we present a suitable algorithm to compute the associated Legendre functions Pn(m)(z), which includes as especial case the Legendre polynomials Pn(z). In fact, Pn(z) = Pn(m)(z) for m = 0 and n integer, since [1]

The formula (4) is computationally stable for the regular solutions Pn(z) of Legendre’s equation (1) using backward recursion. In this case, as the norm of polynomials Pn(z) and Pn-1 (z) is smaller than the norm of Pn-2 (z), there is not round-error accumulations. As there are not cancellations, the finiteness of the mantissa does not influence in the error of each step in the recursion.

Of course, if backward recursion has to be applied by using (4), one needs initial values PN(z) and PN-1(z). The evaluation of these values is very expensive for great values of N using formula (3). Furthermore, the computational implementation of (3) generates additional problems, related to the round-off error and also to accumulation of error in the subtraction of closed values.

The computation of the associated Legendre functions Pn(m) (z) in physical and/or engineering problems results as frequent as those of Legendre polynomials Pn(z), so it is very useful that the algorithm computes both the polynomials and functions. Following the same strategy of [7,9,10] it is possible to find the functions Pn(m)(z) and the polynomials Pn(z) in an unified algorithm.

Solutions of equation (1) [8] are subordinated to the recursion relations


where initial values are given by


However, recursion (7) is computationally unstable backwards. The nature of this instability is the same as in the forward recursion (3) for Legendre's polynomials Pn(z). This is a consequence of the rapid accumulation of error on the computed value Yn from the standard trinomial recursion relation

The associated function of Legendre Pn(m)(z) is analytic on m, -n m n, and for a fixed value of n, the norm Pn(m)(z) monotonically increases with m. Connection between the associated functions of Legendre for the negative and positive values m it is given by the relationship, [1]

From the formulas (8), (9), (11) it follows that

Using (11), (12) and (13) it is possible to find the initial values P(-n)(z), Pn(1-n)(z), which permits to compute Legendre Pn(m)(z) for all mn on the basis of the steady ascending recursion (7).

As a particular case, this formula is also useful to compute Pon(z) which can be readily used to compute Poo(z) backwards in a very fast and accurate way. Moreover, the error of the algorithm can be estimated by e = Po appr (z) - 1 and e1 = P1 appr (z) - z. In Figure 1a and 1b we plot the errors for initial n = 100 and n = 30. Both errors are in the range 10-14, i.e. close to the round-off error. As we can see from Fig. 2, typical time for a run is in the range of seconds for 10.000 points.

A last question: as we want to present this work with a methodological point of view, we only sketch the main problems of Legendre recurrence formulas. We refer to [11] or [14] for more formal comparisons.

Codes developed

The subroutine flgndr computes accurately the values of Legendre's function for any z ∈ (-1,1) and all m integer in -nm n. The maximum possible n is 120. This limitation is due to the small value of , which cannot be accurately described in double precision.

In this subroutine, z is a column vector containing the points where Legendre's functions are going to be calculated, and n is the degree. Output is a matrix, pnm, of size [length(z), 2n + 1] containing the values of Pn(m) for -n m n.

function pnm = flgndr(z,n)
nz = length(z);
pnm = zeros(nz,n+1);

fac = prod(2:n);
sqz2 = sqrt((1.0-z.*z));
hsqz2 = 0.5*sqz2;
ihsqz2 = z./hsqz2;


    pnm(:,1) = 1.0;
    pnm(:,1) = -0.5*sqz2;
    pnm(:,2) = z;
    pnm(:,3) = sqz2;

pnm(:,1) = (1-2*abs(n-2*floor(n/2)))
pnm(:,2) = -pnm(:,1)*n.*ihsqz2;
for mr=1:2*n-1



Using the previous routine, it is possible to compute Pn(z) and Pn-1(z), which allow us to compute Legendre polynomials. The inputs are the same that in the flgndr function, but the output contains the values for all the orders from 0 until n. This n is limited in double precision to 120 approximately.

function pn=plgndr(z,n)


% code if(n==0)

if (n==0)



pnm = flgndr(z,n);

pnm = flgndr(z,n-1);

for nc=n-1:-1:1

    pn(:,nc)=((2*nc+1)*z.*pn(:,nc+1)- nc+1)*pn(:,nc+2))/nc ;




In this work we have developed and described (with pedagogical aim) an algorithm which computes Legendre's functions and polynomials. Testing shows how this method can compute this function with a maximum error of 1e-14. The main limitation is due to the small value of Pn(-n), which can be by-passed using 16 bytes precision. A generaliza- tion of this method to another orthogonal polynomials and functions are currently in progress.


Authors thank Universitat Politècnica de València funding under project SP 20120909.



[1] Abramowitz M, and Stegun IA, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York 1972.         [ Links ]

[2] Zwillinger D, Handbook of Differential Equations, Academic Press, Boston, 1997.         [ Links ]

[3] Olver FWJ, and Smith JM, Associated Legendre functions on the cut, Journal of computational Physics, 1983; 51, 502-518.         [ Links ]

[4] Slater JC, Quamtum Theory of Atomic Structure, Volume I, Mac-Graw-Hill, New York, 1960.         [ Links ]

[5] Jackson JD, Classical Electrodynamics, John Wiley and Sons, Inc., New York, 1998.         [ Links ]

[6] Aksenov EP, Special functions in the celestial mechanics, GlavnayaRedaktsiyaFiziko-MatematicheskojLiteratury, Nauka, Moskva, in Russian, 1986, 320 pages.         [ Links ]

[7] RatisYuL, and Fernández de Córdoba P, A code to evaluate (high order) Bessel functions based on the continued fractions method, Computer Physics Communications 1993; 76, 381-388.         [ Links ]

[8] Segura J, and Gil A, Evaluation of Associated Legendre functions off the cut and parabolic cylinder functions, Electronic transactions on numerical Analysis 1999; 9, 137-146.         [ Links ]

[9] Bastardo JL, Abrahamm Ibrahim S, Fernández de Córdoba P, Urchueguía JF, and RatisYuL, Evaluation of Fresnel integrals based on the continued fractions method, Applied Mathematics Letters, 2005; 18, 23-28.         [ Links ]

[10] Selezneva IA, A rapid algorithm for computing special functions in mathematical physics, Vital problems of contemporary natural science, Samara, Intercollegiate collection of scientific papers dedicated to Alexander IvanovichFedosov's memory, 2004; 1, 22-26.         [ Links ]

[11] Gil A, Segura J, and Temme N, Numerical Methods for special functions. SIAM, Philadelphia, PA, 2007.         [ Links ]

[12] Gil A, Segura J, and Temme N, Basic methods for computing special functions. Recent advances in computational and applied mathematics, 67-121, Springer, Dordrecht, 2011.         [ Links ]

[13] Olver FWJ, Lozier DW, Boisvert RF, and Clark CW, NIST Handbook of Mathematical Functions, U. S. Dept. Commerce, Washington. Cambridge University Press, New York, 2010.         [ Links ]

[14] Schneider I, Segura J, Gil A, Guan X, and Bartschat K, A new Fortran 90 program to compute regular and irregular associated Legendre functions. Comput. Phys. Comm. 181 (2010), no.12, 2091-2097.         [ Links ]

Recibido: 28 enero de 2013
Aceptado: 13 de noviembre de 2013