numerical methods

Calculating the Hermite functions

The Hermite functions appear as the solutions of the quantum mechanical harmonic oscillator. But they have applications in many other fields and applications, e.g., pseudospectral methods. The Hermite functions $h_n(x)$ are defined as
\begin{equation}
\label{eq:h}
h_n(x) = \frac{1}{\sqrt{\sqrt{\pi} 2^n n!}} \mathrm{e}^{-x^2/2} H_n(x) \,,
\end{equation} where $H_n(x)$ denotes the $n$th Hermite polynomial defined via the recurrence relation
\begin{equation}
H_{n}(x) = 2xH_{n-1}(x)-2(n-1)H_{n-2}(x)
\end{equation} with the initial values $H_0(x)=1$ and $H_1(x)=2x$. Although the Hermite functions are quite well behaved, they are rather difficult to calculate especially for large $n$ and/or large $x$.

Calculating the Hermite functions via the definition \eqref{eq:h} fails easily due to numerical overflow or underflow that is caused by the rapid growth or decrease of the individual factors $\frac{1}{\sqrt{\sqrt{\pi} 2^n n!}} \mathrm{e}^{-x^2/2}$ and $H_n(x)$. A partial solution to this problem is to introduce the modified Hermite polynomials defined via the recurrence relation
\begin{equation}
\label{eq:h2}
\tilde H_{n}(x) = \sqrt{\frac{2}{n}}x\tilde H_{n-1}(x)-2\sqrt{\frac{n-1}{n}}\tilde H_{n-2}(x)
\end{equation} with the initial values $\tilde H_0(x)=1/\sqrt[4]{\pi}$ and $\tilde H_1(x)=\sqrt{2}\,x/\sqrt[4]{\pi}$. With these polynomials the Hermite functions become
\begin{equation}
h_n(x) = \mathrm{e}^{-x^2/2} \tilde H_n(x) \,.
\end{equation} Because the modified Hermite polynomials $\tilde H_{n}(x)$ grow much more slowly than the standard Hermite polynomials we got rid of the normalizing factor $\frac{1}{\sqrt{\sqrt{\pi} 2^n n!}}$.

But this does not solve all problems of underflow and overflow. For $x>38$ the factor $\mathrm{e}^{-x^2/2}$ is smaller than any number that can be represented as a double precision floating point number. A very robust remedy to numerical underflow was sketched in BIT Numerical Mathematics, Vol. 49, pp 281-295. The basic idea is to keep the magnitude of the modified Hermite polynomials on the order of one during their calculation via the recurrence relation \eqref{eq:h2} by introducing a suitable normalizing factor. During the calculation one has to keep track of the sum of the logarithms of these normalizing factors. Then the factor $\mathrm{e}^{-x^2/2}$ has to be modified by this sum when the value of the Hermite function is calculated. A Python implementation of this algorithm is shown below. This algorithm allows to calculate the values of high-order Hermite functions even for quite large arguments.

from numpy import *
from pylab import *
 
def h(n, x):
    if n==0:
        return ones_like(x)*pi**(-0.25)*exp(-x**2/2)
    if n==1:
        return sqrt(2.)*x*norm*pi**(-0.25)*exp(-x**2/2)
    h_i_2=ones_like(x)*pi**(-0.25)
    h_i_1=sqrt(2.)*x*pi**(-0.25)
    sum_log_scale=zeros_like(x)
    for i in range(2, n+1):
        h_i=sqrt(2./i)*x*h_i_1-sqrt((i-1.)/i)*h_i_2
        h_i_2, h_i_1=h_i_1, h_i
        log_scale=log(abs(h_i)).round()
        scale=exp(-log_scale)
        h_i=h_i*scale
        h_i_1=h_i_1*scale
        h_i_2=h_i_2*scale
        sum_log_scale+=log_scale
    return h_i*exp(-x**2/2+sum_log_scale)

The Hermite function $h_{800}(x)$. A direct calculation via the definition \eqref{eq:h2} of the Hermite function in terms of the modified Hermite polynomials $\tilde H_n(x)$ fails for $x>38$ due to numerical underflow as shown in the upper part. The shown Python code, which avoids underflow, has been utilized to produce the lower part.
The Hermite function $h_{800}(x)$. A direct calculation via the definition \eqref{eq:h2} of the Hermite function in terms of the modified Hermite polynomials $\tilde H_n(x)$ fails for $x>39$ due to numerical underflow as shown in the upper part. The shown Python code, which avoids underflow, has been utilized to produce the lower part.

I would like to thank Randolf Beerwerth for drawing my attention to this problem.

Leave a Reply

Your email address will not be published. Required fields are marked *