Possible bug in gen_legendre_P (associated Legendre polynomials)

[James Womack]

I would appreciate some help with determining whether I have found a bug, or am simply misusing or misunderstanding some code.

I am using the gen_legendre_P function (an instance of Func_assoc_legendre_P from sage/functions/orthogonal_polys.py) to evaluate associated Legendre polynomials in Sage Math 8.1.

There appears to be a discrepancy in the results I obtain, depending on whether I use gen_legendre_P.eval_poly() or directly call gen_legendre_P() in some cases. I think this is because the _eval_ method first tries to call the _eval_special_values_ method, before using eval_poly.

With this input

x = SR.var('x')
print gen_legendre_P.eval_poly(1,1,x)
print gen_legendre_P(1,1,x)
print gen_legendre_P.eval_poly(1,1,0.5)
print gen_legendre_P(1,1,0.5)

I obtain

-sqrt(-x^2 + 1)
sqrt(x^2 - 1)
-0.866025403784439
5.30287619362453e-17 + 0.866025403784439*I

The result from eval_poly agrees with Mathematica, i.e.

LegendreP[1, 1, 0.5]
-0.866025

This was causing me difficulties because I was using gen_legendre_P(l,m,cos(theta)) to construct real spherical harmonics, but finding that for certain l, m combinations, the resulting spherical harmonics were not real. 

Based on the above output, it seems to me that gen_legendre_P.eval_poly(1,1,cos(theta)) will always be real while gen_legendre_P(1,1,cos(theta)) will be complex (unless |cos(theta)| = 1), since cos(theta) is in the interval [-1,1].

Looking at the code for Func_associ_legendre_P._eval_special_values_, I suspect the culprit is the n == m case, which returns

factorial(2*m)/2**m/factorial(m) * (x**2-1)**(m/2)

Interestingly, it looks like this agrees with 14.7.15 at https://dlmf.nist.gov/14.7.E15, but differs from the definition on Wolfram Mathworld (http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html, Eq. 10), where it is (1-x**2)**(m/2).

I am not an expert in orthogonal polynomials, so perhaps I am missing a subtlety here which means that we can expect this kind of discrepancy? I would be grateful if more experienced mathematicians or Sage Math developers could weigh in on this.

If it is indeed a bug, I am happy to file a report.

[Ralf Stephan]

arb agrees here:

sage: CBF(1/2).legendre_P(1,1)
[-0.8660254037844386 +/- 5.90e-17]

So I'd suggest using complex balls for your numerics until the bug is fixed.

Thanks,

P.S. Still someone should contact DLMF with the right arguments.

[Samuel Lelievre]

Ralf wrote:
> Thanks,
> P.S. Still someone should contact DLMF with the right arguments.

I just emailed them with cc to sage-devel.

[James Womack]

Thanks. I am waiting for an account on Sage trac, then I will submit a bug report.

[Howard Cohl]

There's nothing wrong with the formula. The Legendre function in the DLMF is for arguments greater than 1, and is not valid for arguments less than one. For arguments less than one the correct formula is

P_m^m(x)=(-1)^m (2m)!/(2^m m!) (1-x^2)^(m/2).

Both of these are easy to derive using the well-known formulae for P_\nu^{-\nu} and {\sf P}_\nu^{-\nu} and the connection formulas which relate P_{\nu}^{-m} to P_{\nu}^m, and for Ferrers functions. See http://dlmf.nist.gov/14.5.iv and https://dlmf.nist.gov/14.9.
Where P is the associated Legendre function of the first kind, and {\sf P} is the Ferrers function of the first kind.

[James Womack]

Thanks. If that is the case, then presumably this is a bug in Sage Math and Func_assoc_legendre_P should distinguish the special cases for n == m when x > 1 or x < 1 when evaluating associated Legendre polynomials.

Would you be able to clarify the distinction between Ferrers functions of the first kind and associated Legendre functions for a non-expert? Wolfram Mathworld seems to suggest that they are the same: http://mathworld.wolfram.com/FerrersFunction.html

[Howard Cohl]

The Ferrers functions are defined on the real segment (-1,1).
The associated Legendre functions are in general defined on the Complex plane except for the ray (-\infty,1].

Typically Ferrers functions are written with argument x=\cos\theta, |x|<1 and associated Legendre functions are written with argument z=\cosh\eta, z>1, or something to that effect. One can obtain the Ferrers functions through a limiting process from the associated Legendre functions by taking an appropriately weighted average of x+i0, and x-i0, ,where x\in(-1,1).

If you read the DLMF chapter on associated Legendre functions, then all this is explained quite clearly. The DLMF uses different symbols for the associated Legendre functions vs. the Ferrers functions. An italic P and Q for the associated Legendre functions and a sans serif P and Q for the Ferrers functions. Their definitions and properties are all given in that DLMF chapter. In Abramowitz and Stegun, the same symbol is used for both functions, but you can tell they are using which by the argument z (for associated Legendre) vs. x (for Ferrers functions).

If you have any other questions, please free to email me hcohl -at- nist dot gov.

[Howard Cohl]

Oh, by the way, Wolfram Mathworld is just completely wrong on this page you referenced.
There is a huge difference between the two functions.
Also, there is no such thing as an associated Legendre polynomial.
There is a Legendre polynomials, but if you take the degrees and orders to be integers for associated Legendre functions, they don't always end up being polynomials, as you can see from the formulas which I showed earlier in this thread.

On Thursday, March 22, 2018 at 8:49:16 AM UTC-7, James Womack wrote:

[James Womack]

I have created a ticket on Sage trac for this issue: https://trac.sagemath.org/ticket/25034

As I mention in the ticket, I think that this issue raises a question as to whether the Func_assoc_legendre_P class is correctly defined, given that it now seems to cover both the Ferrers and associated Legendre functions. Given that the intervals on which the functions are defined do not overlap, maybe it makes sense to have a single class for both, but to give it amore generic name?

[Ralf Stephan]

On Tuesday, March 27, 2018 at 11:20:02 AM UTC+2, James Womack wrote:

I have created a ticket on Sage trac for this issue: https://trac.sagemath.org/ticket/25034

Thanks.

 

As I mention in the ticket, I think that this issue raises a question as to whether the Func_assoc_legendre_P class is correctly defined, given that it now seems to cover both the Ferrers and associated Legendre functions. Given that the intervals on which the functions are defined do not overlap, maybe it makes sense to have a single class for both, but to give it amore generic name?

I think it will suffice for now to put the fact in the documentation.

[Eric Gourgoulhon]

Hi,

Le mardi 27 mars 2018 14:46:52 UTC+2, Ralf Stephan a écrit :

I think it will suffice for now to put the fact in the documentation.

I am afraid this is not sufficient: a consequence of this bug is that Sage gives a silly answer for something as elementary as the spherical harmonic Y_1^1(theta, phi), see 
https://trac.sagemath.org/ticket/25034#comment:3

Best wishes,

Eric.

[Michael Jung]

Unfortunately, I am not familiar with the details either.

Nevertheless, I have made the proposed change for the interval -1<x<1 and for now left the name of associated Legendre polynomials even if that might not be rigorously correct. But at least this terminology is in line with https://mathworld.wolfram.com/AssociatedLegendrePolynomial.html. Most importantly, this also solves the problem with spherical harmonics, which was the a highly requested fix.

As for resolving the convention conflicts in a rigorous manner I propose a follow-up ticket: https://trac.sagemath.org/ticket/31637.

Best,
Michael