Series expansion in terms of Hermite polynomials

Question

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?

Answer 1

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[#]], #] & @
     PadRight[CoefficientList[p, x], Length[basis]]

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.

Answer 2

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} \, \mathrm dx.$$ Thus, the $n$-th coefficient can be computed using the inner product of your polynomial with the $n$-th normalized Hermite polynomial.

Example:

p[x_] = 1 + x + x^2 + x^3;
coeffs = Table[
  Integrate[HermiteH[n, x]*p[x]*Exp[-x^2], {x, -∞, ∞}]/
   Integrate[HermiteH[n, x]^2*Exp[-x^2], {x, -∞, ∞}],
  {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 *)

Answer 3

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 *)

Answer 4

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.)

Answer 5

SolveAlways can find coordinates of a polynomial with respect to any given basis. You set up an equation, setting the given polynomial equal to a linear combination of your basis polynomials. This approach will work generally with any polynomial that is a linear combination of a given set of (linearly independent) polynomials.

poly = 1 + x + 3 x^2 + 7 x^3;       (* Jens' example *)

params = Table[C[i], {i, 0, Exponent[poly, x]}];
basis = HermiteH[Range[0, Exponent[poly, x]], x];

coeff = params /. First@SolveAlways[poly == params.basis, x]
(*  {5/2, 23/4, 3/4, 7/8}  *)

Update: Extension to ChangeBasis

Here is a version of the Maple ChangeBasis that works on polynomials (but not on infinite sums). It is based on the above idea.

ClearAll[changeBasis];
changeBasis[poly_, basisfamilies : {_, _} ..] := 
 Module[{vars, families, degrees, params, basis, coeff},
  vars = {basisfamilies}[[All, 1]];
  families = {basisfamilies}[[All, 2]];
  degrees = Exponent[poly, vars];
  params = Outer[C, ##] & @@ (Range[0, #] &) /@ degrees;
  basis =
   With[{bfn = Times @@ MapThread[#1[#3, #2] &,
                               {families, vars, Array[Slot, Length[vars]]}]},
    Outer[HoldForm[bfn] &, ##] & @@ (Range[0, #] &) /@ degrees
    ];
  coeff = params /. First@SolveAlways[
      poly == Total[params*ReleaseHold@basis, -1], vars];
  With[{res = Total[coeff*basis, -1]},
   Hold[res] /. HoldForm[v_] :> v
   ]
  ]

Some examples from the Maple docs. GegenbauerC[n, 1, x] is automatically simplified to ChebyshevU[n, x], which affect the form of the output.

changeBasis[1 + 2 x + 3 x^3, {x, GegenbauerC[#1, 1, #2] &}]
ReleaseHold[%] // Expand
GegenbauerC[n, 1, x]        (* Same as ChebyshevU[n, x] *)
(* 
  Hold[ChebyshevU[0, x] + 7/4 ChebyshevU[1, x] + 3/8 ChebyshevU[3, x]]

  1 + 2 x + 3 x^3
  ChebyshevU[n, x]
*)

Multivariable function:

changeBasis[y^2 + x^2 - 1, {x, LaguerreL[#1, 1, #2] &}, {y, ChebyshevU}]
ReleaseHold[%] // Expand
(*
  Hold[21/4 (ChebyshevU[0, y] LaguerreL[0, 1, x]) + 
      1/4 ChebyshevU[2, y] LaguerreL[0, 1, x] - 
      6 (ChebyshevU[0, y] LaguerreL[1, 1, x]) + 
      2 (ChebyshevU[0, y] LaguerreL[2, 1, x])]

  -1 + x^2 + y^2
*)

Answer 6

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