Implemented some of the classical orthogonal polynomials

[jagot]

Hi!

I implemented some simple functions for creating the classical orthogonal polynomials
through their respective recurrence relations, along with tests that they match the first
few listed on e.g. Wikipedia.

[pwl]

Hi,
Thanks, I was always missing the orthogonal polynomials from this package. Some time ago I was thinking about adding them and what came to my mind was that we could exploit the generic programming by defining e.g. Lagendre polynomials as

LegendreP(n::Int,x) = LegendrePgen(n::Int,x)[1]

function LegendrePgen{T}(n::Int,x::T)
  if n == 0
    return one(T), zero(T)
  else
    Pn1, Pn2 = LegendrePgen(n-1,x)
    return x*Pn1*(2n-1)/n-Pn2*(n-1)/n, Pn1
  end
end

Then you can evaluate the polynomial at point x immediately by calling

LegendreP(3,3.0) => 63.0

or you can generate the Polynomial by calling

LegendreP(3,Poly([1,0])) => Poly(2.5x^3 - 1.5x)

This gives you the advantage of not having to construct the polynomial type when you only want to evaluate it at a single point.

However, I had some trouble implementing the general function for generating any polynomial and I finally gave up. Do you think it would make sense to incorporate this approach to your pull request?

[jagot]

Hi,

I was thinking rather to have separate functions like,

P(n::Integer) -> Poly([…])
legendre(n::Integer, x) = polyval(P(n), x)

This way, you don’t have to pass a Poly([0,1]) each time you want the polynomial. What do you think about that?

Cheers,
Stefanos

On 16 Mar 2015, at 19:26, Paweł Biernat [email protected] wrote:

Hi,
Thanks, I was always missing the orthogonal polynomials from this package. Some time ago I was thinking about adding them and what came to my mind was that we could exploit the generic programming by defining e.g. Lagendre polynomials as

LegendreP(n::Int,x) = LegendrePgen(n::Int,x)[1]

function LegendrePgen{T}(n::Int,x::T)
if n == 0
return one(T), zero(T)
else
Pn1, Pn2 = LegendrePgen(n-1,x)
return x_Pn1_(2n-1)/n-Pn2*(n-1)/n, Pn1
end
end
Then you can evaluate the polynomial at point x immediately by calling

LegendreP(3,3.0) => 63.0
or you can generate the Polynomial by calling

LegendreP(3,Poly([1,0])) => Poly(2.5x^3 - 1.5x)
This gives you the advantage of not having to construct the polynomial type when you only want to evaluate it at a single point.

However, I had some trouble implementing the general function for generating any polynomial and I finally gave up. Do you think it would make sense to incorporate this approach to your pull request?

—
Reply to this email directly or view it on GitHub #18 (comment).

[pwl]

Yeah, passing Poly([0,1]) can be cumbersome at times, but it's easy enough to write a wrapper function around it. I also thought it would be nice to have a function returning Poly([1,0]) in the same fashion as there are one() and zero() functions for polynomials. Although I have no idea what name it should have.

The main issue I had with your approach (in fact this was also my approach at first) was that I needed only the values of some low order polynomials, and allocating space needed for Poly types at every iteration was significantly slowing my program down. I could also cache the results, but I thought the re-usability of my approach was better.

Another problem is generating the polynomials using types other than Float64. Float64 coefficients may not give you enough precision to accurately generate values for some high order polynomials. Having a function that accepts Poly{T}() as an argument would give you the flexibility of choosing the type of coefficients. I remember having to use polynomials built on top of Rational type to alleviate the lack of precision.

I think the problem with precision is even more important than the performance, especially when you give people the tool to automatically generate high order polynomials.

[jagot]

I found a way of implementing a common generator function as a macro. Specifying the polynomials then became simple one-liners. However, I have not figured out a nice way of accommodating the more general polynomials, such as the Jacobi, and the associated Legendre and Laguerre polynomials, who all need extra parameters.

[jagot]

Would you like me to add something more/change something, before this PR can be merged?