# Series expansion in terms of Hermite polynomials

I am trying to expand a polynomial in terms of orthogonal polynomials (in my case, Hermite). Maple has a nice built-in function for this, ChangeBasis.

Is there a similar function in Mathematica? And if not, where should I look for the algorithm?

-

The Hermite polynomials are orthogonal with respect to the inner product $$\langle f,g \rangle = \int_{-\infty}^{\infty} f(x)g(x)e^{-x^2} \, dx.$$ Thus, the nth coefficient can be computed using the inner product of your polynomial with the nth normalized Hermite polynomial.

Example:

p[x_] = 1 + x + x^2 + x^3;
coeffs = Table[
Integrate[HermiteH[n, x]*p[x]*Exp[-x^2], {x, -Infinity, Infinity}]/
Integrate[HermiteH[n, x]^2*Exp[-x^2], {x, -Infinity, Infinity}],
{n, 0, 3}]
(* Out: {3/2, 5/4, 1/4, 1/8} *)

coeffs.Table[HermiteH[n, x], {n, 0, 3}] // Expand
(* Out: 1 + x + x^2 + x^3 *)

-

For polynomials, you don't need to do any integrals to find the expansion. Take a polynomial p and a list basis containing the basis functions. Then define a function that takes these two, identifies the variable x, and solves for the coefficients in basis that make the two polynomials equal in terms of their CoefficientLists:

expandPoly[p_, basis_, x_] :=
# /. First@Solve[CoefficientList[#.basis, x] == #2, #] &[
Array["a", Length[#]], #] & @

expandPoly[1 + x + 3 x^2 + 7 x^3, HermiteH[Range[4] - 1, x], x]

(* ==> {5/2, 23/4, 3/4, 7/8} *)


Edit

In response to belisarius: if you already know that you're only interested in a basis of HermiteH, you could incorporate that into the function and do away with the specification of the variable basis as follows:

expandPoly[p_, x_] := # /.
First @ Solve[
CoefficientList[#.HermiteH[Range[Length[#]] - 1, x],
x] == #2, #] &[Array["a", Length[#]], #] & @ CoefficientList[p, x]

expandPoly[1 + x + 3 x^2 + 7 x^3, x]

(* ==> {5/2, 23/4, 3/4, 7/8} *)


Edit 2

With the general function given as the first solution above, you can specify any set of polynomials that is known to form a basis for degree n or larger. This means the basis functions don't have to be orthogonal polynomials at all.

-
(+1) Very nice. I was just doing similar thing with another problem, but didn't think of it here. –  Michael E2 Apr 14 '13 at 22:02
How to get rid of that Range@4? –  belisarius Apr 14 '13 at 22:08
@belisarius See edit - I hardcoded the HermiteH part into the function in case you're only interested in those basis functions. Then the counting of terms is automatic. In the first version, I wanted to keep the list basis deliberately general so you can use other polynomials there, too. –  Jens Apr 14 '13 at 22:13
Ahh... Algebra. I accidentally studied analysis in graduate school. :) –  Mark McClure Apr 14 '13 at 22:17
If the expansion is known finite, one should be able to manage the upper bound automagically. –  belisarius Apr 14 '13 at 23:08

The inner product for the Hermite polynomials, $$\langle f, g\rangle \int_{-\infty}^{\infty} f(x)\,g(x)\,e^{-x^2}\;dx\,,$$ has nice formulas for power functions (where $n=a+b$) and for the Hermite polynomials: \begin{align} \langle x^a, x^b \rangle = \langle x^n, 1\rangle &= \frac{1}{2} \left((-1)^n+1\right)\, \Gamma \left(\frac{n+1}{2}\right)\cr \langle H_n(x), H_n(x) \rangle &= \sqrt{\pi}\,2^n n! \cr \end{align}

These can be used to give a quick change of basis function for polynomials.

hermiteIP[f_, g_, x_] := With[{coeff = CoefficientList[f g, x]},
coeff.Table[1/2 (1 + (-1)^(-1 + n)) Gamma[n/2], {n, Length@coeff}]];

hermiteExpand[poly_, var_] /; PolynomialQ[poly, var] :=
Sum[hermiteIP[poly, HermiteH[n, var], var] H[n, var]/(Sqrt[Pi] 2^n n!),
{n, 0, Exponent[poly, var]}]


I used H[n, x] as a place holder for HermiteH[n, x].

hermiteExpand[(1 + x)^5, x]
(* 39/4 H[0, x] + 95/8 H[1, x] + 25/4 H[2, x] + 15/8 H[3, x] +
5/16 H[4, x] + 1/32 H[5, x]  *)

hermiteExpand[(1 + x)^5, x] /. H -> HermiteH
(* 39/4 H[0, x] + 95/8 H[1, x] + 25/4 H[2, x] + 15/8 H[3, x] +
5/16 H[4, x] + 1/32 H[5, x]  *)

% // Factor
(* (1 + x)^5 *)

-
Probably the best answer. It uses the orthogonality and speeds things up using the closed form. +1 –  Mark McClure Apr 14 '13 at 22:24

Whenever I want to convert some polynomial expressed with respect to a certain basis in terms of another polynomial basis. my go-to algorithm is Salzer's algorithm. It's rather fast, since it relies only on recurrences. Here's a specialization of that algorithm for the case of monomial-Hermite conversion:

monomialToHermite[cofs_?VectorQ] := Module[{n = Length[cofs] - 1, a},
a[0, 0] = cofs[[n + 1]]; a[0, 1] = cofs[[n]];
a[1, 1] = cofs[[n + 1]]/2;
Do[
a[0, k + 1] = cofs[[n - k]] + a[1, k];
Do[
a[m, k + 1] = (m + 1) a[m + 1, k] + a[m - 1, k]/2,
{m, k - 1}];
a[k, k + 1] = a[k - 1, k]/2; a[k + 1, k + 1] = a[k, k]/2,
{k, n - 1}];
Table[a[m, n], {m, 0, n}]]


The algorithm as I presented it here uses an implicit two-dimensional array, a, to clearly show off the recurrence. The algorithm can be easily rewritten so that it uses only a pair or so of one-dimensional arrays, but I'll leave out that version for now.

Here's a test of Salzer's method:

monomialToHermite[{1, 1, 3, 7}]
{5/2, 23/4, 3/4, 7/8}

{1, 1, 3, 7}.x^Range[0, 3] == {5/2, 23/4, 3/4, 7/8}.HermiteH[Range[0, 3], x] // Expand
True

CoefficientList[(1 + x)^5, x] // monomialToHermite
{39/4, 95/8, 25/4, 15/8, 5/16, 1/32}

%.HermiteH[Range[0, 5], x] == (1 + x)^5 // Expand
True


(Other instances where I used Salzer's algorithm include this and this.)

-
(Of course, the algorithm is easily implemented even in languages that don't have symbolic capabilities.) –  Ｊ. Ｍ. Apr 15 '13 at 0:44
Late answers are usually underestimated, +2. –  Artes Apr 15 '13 at 19:29
Can't argue with that. Thanks @Artes! –  Ｊ. Ｍ. Apr 16 '13 at 5:13

As another variation, here's another method based on repeated greedy division:

Reap[Fold[Block[{q, r}, {q, r} = PolynomialQuotientRemainder[#1, #2, x]; Sow[q]; r] &,
x^3 + x^2 + x + 1, HermiteH[Range[3, 0, -1], x]]][[-1, 1]] // Reverse
{3/2, 5/4, 1/4, 1/8}


Check:

%.HermiteH[Range[0, 3], x] == x^3 + x^2 + x + 1 // Expand
True

-