Fourier Series


This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0340. It is primarily for students who have some experience using Mathematica. If you have never used Mathematica before and would like to learn more of the basics for this computer algebra system, it is strongly recommended looking at the APMA 0330 tutorial. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0340
Return to the main page for the course APMA0330
Return to Part V of the course APMA0340
Introduction to Linear Algebra
Charles Hermite

 

The Hermite polynomials, conventionally denoted by \( H_n (x) , \) were introduced in 1859 by Pafnuty Chebyshev. Later, in 1864 they were studied by a French mathematician Charles Hermite (1822--1901). Therefore, they sometimes are called Chebyshev--Hermite polynomials.

 

Hermite Polynomials


The Hermite polynomials \( H_n (x) \) can be defined either by the exponential generating function

\[ \exp \left\{ 2xt - t^2 \right\} = \sum_{n\ge 0} H_n (x) \, \frac{t^n}{n!} \]
or by Rodrigues formula:
\[ H(x) = (-1)^n e^{x^2} \frac{{\text d}^n}{{\text d} x^n} e^{-x^2} = e^{x^2 /2} \left( x- \frac{\text d}{{\text d} x} \right)^n e^{-x^2 /2} , \qquad n=0,1,2,\ldots . \]
The Hermite polynomials can be written explicitly as
\[ H_n (x) = n! \, \sum_{k=0}^{\left\lfloor n/2 \right\rfloor} \frac{(-1)^k}{k! \left( n-2k \right)!!} \left( 2x \right)^{n-2k} , \qquad n=0,1,2,\ldots . \]
Symmetry:
\[ H_n (-x) = (-1)^n \,H_n (x) , \qquad n=0,1,2,\ldots . \]
Recurrence relation:
\[ H_{n+1} (x) = 2x \,H_n (x) - 2n\,H_{n-1} (x), \qquad n=1,2,\ldots ; \qquad H_0 (x) =1, \quad H_1 (x) =2x . \]
Derivatives:
\[ H_{n}' (x) = 2n\,H_{n-1} (x), \qquad n=1,2,\ldots . \]
Hermite polynomials are eigenfunctions of the Sturm--Liouville problem for the Hermite differential equation on the infinite interval \( (-\infty , \infty ) : \)
\[ y'' -2x\,y' + \lambda\,y =0 , \qquad \mbox{or} \qquad \frac{{\text d}}{{\text d}x} \left( e^{-x^2} \, \frac{{\text d}y}{{\text d}x} \right) + \lambda e^{-x^2}\, y =0, \]
subject to the condition that its solutions grow at infinity no faster than a polynomial, corresponding to eigenvalues \( \lambda_n = 2n \ge 0, \) an integer.

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} \,{\text d}x . \] In particular,

\[ \int_{-\infty}^{\infty} e^{-x^2} H_m (x)\,H_n (x) \,{\text d} x = \begin{cases} 2^n n! \,\sqrt{\pi} , & \quad \mbox{if $n=m$}, \\ 0 , & \quad \mbox{if $n\ne m$}. \end{cases} \]

The Hermite polynomials evaluated at zero argument are called Hermite numbers.

\[ H_n (0) = \begin{cases} (-1)^{n/2} 2^{n/2} \,(n-1)!! , & \quad \mbox{if $n$ is even}, \\ 0 , & \quad \mbox{if $n$ is odd,} \end{cases} \qquad n=0,1,2,\ldots \]

One can define the Hermite functions (wave functions):

\[ \psi_n (x) = \left( 2^n n! \,\sqrt{\pi} \right)^{-1/2} e^{-x^2 /2} \,H_n (x) \qquad n=0,1,2,\ldots \]
that form a complete orthonormal system in \( L^2 (-\infty , \infty ) , \) and are orthonormal:
\[ \int_{-\infty}^{\infty} \psi_n (x) \,\psi_m (x) \,{\text d} x = \delta_{n,m} = \begin{cases} 1 , & \quad \mbox{if $n=m$}, \\ 0 , & \quad \mbox{if $n\ne m$}. \end{cases} \]
The Hermite functions are the eigenfunctions of the elliptic operator \( - \frac{{\text d}^2}{{\text d} x^2} + x^2 . \) These functions are solutions to the differential equation that involves a quantum mechanical, simple harmonic oscillator:
\[ \psi''_n (x) + \left( 2n +1-x^2 \right) \psi_n (x) =0, \qquad n=0,1,2,\ldots \]

