Special case for gen_legendre_P with n == m gives incorrect results for x in (-1,1)

[60879728-ecca-41e5-88fc-bb29a80f6c41]

The current implementation of gen_legendre_P is not accurate. No distinction between Ferrers functions and Legendre functions has been made and some special cases are not correctly implemented. This heavily affects spherical harmonics. For a quick fix, to make spherical harmonics work, I propose to temporarily restrict gen_legendre_P to Ferrers functions. In a follow-up #31637, we will make the distinction complete.

---

Old Description

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 functions / Ferrers functions in SageMath 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

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_assoc_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)

This discrepancy also seems to be present in spherical_harmonic when n == m (an instance of SphericalHarmonic from sage/functions/special.py), which is built using gen_legendre_P.

After discussion in the sage-devel mailing list it appears that this is because the n == m case in _eval_special_values_ is based on https://dlmf.nist.gov/14.7#E15, but this is not defined in (-infinity,1].

For the spherical harmonics, where the argument x = cos(theta), x will always be in the range [-1, 1 ], where special case used in _eval_special_values_ is not defined.

On the sage-devel mailing list, Howard Cohl suggested that the correct formula for x in [-1, 1] is

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

See: https://groups.google.com/d/msg/sage-devel/IDtiGF6HB28/QWwnAeLJBAAJ

According to Howard Cohl this is formally a Ferrers function (defined on (-1,1) ), rather than an associated Legendre polynomial. However, the existing code for Func_assoc_legendre_P does not seem to make any distinction between Ferrers and associated Legendre functions.

My proposed fix would be to have Func_assoc_legendre_P._eval_special_values_ choose between two n == m special cases, based on whether -1 <= x <= 1 (above expression) or > 1 (current expression).

This raises the question of whether Func_assoc_legendre_P is correctly defined, as at present it would seem to cover both Ferrers functions and associated Legendre functions.

In my experience with the physics/chemistry literature, the spherical harmonics are universally defined in terms of "associated Legendre functions", even though the argument is x = cos(theta). DLMF suggests these are defined in terms of Ferrers functions of the first kind (https://dlmf.nist.gov/14.30.E1). Wolfram Mathematica does not seem to distinguish. Possibly it is worth flagging in the docstring for Func_assoc_legendre_P that the class seems to cover both functions.

CC: @rwst @slel @fredrik-johansson @tscrim

Component: misc

Keywords: legendre, special function, spherical harmonic

Author: Michael Jung

Branch: 0b14d02

Reviewer: Eric Gourgoulhon, Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/25034

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

c746088

Trac #25034: add abs tol

[mjungmath]

comment:46

I hope this will make it in Sage 9.3.

[egourgoulhon]

comment:47

Thanks again! It's great that Sage has decent spherical harmonics at last!

Travis, do you agree with the positive review?

[tscrim]

comment:49

You have renamed the public function eval_poly, so this need a deprecated_function_alias.

As per comment:39, I consider this to be a regression because you removed _eval_special_values_ from _evalf_:

sage: gen_legendre_P(x,x,.2)                                                                              
gen_legendre_P(x, x, 0.200000000000000)

versus before:

(-0.960000000000000)^(1/2*x)*factorial(2*x)/(2^x*factorial(x))

[mjungmath]

comment:50

Replying to @tscrim:

You have renamed the public function eval_poly, so this need a deprecated_function_alias.

Right!

As per comment:39, I consider this to be a regression because you removed _eval_special_values_ from _evalf_:

sage: gen_legendre_P(x,x,.2)                                                                              
gen_legendre_P(x, x, 0.200000000000000)

versus before:

(-0.960000000000000)^(1/2*x)*factorial(2*x)/(2^x*factorial(x))

This has nothing to do with removing _eval_special_values_. In fact:

sage: assume(x, 'integer')                                                      
sage: gen_legendre_P(x, x, .2)                                                  
0.960000000000000^(1/2*x)*(-1)^x*factorial(2*x)/(2^x*factorial(x))

It is the line

+        if m.is_integer() and n.is_integer():
             ...
             if m == n:

that makes sure that this formula is only used when x is an integer.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c746088 to 82a447c

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

82a447c

Trac #25034: eval_poly deprecation + doc improved

[mjungmath]

comment:52

Here we go. Moreover, I have added some examples in the visible doc that show some floating number examples like the one above.

[tscrim]

comment:54

Okay, thanks for the explanation. This will be good. Just a few more minor things:

• Case|m| > |n|for integers: needs a ::
• Add a period at the end of the .. MATH:: in eval_gen_poly.
• Some formatting: Evaluate the Ferrers function `P(n, m, x)` for `m` and `n` being concrete
• I would move the paragraph These expressions ... `where |m| \leq |n|`. to the main part of the documentation so it is more visible.

Only the first is a must-do, but the rest are good while-we-are-at-it things.

Once these are done, it will be a positive review from me.

[DaveWitteMorris]

comment:55

I think references should go in the master bibliography file (src/doc/en/reference/references/index.rst), instead of the docstring. See #27497.

[mjungmath]

comment:56

Replying to @tscrim:

• I would move the paragraph These expressions ... `where |m| \leq |n|`. to the main part of the documentation so it is more visible.

I agree. However, the documentation of the whole file needs to be refactored imo (see #31636). Thus, I'll leave it there for now.

As for the rest: thank you for taking such a thorough look! :)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

0cda09b

Trac #25034: formatting + typos

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 82a447c to 0cda09b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 0cda09b to 0b14d02

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

0b14d02

Trac #25034: move reference to master bib

[mjungmath]

comment:59

Replying to @DaveWitteMorris:

I think references should go in the master bibliography file (src/doc/en/reference/references/index.rst), instead of the docstring. See #27497.

Is that alright?

Branch pushed to git repo; I updated commit sha1. New commits:

0b14d02

Trac #25034: move reference to master bib

[mjungmath]

comment:60

Patchbot is green. :)

[mkoeppe]

comment:61

Looks like all sufficient conditions for positive review are satisfied.

[vbraun]

Changed branch from public/legendre_25034 to 0b14d02

[mjungmath]

comment:63

Thank you Matthias! :)

[mjungmath]

Changed commit from 0b14d02 to none