Where can I find the most widely adopted implementation of Legendre polynomials?

Apart from an implementation hidden somewhere in Jacobi, I could not find where this special, but quite common, function was implemented. Am I missing the obvious location? Maybe it belongs in SpecialFunctions?

2 Likes

do either of these help?

https://github.com/milthorpe/SphericalHarmonics.jl

3 Likes

I wouldn’t call the function legendre listed in the README of https://github.com/pjabardo/Jacobi.jl “hidden away”.

2 Likes

just re-found that – nice package … I added it above

2 Likes

It requires unnecessary allocations, but you can also use ApproxFun:

f = Fun(Legendre(), [zeros(k-1);1])
f(0.1)
1 Like

Thanks everyone for the overview. SphericalHarmonics comes closest to where I expected legendre to reside. In striving for a minimal amount of dependencies I would be a bit adverse to pulling in large complex packages like GSL and Jacobi if I only wanted a single function from them.

1 Like

Is there now a standard package for associated Legendre polynomials in Julia? This package SphericalHarmonics.jl has not been updated in over 3 years, and the relationship between associated Legendre polynomials and Jacobi polynomials is not exactly straightforward.

Associated Legendre polynomials are necessary for Spherical Harmonics, which are a corner stone of the mathematics of waves.

3 Likes

There’s an implementation buried here: https://github.com/JuliaApproximation/FastTransforms.jl/blob/master/src/SphericalHarmonics/sphfunctions.jl

1 Like

Thanks @dlfivefifty, that is quite buried!

For now, I’ve gone with GSL.jl. Am writing some basic tests in my package EffectiveWaves.jl if anyone wants to see how to use GSL.jl

1 Like

The easiest to use implementation I’ve found is in CMB.jl at https://github.com/jmert/CMB.jl/blob/master/src/legendre.jl , also it looks to be well written.

1 Like

I’ve cleaned it up and put it into one of my codes, but was thinking of merging it back once I’ve optimised and stabilised it. Let me know if there is any interest. This would be a much more general Spherical Harmonics package though which also computes cg-coefficients, etc.

3 Likes

interested

Maybe it’s a bit late, but if you are interested in evaluating Legendre polynomials (and other classical orthogonal polynomials), you can also use my package PolynomialBases.jl. There is a brief tutorial given as Notebook, showing how to evaluate some orthogonal polynomials. However, it does not contain general associated Legendre polynomials.

3 Likes

To add to the noise: there is the work-in-progress ContinuumArrays.jl which builds on QuasiArrays.jl. The idea is to treat bases as “quasi-arrays” where the first axis is continuous. These means we can bootstrap array terminology:

julia> using ContinuumArrays

julia> P = Legendre(); # Legendre polynomials

julia> size(P) # uncountable ∞ x countable ∞
(ContinuumArrays.AlephInfinity{1}(), ∞)

julia> axes(P) # essentially (-1..1, 1:∞)
(QuasiArrays.Inclusion{Float64,DomainSets.ChebyshevInterval{Float64}}(-1.0..1.0 (Chebyshev)), OneToInf())

julia> P[0.1,1:10] # [P_0(0.1), …, P_9(0.1)]
10-element Array{Float64,1}:
  1.0                
  0.1                
 -0.485              
 -0.14750000000000002
  0.3379375          
  0.17882875         
 -0.2488293125       
 -0.19949294375000004
  0.180320721484375  
  0.21138764183593753

julia> @time P[range(-1,1; length=10_000), 1:10_000]; # construct 10_000^2 Vandermonde matrix
  1.624796 seconds (10.02 k allocations: 1.491 GiB, 6.81% gc time)

This also works for associated Legendre polynomials as weighted Ultraspherical polynomials:

julia> associatedlegendre(m) = ((-1)^m*prod(1:2:(2m-1)))*(UltrasphericalWeight((m+1)/2).*Ultraspherical(m+1/2))
associatedlegendre (generic function with 1 method)

julia> associatedlegendre(2)[0.1,1:10]
10-element Array{Float64,1}:
   2.9699999999999998
   1.4849999999999999
  -6.9052500000000006
  -5.041575          
  10.697754375       
  10.8479361375      
 -13.334647528125    
 -18.735466024687497 
  13.885467170308594 
  28.220563705988674 

It’s very experimental but happy to have help. The long term aim is to replace the internals of ApproxFun to give better performance for spherical harmonics, etc.

6 Likes

Does this use recursion to save time when values for P_l at consecutive l are asked for?

Yes of course: it would never make to 10k otherwise.

2 Likes

Also interested!

Hah, I didn’t see this topic come through, but that’s my package! Thanks for the kind words!

I’m glad somebody has found it useful. Unfortunately, I haven’t made any improvements to it in quite some time — I’m trying to finish my PhD, and Julia is only a hobby for me, so development has taken a back seat until I graduate.

1 Like