symbolic Legendre / associated Legendre functions / polynomials

[rwst]

Defect because there is no Sage binding for a result returned by Maxima:

sage: hypergeometric([-1/2,1/2],[2],4).simplify_hypergeometric()
1/2*I*sqrt(3)*assoc_legendre_p(1/2, -1, -5/3)
sage: assoc_legendre_p
...
NameError: name 'assoc_legendre_p' is not defined

Existing numeric functions are legendre_P, legendre_Q, gen_legendre_P, and gen_legendre_Q which correspond to P(n,x) / Q(n,x) and associated Legendre P(n,m,x) / Q(n,m,x).

They should be made symbolic. FLINT has fast code for P(n,x).

• PQ(n,x)
  • https://en.wikipedia.org/wiki/Legendre_polynomials
  • http://mathworld.wolfram.com/LegendrePolynomial.html
• PQ(n,m,x)
  • https://en.wikipedia.org/wiki/Associated_Legendre_polynomials
  • http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html
• functions PQ
  • https://en.wikipedia.org/wiki/Legendre_function
  • http://mathworld.wolfram.com/LegendreFunctionoftheFirstKind.html
  • http://mathworld.wolfram.com/LegendreFunctionoftheSecondKind.html

See also http://ask.sagemath.org/question/27230/problem-with-hypergeometric/

Component: symbolics

Keywords: orthogonal

Author: Ralf Stephan, Stefan Reiterer

Branch/Commit: 7238198

Reviewer: Marc Mezzarobba, Travis Scrimshaw

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

[sagetrac-maldun]

comment:1

Hi!

Good to see that someone else is working on making the orthogonal polynomials symbolic, since my research interests shifted heavily in the past years.

A good read on Legendre polynomials is also the bible for ortho polys: Abramowitz and Stegun http://people.math.sfu.ca/~cbm/aands/page_331.htm

There you will find also much information on special values and other properties.
I currently have some time left and can help with one thing or the other if you like,

[rwst]

comment:2

Fine! See also http://trac.sagemath.org/wiki/symbolics/functions

[rwst]

comment:3

OK I have a question. What is the equivalent recursive algorithm to https://github.com/sagemath/sage-prod/blob/master/src/sage/functions/orthogonal_polys.py#L812-834
for Legendre polynomials?

The link is valid as long #16812 is not merged.

[sagetrac-maldun]

comment:4

Replying to @rwst:

OK I have a question. What is the equivalent recursive algorithm to https://github.com/sagemath/sage-prod/blob/master/src/sage/functions/orthogonal_polys.py#L812-834
for Legendre polynomials?

The link is valid as long #16812 is not merged.

Hi!

You won't have luck to find an equivalent recursion algorithm for Legendre Polynomials, since the recursion algorithm for Chebyshev Polynomials uses the fact that cheby polynomials are cosines in disguise, and thus one is able to build Cheby polyis in O(log N) time. For Legendre polynomials you have to use the classic recursion formula given in https://en.wikipedia.org/wiki/Legendre_polynomials#Recursive_definition

[sagetrac-maldun]

comment:5

Maybe this could help you:

I already implemented all Orthopolys one time: https://github.com/sagemath/sage-prod/files/10650394/trac_9706_ortho_polys.patch.gz

they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)

[rwst]

comment:6

Some timings for P(n,z):

sage: legendre_P(100,2.5)
6.39483750487443e66
sage: timeit('legendre_P(100,2.5)')
25 loops, best of 3: 21 ms per loop

sage: from mpmath import legenp 
sage: legenp(100,0,2.5)
mpf('6.3948375048744286e+66')
sage: timeit('legenp(100,0,2.5)')
625 loops, best of 3: 97.2 µs per loop

sage: from scipy.special import eval_legendre
sage: eval_legendre(int(100),float(2.5))
6.3948375048744324e+66
sage: timeit('eval_legendre(int(100),float(2.5))')
625 loops, best of 3: 7.62 µs per loop

sage: eval_legendre(int(10^5),float(1.00001))
3.1548483029540554e+192
sage: timeit('eval_legendre(int(10^6),float(2.5))')
25 loops, best of 3: 11.8 ms per loop
sage: eval_legendre(int(10^6),float(2.5))
inf

