Opened 3 years ago
Last modified 3 months ago
#25034 new defect
Special case for gen_legendre_P with n == m gives incorrect results for x in (-1,1)
Reported by: | jcwomack | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-9.3 |
Component: | misc | Keywords: | legendre, special function, spherical harmonic |
Cc: | rws, slelievre, fredrik.johansson | Merged in: | |
Authors: | Reviewers: | ||
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
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.
Change History (8)
comment:1 Changed 3 years ago by
- Description modified (diff)
comment:2 Changed 3 years ago by
- Cc rws slelievre added
- Description modified (diff)
comment:3 Changed 3 years ago by
comment:4 follow-up: ↓ 5 Changed 12 months ago by
For reference, here are some related issues reported by users:
comment:5 in reply to: ↑ 4 ; follow-up: ↓ 6 Changed 12 months ago by
comment:6 in reply to: ↑ 5 Changed 7 months ago by
For reference, here are some related issues reported by users:
One more:
Yet one more:
comment:7 Changed 7 months ago by
- Cc fredrik.johansson added
- Description modified (diff)
- Milestone changed from sage-8.2 to sage-9.2
comment:8 Changed 3 months ago by
- Milestone changed from sage-9.2 to sage-9.3
A (quite severe IMHO) consequence of this bug is
which is plain wrong: the term
sqrt(cos(theta)^2-1)
(which is imaginary fortheta
real) should actually besin(theta)
.