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.
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-devel+unsubscr...@googlegroups.com.
To post to this group, send email to sage-devel@googlegroups.com.
Visit this group at https://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.
