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

Specific Domains Numerics
#1

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
#2

do either of these help?

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

3 Likes
#3

I wouldnā€™t call the function legendre listed in the README of https://github.com/pjabardo/Jacobi.jl ā€œhidden awayā€.

2 Likes
#4

just re-found that ā€“ nice package ā€¦ I added it above

2 Likes
#5

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

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

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
#7

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
#8

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

1 Like
#9

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
#10

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
#11

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
#12

interested

#13

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
#14

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
#15

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

#16

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

2 Likes
#17

Also interested!

#18

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