while legenp will already bail out at 10⁵ because of F convergence issues.

[rwst]

comment:7

With P(n,x) symbolics and algebra, Pari is much better than Maxima

sage: P.<t> = QQ[]
sage: timeit('legendre_P(1000,t)')
5 loops, best of 3: 2.8 s per loop
sage: timeit('pari.pollegendre(1000,t)')
625 loops, best of 3: 366 µs per loop

[rwst]

Branch: u/rws/symbolic_legendre___associated_legendre_functions___polynomials

[rwst]

Commit: 0f86b77

[rwst]

comment:9

This is a proof of concept patch, and one can already use legendre_P and see from that and the code how the other three functions will look like. So, now would be a good time for fundamental criticism 8)

---

New commits:

50da8c5	16813: skeleton P(n,x)
e74539b	16813: P(n,x) refined, documentation
0f86b77	16813: fixes for doctest failures

[rwst]

comment:10

Replying to @sagetrac-maldun:

A good read on Legendre polynomials is also the bible for ortho polys: Abramowitz and Stegun http://people.math.sfu.ca/~cbm/aands/page_331.htm

This appears outdated, it is replaced by http://dlmf.nist.gov/14

[sagetrac-maldun]

comment:11

Replying to @rwst:

Replying to @sagetrac-maldun:

A good read on Legendre polynomials is also the bible for ortho polys: Abramowitz and Stegun http://people.math.sfu.ca/~cbm/aands/page_331.htm

This appears outdated, it is replaced by http://dlmf.nist.gov/14

You can't call a source outdated, which still covers information that the newer source doesn't. I checked your link, and some things from A&S are missign e.g. explicit representation of Legendre Polynomials with their polynomial coefficients.
And on another note: I don't see much harm in citing a classic work on this topic ...

[rwst]

comment:12

Replying to @sagetrac-maldun:

I already implemented all Orthopolys one time: https://github.com/sagemath/sage-prod/files/10650394/trac_9706_ortho_polys.patch.gz

they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)

I am not sure about the derivatives. For P(3,2,x).diff(x) I get -45*x^2 + 15 (Wolfram agrees) while with your formula (lines 2377-2395 of the patch) I get (after simplification) -45*x^2 - 15.

Update: what's your reference there?

[sagetrac-maldun]

comment:13

Replying to @rwst:

Replying to @sagetrac-maldun:

I already implemented all Orthopolys one time: https://github.com/sagemath/sage-prod/files/10650394/trac_9706_ortho_polys.patch.gz

they only would need cleanup/restructuring. Maybe you can reuse some of the implemented methods (like recursions and derivatives)

I am not sure about the derivatives. For P(3,2,x).diff(x) I get -45*x^2 + 15 (Wolfram agrees) while with your formula (lines 2377-2395 of the patch) I get (after simplification) -45*x^2 - 15.

Update: what's your reference there?

It seems you are right. from Gradshteyn-Ryzhik p.1004 formula 8.731-1 we have the relation

P(n,m,x).diff(x) = ((n+1-m)*P(n+1,m,x)-(n+1)*x*P(n,m,x))/(x**2-1)

The same relation holds for gen_legendre_Q

I suppose that's an copy/paste/rewrite mistake from my side.

[rwst]

comment:14

Also, your recursive functions for Q(n,x) and Q(n,m,x) appear to be wrong:

sage: legendre_Q.eval_recursive(2,x).subs(x=3)
13/2*I*pi + 13/2*log(2) - 9/2
sage: legendre_Q.eval_recursive(2,x).subs(x=3).n()
0.00545667363964419 + 20.4203522483337*I
sage: legendre_Q(2,3.)
0.00545667363964451 - 20.4203522483337*I

The latter result from mpmath is supported by Wolfram.

As to Q(n,m,x):

