compute hermite polynomials

[Andrei Berceanu]

Hi,

Are Hermite polynomials (http://en.wikipedia.org/wiki/Hermite_polynomials) implemented in Julia? Is there an easy way to compute Hn(x)?

//A

[Jiahao Chen]

>  Is there an easy way to compute Hn(x)?

Do you mean to evaluate a given Hermite polynomial of order n at a value x?

[Andrei Berceanu]

Yes, exactly, in order to generate plots like http://en.wikipedia.org/wiki/Hermite_polynomials#mediaviewer/File:Hermite_poly_phys.svg

//A

[Jiahao Chen]

All you need to do is compute the polynomials using the recurrence relations in the Wikipedia article. There are more sophisticated ways to evaluate Hermite polynomials but for plotting a few of the lower-order polynomials computing them directly from the definition should suffice.

Thanks,

Jiahao Chen
Staff Research Scientist
MIT Computer Science and Artificial Intelligence Laboratory

[Andras Niedermayer]

I was looking for Hermite polynomials and haven't found any code. I have some (very unpolished) code.

I haven't made a public package yet, since it needs to be improved (especially in terms of efficiency, also documentation). Unfortunately, I'm unlikely to have time for this in the near future, so I'll just post a link to a gist:

https://gist.github.com/afniedermayer/57873094430e8ddb201c

I mainly used it with the output of the ODE.jl.

I hope this is a useful starting point...

Best,

Andras

[Andrei Berceanu]

Jiahao, it's not just for reproducing the Wikipedia figure. I will need to compute higher orders as well, i.e. Hn for n = 48. I was just wondering if there was anything already implemented in Julia. Meanwhile I found this: https://github.com/daviddelaat/Orthopolys.jl

[Jiahao Chen]

The entire point of orthogonal polynomials is that they can be evaluated very efficiently using their recurrence relations. I would make use of this fact to the extent possible allowed by numerical roundoff.

[Andras Niedermayer]

Sorry, I meant Cubic Hermite Interpolation. Now I see you're looking for Hermite polynomials.

On Monday, February 2, 2015 at 4:50:00 PM UTC+1, Andras Niedermayer wrote:

[Andrei Berceanu]

Andras, no worries :) Now I understand why I couldn't find the polynomials in your gist!

//A

[Andrei Berceanu]

I came up with this, so far

function compute_hermite_polynomial(n)
    P = Poly([1])
    const x = Poly([0; 1])                                                                                
    for i = 1:n
        P = 2x*P - polyder(P)
    end
    P
end

[Jiahao Chen]

This code computes the coefficient of the leading order term, not the value of the polynomial.

Thanks,

Jiahao Chen
Staff Research Scientist
MIT Computer Science and Artificial Intelligence Laboratory

[Andrei Berceanu]

That depends on weather one uses the Polynomial or Polynomials package. I use the latter (which has in fact superseded the former). The main difference is the ordering of the coefficients.

[Steven G. Johnson]

On Monday, February 2, 2015 at 11:17:36 AM UTC-5, Jiahao Chen wrote:

The entire point of orthogonal polynomials is that they can be evaluated very efficiently using their recurrence relations. I would make use of this fact to the extent possible allowed by numerical roundoff. 

(In a serious package, I would tend to want to see additional tricks, e.g hard-coding the low-order base cases with @evalpoly, combining multiple recurrence stages to coarsen the loop body, etc.)

[Jiahao Chen]

Thanks for catching that. But even with Polynomials, the code computes the coefficients of the polynomial in the monomial basis, not the value of the polynomial at a given x. Using the former to compute the latter is dangerous because it is numerically unstable.

julia> Hermite48(x) = 20007974164906320568399715106816000000-960382759915503387283186325127168000000x^2+7362934492685525969171095159308288000000x^4-21597941178544209509568545800637644800000x^6+32396911767816314264352818700956467200000x^8-28797254904725612679424727734183526400000x^10+16580237672417776997244540210590515200000x^12-6559214903374065625283554369024819200000x^14+1858444222622651927163673737890365440000x^16-388694216496240925942729147794063360000x^18+61372771025722251464641444388536320000x^20-7439123760693606238138356895580160000x^22+700787020934904935476801736540160000x^24-51750426161346826004440743621427200x^26+3011929564946111566396022115532800x^28-138479520227407428340046993817600x^30+5025466459865592157501705420800x^32-143328811689571255828925644800x^34+3185084704212694573976125440x^36-54368444453132766554357760x^38+697031339142727776337920x^40-6476481664508504309760x^42+41077050726269583360x^44-158751886864809984x^46+281474976710656x^48

Hermite48 (generic function with 1 method)

julia> Hermite48(2.0)

1.13060295173709763271490845506831122432e+38 with 256 bits of precision

julia> Hermite48(2) #Integer overflow even with BigInt

-42115274364030603711130058827879640052269056

[Jiahao Chen]

And in fact the coefficients you obtained with the code are not correct because of integer overflow. (I got these coefficients from Mathematica.)

julia> compute_hermite_polynomial(48)

Poly(8629227277223198720 - 8374539685103403008x^2 + 2715657340094251008x^4 + 6791467061357838336x^6 + 8259543481672794112x^8 + 9055289415143784448x^10 + 376270965132230656x^12 + 1067414392280186880x^14 + 8920937959042056192x^16 - 4036022274615148544x^18 + 6462554277163302912x^20 - 783339912383430656x^22 + 7024739932564357120x^24 + 6859947603694452736x^26 - 5181745414033899520x^28 - 5062547358466703360x^30 - 634480881389535232x^32 + 5903955228692054016x^34 - 3820547819823955968x^36 - 5156551204595040256x^38 + 2667573538658975744x^40 - 1674494636451692544x^42 + 4183562578850480128x^44 - 158751886864809984x^46 + 281474976710656x^48)