id summary reporter owner description type status priority milestone component resolution keywords cc merged author reviewer upstream work_issues branch commit dependencies stopgaps
25034 Special case for gen_legendre_P with n == m gives incorrect results for x in (-1,1) jcwomack "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 [[https://groups.google.com/d/msg/sage-devel/IDtiGF6HB28/ErLsqI1eBAAJ|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." defect closed blocker sage-9.3 misc fixed legendre, special function, spherical harmonic rws slelievre fredrik.johansson tscrim Michael Jung Eric Gourgoulhon, Travis Scrimshaw N/A 0b14d02cbccf9d46149589acba7ee05dc2f12397