sage: gen_legendre_Q(2,1,x).subs(x=3)
-1/8*sqrt(-2)*(72*I*pi + 72*log(2) - 50)
sage: gen_legendre_Q(2,1,x).subs(x=3).n()
39.9859464434253 + 0.0165114736149170*I
sage: gen_legendre_Q(2,1,3.)
-39.9859464434253 + 0.0165114736149193*I

Again, Wolfram supports the latter value from mpmath (symbolic as (25 i)/(2 sqrt(2))-18 i sqrt(2) ((log(4))/2+1/2 (-log(2)-i pi))).

[rwst]

comment:15

OK, I resolved it by using conjugate() on every logarithm in the Q(n,x) algorithms (on which the Q(n,m,x) recurrence is based, too).

Update: it however makes symbolic work tedious and differentiation impossible, at the moment.

See also https://groups.google.com/forum/?hl=en#!topic/sage-support/bEMPMEYeZKU on derivatives of conjugates in Sage.

[sagetrac-maldun]

comment:16

Thanks for resolving this issue! I suppose I wasn't careful enough with complex arguments. But in my defense: I hadn't time to test this codes well enough when I wrote them ... but hopefully they give some useful informations.

concerning complex conjugation: I hope my answer on the mailing list give some clues: https://groups.google.com/forum/?hl=en#!topic/sage-support/bEMPMEYeZKU

Replying to @rwst:

OK, I resolved it by using conjugate() on every logarithm in the Q(n,x) algorithms (on which the Q(n,m,x) recurrence is based, too).

Update: it however makes symbolic work tedious and differentiation impossible, at the moment.

See also https://groups.google.com/forum/?hl=en#!topic/sage-support/bEMPMEYeZKU on derivatives of conjugates in Sage.

[rwst]

comment:17

Replying to @sagetrac-maldun:

Thanks for resolving this issue!

Unfortunately, while it would be easy to resolve this numerically, the instances of conjugate() introduced in the recurrence will make symbolic results from the recurrence unreadable and, in case of derivatives, impossible to use. A different way of computing the recurrences is needed, one which does away with usage of conjugate().

[sagetrac-maldun]

comment:18

Hi!

We have several possible ways out of this:

