[Scipy-tickets] [SciPy] #1296: associated legendre polynomials in scipy.special fail for calculations with degree higher than approx. 22
SciPy Trac
scipy-tickets@scipy....
Sat Oct 9 16:18:33 CDT 2010
#1296: associated legendre polynomials in scipy.special fail for calculations with
degree higher than approx. 22
------------------------------+---------------------------------------------
Reporter: josip.mihaljevic | Owner: pv
Type: defect | Status: new
Priority: high | Milestone: 0.9.0
Component: scipy.special | Version: 0.7.0
Keywords: legendre |
------------------------------+---------------------------------------------
Comment(by warren.weckesser):
(It looks like you have 'order' and 'degree' reversed in your description.
The first argument of lpmv is the order, and the second is the degree.
Your title of the bug report matches this usage.)
It does appear that lpmv is not a very good implementation. I have
attached a file that contains two different implementations of the
associated Legendre polynomial. One is based on a recursion formula (see
Numerical Recipes for a similar implementation), and one is based on an
integral representation. Both are much more accurate than lpmv. (The
first two arguments of assoc_legendre_p and assoc_legendre_p2 are the
degree then the order.)
Order 0, degree 16:
{{{
In [325]: lpmv(0, 16, 0.001)
Out[325]: 0.19635390721305157
In [326]: assoc_legendre_p(16, 0, 0.001)
Out[326]: 0.19635390806272032
In [327]: assoc_legendre_p2(16, 0, 0.001)
Out[327]: (0.19635390806272018, 2.3468911764622238e-15)
}}}
Maxima's assoc_legendre_p:
{{{
assoc_legendre_p(16, 0, 0.001);
interval(.1963539080627201,5.223692285805479*10^-15)
}}}
Order 0, degree 22:
{{{
In [328]: lpmv(0, 22, 0.001)
Out[328]: -0.16814844246255234
In [329]: assoc_legendre_p(22, 0, 0.001)
Out[329]: -0.16814554527766912
In [330]: assoc_legendre_p2(22, 0, 0.001)
Out[330]: (-0.16814554527766906, 4.7065059058218157e-14)
}}}
Maxima:
{{{
assoc_legendre_p(22, 0, 0.001);
interval(-.1681455452776688,5.982697311505185*10^-15)
}}}
Order 0, degree 40 (lpmv is garbage):
{{{
In [331]: lpmv(0, 40, 0.001)
Out[331]: -21449.235366821289
In [332]: assoc_legendre_p(40, 0, 0.001)
Out[332]: 0.1252678976534484
In [333]: assoc_legendre_p2(40, 0, 0.001)
Out[333]: (0.12526789765344834, 1.1965288185813496e-13)
}}}
Maxima:
{{{
assoc_legendre_p(40, 0, 0.001);
interval(.1252678976534483,7.666725305397109*10^-15)
}}}
--
Ticket URL: <http://projects.scipy.org/scipy/ticket/1296#comment:1>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
More information about the Scipy-tickets
mailing list