Example: Quantum Mechanical Simple Harmonic Oscillator. The wave function \( \psi_n (x) \) are eigenfunctions of the differential operator: \( \displaystyle \left( \frac{{\text d}^2}{{\text d} x^2} - x^2 \right) \psi_n (x) = -(2n+1)\,\psi_n (x) . \)

We introduce a lowering (or annihilation) operator \( \hat{a} \) and a raising, or creation, operator \( \hat{a}^{\ast} \) according to the formulas:

\[ \hat{a} = \frac{1}{\sqrt{2}} \left( x + \frac{{\text d}}{{\text d}x} \right) \qquad \mbox{and} \qquad \hat{a}^{\ast} = \frac{1}{\sqrt{2}} \left( x - \frac{{\text d}}{{\text d}x} \right) . \]
that are adjoint to each other. They obey the commutation relations
\[ \left[ \hat{a} , \hat{a}^{\ast} \right] = \left[ \frac{{\text d}}{{\text d}x} , x \right] =1, \qquad \mbox{and} \qquad \left[ \hat{a} , \hat{a} \right] = 0 = \left[ \hat{a}^{\ast} , \hat{a}^{\ast} \right] . \]
Let us introduce the Hermitian number operator \( N= \hat{a}^{\ast} \hat{a} \) and the Hamiltonian
\[ {\cal H} = N + \frac{1}{2} = \hat{a}^{\ast} \hat{a} + \frac{1}{2} = \frac{1}{2} \left( \hat{a}^{\ast} \hat{a} + \hat{a} \hat{a}^{\ast} \right) . \]
Let \( |n\rangle \) be an eigenfunction of \( {\cal H} : \ {\cal H} \,|n\rangle = \lambda_n |n\rangle , \) whose eigenvalue λn is unknown at this point. Now we prove the key property that N has nonnegative integer eigenvalues:
\[ N \, |n\rangle = \left( \lambda_n - \frac{1}{2} \right) |n\rangle = n\, |n\rangle , \qquad n=0,1,2,\ldots , \]
that is, \( \lambda_n = n + 1/2 . \) Since \( \hat{a} |n\rangle \) is complex conjugate to \( \langle n | \hat{a}^{\ast} , \) the normalization integral \( \langle n | \hat{a}^{\ast} \, \hat{a} |n\rangle \ge 0 \) and is finite. From
\[ \left( \hat{a} |n\rangle \right)^{\ast} \hat{a} |n\rangle = \langle n | \hat{a}^{\ast} \, \hat{a} |n\rangle = \left( \lambda_n - \frac{1}{2} \right) \ge 0 , \]
we see that N has nonnegative eigenvalues.

Using commutation relations