1. avoid recursion for symbolic argument for Legendre_Q and use another library (maxima, flint, sympy ... ) for evaluation.
2. Let it be, but avoid it as default method.
3. Maybe more elegantly: There is a more closed relation for legendre_Q (but it's not really a recursion):

Q(n,z) = ½P(n,z) ln((z+1)/(z-1)) - W(n-1,z)

with

W(n,z) = Σ_{k=1}^n (1/k) P(k-1,z) P(n-k,z)

(Gradshteyn-Ryzhik p 1019f)

Maybe we could find an recursion for W(n,z)

Hope this could be of some use

Edit there should be a recursion since W(n,z) is the convolution of aseries of holonomic functions, and if I remeber correctly there is an theorem saying that convolutions of holonomic functions are also holonomic, thus should have a recursion.

Edit: The above mentioned relation can also be found here: http://people.math.sfu.ca/~cbm/aands/page_334.htm

[rwst]

comment:19

sage: from ore_algebra import *
sage: def W(n):
    return sum([(1/k)*legendre_P(k-1,t)*legendre_P(n-k,t) for k in range(1,n+1)])
....: 
sage: R.<t> = QQ['t']
sage: guess([W(n) for n in range(1,10)], OreAlgebra(R['n'], 'Sn'))
(-n - 3)*Sn^2 + (2*t*n + 5*t)*Sn - n - 2

[sagetrac-maldun]

comment:20

Nice! I didn't know that sage already supports Ore Algebras.
It appears that my holonomic function package on Mathematica stopped to work.

[rwst]

comment:21

I always need some time to figure out the final form (that offset of n!) but: n*W_n = (2tn-t)*W_n-1 - (n-1)W_n-2 (W₀=0, W₁=1).

[rwst]

comment:22

However, this will yield the same result unless the log has conjugate associated with it. This shows your recurrence is actually correct but numerical results derived from it by substitution may need a warning about the log branch. I know not enough about calculus, maybe I'll ask again, this time on sage-devel, if there should be a switch for log() to select the branch in case of numerical evaluation. What do you think?

[rwst]

comment:23

I think there must be another different formula, because Wolfram has this for Q(2,x):

sage: ((3*x^2-1)/2*(log(x+1)-log(1-x))/2-3*x/2).subs(x=3)
-13/2*I*pi + 13/2*log(4) - 13/2*log(2) - 9/2
sage: ((3*x^2-1)/2*(log(x+1)-log(1-x))/2-3*x/2).subs(x=3).n()
0.00545667363964419 - 20.4203522483337*I

which yields the correct value without use of conjugate.

The first few Q(n,x) from Wolfram are:

Q(0,x) = 1/2 log(x+1)-1/2 log(1-x)
Q(1,x) = x (1/2 log(x+1)-1/2 log(1-x))-1
Q(2,x) = 1/2 (3 x^2-1) (1/2 log(x+1)-1/2 log(1-x))-(3 x)/2
Q(3,x) = -(5 x^2)/2-1/2 (3-5 x^2) x (1/2 log(x+1)-1/2 log(1-x))+2/3

which makes clear that instead of log((1+x)/(1-x)).conjugate() we should just use log(1+x)-log(1-x), of course. Oh well.

[sagetrac-maldun]

comment:24

Oh yeah it's again the non uniqueness of the representation of the complex logarithm

sage: log((x+1)/(1-x)).subs(x=3)
I*pi + log(2)
sage: (log(x+1)-log(1-x)).subs(x=3).simplify_log()
-I*pi + log(2)
sage: log((x+1)/(1-x)).subs(x=3).conjugate()
-I*pi + log(2)

confusing as hell ...

I think Wolfram uses the log(1+x)-log(1-x) representation simply by the fact that it is independent of the branch in the following sense:
Let log(x) = ln|x| + iarg(x) + 2kπi and log(y) = ln|y| + iarg(y) + 2kπi
then

log(x) - log(y) = ln|x| + i*arg(x) + 2kπi- ln|y| + i*arg(y) + 2kπi = 
= ln|x/y| + i*(arg(x) - arg(y)) + 0 

I.e. if we have the same branch on the logarithm the module of 2kπi cancels out.

That means the formula isn't exactly wrong, it uses simply a different branch of the logarithm. But the representation of log as difference saves us indeed a lot of trouble, and as showed above is independent of the branch we use.

Nevertheless I think we should stick to the recursion with W(n,x), because from a computational view it is a lot better since:

1. The computational complexity is the same (solving a two term recursion)

2. we save computation time since we don't have to simplify expressions containing logarithms but only polynomials which are much simpler to handle and expand.

[sagetrac-git]

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

0b67667	16813: Q(n,x) implementation
06df7c8	16813: implement P(n,m,x)
7bd4a70	Merge branch 'develop' into t/16813/symbolic_legendre___associated_legendre_functions___polynomials
74ca8ea	16813: implement Q(n,m,x); fixes, doctests

[sagetrac-git]

Changed commit from 60b8649 to 7238198

[sagetrac-git]

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

cfccbb7	Merge branch 'develop' into t/16813/16813-1
7238198	fix doctest

[rwst]

comment:90

Review needed for the last two commits.

[tscrim]

Changed reviewer from Marc Mezzarobba to Marc Mezzarobba, Travis Scrimshaw

[tscrim]

comment:92

One failing doctest, see patchbot.

[rwst]

comment:93

Does not fail here. Note the patchbot is running 7.5beta1 not 2, and there were gegenbauer changes with pynac-0.7.0. I'll try to trigger the bots.

[tscrim]

comment:94

I ended up accidentally trying this on my beta1 as well. All is good.

[vbraun]

comment:95

The dependency is invalid/wontfix. What are you depending on?

[rwst]

comment:96

Replying to @vbraun:

The dependency is invalid/wontfix. What are you depending on?

Actually I don't recall and couldn't find out how floor(x,hold=True) was fixed. That bug was the reason for the random_expr() doctest fail in this and some other ticket. I'll remove the dependency.

[rwst]

Changed dependencies from #19464 to none

[vbraun]

Changed branch from u/rws/16813-1 to 7238198