haskellmathrecursionhermite

Haskell Hermite polynomials implementation

Haskell allows to represent recurrent functions in a very concise way. For example, infinite list, that contains Fibonacci numbers can be defined as follows:

fibs :: [Integer]
fibs = 1 : 1 : zipWith (+) fibs (tail fibs)

I am dealing with 'probabilists' Hermite polynomials, which have the following recursion relation:

enter image description here

What would be the optimal way to construct the infinite list of n-th Hermite polynomials for given x?

Solution

  • We can write it as:

    hermite :: (Enum n, Num n) => n -> [n]
hermite x = s
    where s@(_:ts) = 1 : x : zipWith3 (\hn2 hn1 n1 -> x*hn1 - n1*hn2) s ts [1..]

    where the first items 1 : x : ... are the first elements of the hermite (you can fill in other values).

    For the next one, we zip the original values s (so that starts with H0), the tail ts of s (that starts with H1) and the index (that starts with 2, 3, ...) and perform an operation x*hn1 - x*hn2 (nh1 stands for Hn-1, and nh2 stands for Hn-2), and so we calculate the next element each time.

    The first 11 values for x = 0.75 are:

    Prelude> take 11 (hermite 0.75)
[1.0,0.75,-0.4375,-1.828125,-5.859375e-2,7.2685546875,5.744384765625,-39.30303955078125,-69.68797302246094,262.1583366394043,823.8105096817017]

    So the first value is 1, the second x, the third one x*x-2, the fourth one x*x*x-2*x-3*x, and so on.

    That being said, if I recall correctly, the recursion formula of the Hermite polynomials is:

    Hn(x) = 2×x×Hn-1(x)-2×(n-1)Hn-2(x)

    instead of the one quoted in the question.

    In that case the formula is thus:

    hermite :: (Enum n, Num n) => n -> [n]
hermite x = s
    where s@(_:ts) = 1 : 2 * x : zipWith3 helper s ts [1..]
          helper hn2 hn1 n1 = 2 * (x * hn1 - n1 * hn2)

    Then the first 11 values are:

    Prelude> take 11 (hermite 0.75)
[1.0,1.5,0.25,-5.625,-9.9375,30.09375,144.515625,-144.3515625,-2239.74609375,-1049.994140625,38740.4384765625]

    Which is correct acording to this Wolfram article:

    H0 = 1

    H1 = 2*x

    H2 = 4˙x2 - 2

    H3 = 8˙x3 - 4˙x

    H4 = 16˙x4 - 48˙x2 + 12

    Which maps exactly on the values we obtained:

    Prelude> let x = 0.75 in [1,2*x,4*x*x-2,8*x*x*x-4*x,16*x*x*x*x-48*x*x+12]
[1.0,1.5,0.25,0.375,-9.9375]