Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question

5 Answers 5

up vote 9 down vote accepted

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 *)
share|improve this answer

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.

share|improve this answer
    
(+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 *)
share|improve this answer
    
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.)

share|improve this answer
    
(Of course, the algorithm is easily implemented even in languages that don't have symbolic capabilities.) –  J. M. 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! –  J. M. 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
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.