## Services on Demand

## Journal

## Article

## Indicators

- Cited by SciELO
- Access statistics

## Related links

- Cited by Google
- Similars in SciELO
- Similars in Google

## Share

## 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

**MATEMÁTICAS**

**A CODE TO CALCULATE HIGH ORDER LEGENDRE POLYNOMIALS AND FUNCTIONS**

**UN CÓDIGO PARA EL CÁLCULO DE POLINOMIOS Y FUNCIONES DE LEGENDRE DE ALTO ORDEN**

**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

**ABSTRACT**

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.

**RESUMEN**

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.

**Introduction**

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 *P _{n}*(

*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 *P _{0}*(

*z*) = 1 and

*P*(

_{1}*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

*P*(

_{n}*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 *P _{n}*

^{(m)}(

*z*), which includes as especial case the Legendre polynomials

*P*(

_{n}*z*). In fact,

*P*(

_{n}*z*) =

*P*

_{n}^{(m)}(

*z*) for

*m*= 0 and

*n*integer, since [1]

The formula (4) is computationally stable for the regular solutions *P _{n}*(

*z*) of Legendre’s equation (1) using backward recursion. In this case, as the norm of polynomials

*P*(

_{n}*z*) and

*P*(

_{n-1}*z*) is smaller than the norm of

*P*(

_{n-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 *P _{N}*(

*z*) and

*P*(

_{N-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 *P _{n}*

*(*

^{(m)}*z*) in physical and/or engineering problems results as frequent as those of Legendre polynomials

*P*(

_{n}*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

*P*

_{n}*(*

^{(m)}*z*) and the polynomials

*P*(

_{n}*z*) in an unified algorithm.

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

and

where initial values are given by

and

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 *P _{n}*(

*z*). This is a consequence of the rapid accumulation of error on the computed value

*Y*from the standard trinomial recursion relation

_{n}

The associated function of Legendre *P _{n}^{(m)}*(

*z*) is analytic on

*m*,

*-n*≤

*m*≤

*n*, and for a fixed value of

*n*, the norm

*P*

_{n}*(*

^{(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*),

*P*

_{n}*(*

^{(1-n)}*z*), which permits to compute Legendre

*P*

_{n}*(*

^{(m)}*z*) for all

*m*≤

*n*on the basis of the steady ascending recursion (7).

As a particular case, this formula is also useful to compute *P*^{o}* _{n}*(

*z*) which can be readily used to compute

*P*

^{o}_{o}(

*z*) backwards in a very fast and accurate way. Moreover, the error of the algorithm can be estimated by e =

*P*

_{o}

*(*

^{appr}*z*) - 1 and e

_{1}=

*P*

_{1}

*(*

^{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 *-n* ≤ *m *≤ *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*), 2*n* + 1] containing the values of *P _{n}*

*for*

^{(m) }*-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;

if(n==0)

- pnm(:,1) = 1.0;

return

if(n==1)

- pnm(:,1) = -0.5*sqz2;

pnm(:,2) = z;

pnm(:,3) = sqz2;

return

pnm(:,1) = (1-2*abs(n-2*floor(n/2)))

*hsqz2.^n/fac;

pnm(:,2) = -pnm(:,1)*n.*ihsqz2;

for mr=1:2*n-1

- pnm(:,mr+2)=(mr-n).*ihsqz2.

end

end

Using the previous routine, it is possible to compute *P _{n}*(

*z*) and

*P*

*(*

_{n-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)

pn=zeros(length(z),n+1);

% code if(n==0)

if (n==0)

- pn(:,1)=1;

return

end

pnm = flgndr(z,n);

pn(:,n+1)=pnm(:,n+1);

pnm = flgndr(z,n-1);

pn(:,n)=pnm(:,n);

for nc=n-1:-1:1

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

end

**Conclusion**

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 *P _{n}*

^{(-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.

**Acknowledgements**

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

**Bibliografía**

[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