\[ \left[ N, \hat{a}^{\ast} \right] = \hat{a}^{\ast} , \qquad \left[ N, \hat{a} \right] = -\hat{a} , \]
and \( \hat{a} \, \hat{a}^{\ast} = N + \left[ \hat{a}, \hat{a}^{\ast} \right] = N+1 , \) we find that
\begin{align*} N \left( \hat{a}^{\ast}\, |n\rangle \right) &= \hat{a}^{\ast} \left( \hat{a} \, \hat{a}^{\ast} \right) |n\rangle = \hat{a}^{\ast} \left( \left[ \hat{a}, \hat{a}^{\ast} \right] + N \right) |n\rangle \\ &= \hat{a}^{\ast} \left( N+1 \right) |n\rangle = \left( \lambda_n + \frac{1}{2} \right)= (n+1) \hat{a}^{\ast}\, |n\rangle = (n+1) \,\hat{a}^{\ast}\, |n\rangle , \\ N \left( \hat{a}\, |n\rangle \right) &= \left( \hat{a}, \hat{a}^{\ast} -1 \right) |n\rangle = \hat{a} \left( N-1 \right) |n\rangle = \left( n-1 \right) |n\rangle . \end{align*}
In other words, N acting on \( \hat{a} \,|n\rangle \) shows that \( \hat{a} \) has lowered the eigenvalue n by one unit; hence it is a lowering (or annihilation) operator. Similarly, N acting on \( \hat{a}^{\ast} \,|n\rangle \) shows that \( \hat{a}^{\ast} \) has raised the eigenvalue n corresponding to \( |n \rangle \)  by one unit, whence its name raising, or creation, operator. Therefore,
\[ \hat{a} |n\rangle \, \sim |n-1 \rangle , \qquad \hat{a}^{\ast} |n\rangle \,\sim\, |n+1 \rangle . \]
Applying \( \hat{a} \) repeatedly, we can reach all lower states, including the lowest, or ground, state \( |0\rangle , \) with eigenvalue λ0. We cannot step lower because \( \lambda_0 \ge 0. \) Hence, \( \hat{a}\,|0\rangle \equiv 0, \) suggesting we construct \( \psi_0 = | 0 \rangle \) from the (factored) first-order ODE
\[ \sqrt{2}\,\hat{a} \, \psi_0 = \left( \frac{\text d}{{\text d}x} +x \right) \psi_0 =0 . \]
Integrating
\[ \frac{\psi_0'}{\psi_0} = -x, \]
we obtain
\[ \ln \, \psi_0 (x) = -\frac{x^2}{2} , \]
where c0 is an integration constant. The solution
\[ \psi_0 (x) = c_0 e^{-x^2 /2} . \]
is a Gaussian that can be normalized, with \( c_0 = \pi^{-1/4} , \) using the error integral. Using ψ0, we find
\[ {\cal H} \,|0\rangle = \left( \hat{a}^{\ast}\hat{a} + \frac{1}{2} \right) |0\rangle = \frac{1}{2}\,|0\rangle . \]
So its energy eigenvalue is \( \lambda_0 = 1/2 , \) and its number eigenvalue is n = 0, confirming the notation \( |0\rangle . \) Applying \( \hat{a}^{\ast} \) repeatedly to \( \psi_0 = |0\rangle , \) all other eigenvalues are confirmed to be \( \lambda_n = n+1/2 . \) Using normalizations, we find
\[ \left\langle \,\vert \hat{a} \,\hat{a}^{\ast} \vert n \right\rangle = \left\langle \,\vert \hat{a}^{\ast} \,\hat{a} +1 \vert n \right\rangle = n+1 , \]
showing
\[ \sqrt{n+1}\,|n+1 \rangle = \hat{a}^{\ast} \vert n \rangle , \qquad \sqrt{n}\,|n-1\rangle = \hat{a} \vert n \rangle . \]
Thus, the excited-state wave functions, ψ1, ψ2, and so on, are generated by the raising operator
\[ |1 \rangle = \hat{a}^{\ast} \vert 0 \rangle = \frac{1}{\sqrt{2}} \left( x-\frac{\text d}{{\text d}x} \right) \psi_0 = \frac{x\sqrt{2}}{\pi^{1/4}} \, e^{-x^2 /2} , \]
yielding
\[ \psi_n (x) = N_n H_n (x) \,e^{-x^2 /2} , \qquad N_n = \pi^{-1/4} \left( 2^n n! \right)^{-1/2} . \]

 

Hermite Series


Hermite series:

\[ f (x) = \sum_{n\ge 0} c_n H_n (x) , \qquad c_n = \frac{1}{2^n n! \,\sqrt{\pi}}\, \int_{-\infty}^{\infty} f(x)\, e^{-x^2} \,H_n (x) \,{\text d} x, \qquad n=0,1,2,\ldots , \]
converges in \( L^p (-\infty , \infty ) \) for \( p \in \left( \frac{4}{3} , 4 \right) . \) The corresponding expansion for Dirac delta function is
\[ \delta (x) = \frac{1}{\sqrt{\pi}}\,e^{-(x^2 + a^2 )/2} \, \sum_{n\ge 0} \frac{1}{2^n \, n!}\, H_n (x)\, H_n (a) . \]

We demonstrate the Hermite expansion for the following polynomial \( p(x) = 1 + x + 3 x^2 + 7 x^3 . \) First, we determine its Hermite coefficients and then expand p(x) into Hermite series.

p[x_] = 1 + x + 3 x^2 + 7 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[1]= {5/2, 23/4, 3/4, 7/8}
coeffs.Table[HermiteH[n, x], {n, 0, 3}] // Expand
Out[2]= 1 + x + 3 x^2 + 7 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[#]], #] & @
PadRight[CoefficientList[p, x], Length[basis]]
expandPoly[1 + x + 3 x^2 + 7 x^3, HermiteH[Range[4] - 1, x], x]
Out[4]= {5/2, 23/4, 3/4, 7/8}

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]
Out[6]= {5/2, 23/4, 3/4, 7/8}

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]
Out[9]= 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
Out[10]= 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
Out[11]= 1/32 (312 H[0, x] + 380 H[1, x] + 200 H[2, x] + 60 H[3, x] + 10 H[4, x] + H[5, x])

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:

Herbert E. Salzer, A recurrence scheme for converting from one orthogonal expansion into another. Communications of the ACM CACM Homepage archive Volume 16 Issue 11, Nov. 1973 Pages 705-707 1. H.E. Salzer, Orthogonal Polynomials Arising in the Numerical Evaluation of Inverse Laplace Transforms, Math. Comp. 9, 164-177 (1955) 2. C.P. Jeffreson, E.-P. Chow, Least squares coefficients for a quadrature formula for Laplace transform inversion, J. Comp. Appl. Math. 4, 53-58 (1978)


poly = 1 - x + 2 x^2 - 5 x^4;
c = CoefficientList[poly, x];
n = Length[c] - 1; Remove[a];
a[0, 2] = c[[n - 1]] + c[[n + 1]]/2;
a[1, 2] = c[[n]]; a[2, 2] = c[[n + 1]]/2;
Do[
   a[0, k + 1] = c[[n - k]] + a[1, k]/2;
   a[1, k + 1] = a[0, k] + a[2, k]/2;
   Do[
      a[m, k + 1] = (a[m + 1, k] + a[m - 1, k])/2,
      {m, 2, k - 1}];
   a[k, k + 1] = a[k - 1, k]/2;
   a[k + 1, k + 1] = a[k, k]/2,
   {k, 2, n - 1}];

ccof = Table[a[m, n], {m, 0, n}]
   {1/8, -1, -3/2, 0, -5/8}

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}]
Out[13]= {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
Out[14]= True
CoefficientList[(1 + x)^5, x] // monomialToHermite
Out[15]= {39/4, 95/8, 25/4, 15/8, 5/16, 1/32}
%.HermiteH[Range[0, 5], x] == (1 + x)^5 // Expand
Out[16]= True

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
Out[18]= {3/2, 5/4, 1/4, 1/8}

Check:

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

Example: We expand \( \cos x \) into Hermite series. Using Mathematica, we find coefficients in its expansion:

\[ \cos (x) = e^{-1/4} \,\sum_{k\ge 0} \frac{(-1)^k}{2^{2k} \, (2k)!} \, H_{2k} (x) .
Then we plot its approximation with 10 terms:
Integrate[Cos[x]*HermiteH[6, x]*Exp[-x^2], {x, -Infinity, Infinity}]
hercos[m_] = Exp[-1/4]* Sum[(-1)^k/(2^(2*k) *Factorial[2*k])*HermiteH[2*k, x], {k, 0, m}]
Plot[{Cos[x], hercos[10]}, {x, -15, 15}, PlotStyle -> {{Thick, Blue}, {Thick, Red}}]

 

Example: Mathematica is helpful to confirm the expansion of a power function into Hermite series:

\[ x^n = \frac{n!}{2^n} \,\sum_{k= 0}^{\left\lfloor n/2 \right\rfloor} \frac{1}{k! \,(n-2k)!} \, H_{n-2k} (x) , \qquad n \in \mathbb{Z}_{+} . \]

 

 

Return to Mathematica page

Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions