Generate high n quantum harmonic oscillator states numerically

Question

How can I generate the higher $n$ quantum harmonic oscillator wavefunction (in position space) numerically? Here, higher means around $n=500$, or say $n=2000$, where $n$ is the $n$th oscillator wavefunction.

---

The position-space wavefunction of the $n$th state involves Hermite polynomials of order n (see Griffith's book for the detailed form of the solution). I can generate wavefunction up to $n\approx 160$ only using 'gsl' library in C because the library doesn't offer Hermite polynomials of order higher than $\approx 170$.

To get rid of this limitation, I found that the higher $n$ wavefunctions can be achieved by using 'raising operator', and then repeatedly acting it on the previous wavefunction to get the next wavefunction, and so on.

Raising operator: $a^{\dagger} = \frac{1}{\sqrt{2}} \left(x - \frac{d}{dx}\right)$. (Consider, constants in the expression to be 1.)

In the expression, the derivative term, i.e., $d/dx$ needs to be evaluated. It can be solved either by the finite-difference method or by using Fourier transform to get 1st derivative (which is more accurate than the finite-difference method).

I wrote a Matlab program for $a^{\dagger}$ (shown below) using the 2nd method, i.e., Fourier transform. After evaluation of states greater than $n=10$, the program is becoming unstable and becoming even worse for higher $n$ (instability plot is shown below in the link). The same issue had appeared when I used the finite-difference method to evaluate the derivative (code is not attached).

Solution to the issue or an alternate approach to the problem will be helpful.

UPDATE: I'm also interested to know the reason for the development of instability while using 'fft' in the code.

sigma = 1.0;
xmin = -10.0;
xmax = 10.0;
npts = 512;
nstates = 14;

dx = (xmax-xmin)/npts;
x = xmin + dx*(0:npts-1);

% -- initial state/wavefunction
psi_init = exp(-0.5*x.^2/sigma^2)*(pi*sigma^2)^(-0.25);

psi = zeros(nstates,npts);             % -- list to store oscillator states
psi(1,:) = psi_init;
for nn=2:nstates  
     psi(nn,:) = raising_psi(psi(nn-1,:),xmin,xmax,npts,sigma); 
end

function adag_fn = raising_psi(previous_fn,a,b,n,sigma)     % -- raising_operator
 dx = (b-a)/n;
 x = a + dx*(0:n-1);

 % -- going into fourier space using 'fft' library 
 fwd_fft = fft(previous_fn);
 k = (2*pi/(b-a))*[0:n/2-1,0,-n/2+1:-1];
 dfk_dx = 1i*k.*fwd_fft;                     % -- 1st derivative in fourier space
 df_dx = ifft(dfk_dx);                       % -- back into position space

 % -- a^dagger acting on the previous state
 adag_fn = (x.*previous_fn/sigma - sigma*real(df_dx))/sqrt(2);
 norm_fn = adag_fn*transpose(adag_fn);
 adag_fn = adag_fn/sqrt(norm_fn);            % -- normalization
end

$n=10$:

$n=11$:

Answer 1

I'm posting my comment as an answer, partly to guard against link rot, and partly because I noticed a couple of minor errors.

I found this at Heiko Bauke's blog, which also refers to a paper by B Bunck BIT Numerical Mathematics, 49, 281 (2009) where more details are given (and yet another method involving evaluation of contour integrals is proposed). So the credit really belongs to them.

The aim is to calculate the normalized Hermite functions $$ h_n(x) = c_n H_n(x) \exp(-x^2/2), \quad c_n = \frac{1}{\pi^{1/4}}\,\frac{1}{\sqrt{2^n \, n!}} $$ and the problem is that at high $n$, the Hermite polynomials $H_n(x)$ become very large, and the normalizing factors $c_n$ become very small, as well as the function $\exp(-x^2/2)$ when $x$ is large. The starting point is the recurrence relation $$ H_n(x) = 2x H_{n-1}(x) -2(n-1) H_{n-2}(x) $$ where $H_0(x)=1$ and $H_1(x)=2x$. The corresponding recurrence relation for the functions $h_n(x)$ is $$ h_n(x) = \sqrt{\frac{2}{n}} \, x \, h_{n-1}(x) - \sqrt{\frac{n-1}{n}} h_{n-2}(x) $$ and the problem of large values partly goes away. Here I think there is a typo in eqn (3) of the blog post (a spurious factor $2$ in the second term) but the typo is not repeated in the code. Then, in the paper, a method is described for handling the evaluation of the functions at large $x$. Basically a scaling factor is progressively introduced at higher $n$ to keep the scaled function within reasonable bounds. This method is mentioned on the blog, and is implemented in the code provided.

To guard against link rot, I'm copying the code here from Bauke's blog, but I emphasize that credit really belongs to Bunck and Bauke. I have made a couple of changes. Firstly, there was a typo in the statement evaluating $h_1$ (an undefined variable norm) which I've corrected. Secondly, I've inserted a guard against taking the log of zero, which can happen when one of the $x$ points coincides exactly with a zero of $h_n(x)$: the most obvious instance is $x=0$ for odd $n$.

def h(n, x):
    import numpy as np
    if n==0:
        return np.ones_like(x)*np.pi**(-0.25)*np.exp(-x**2/2)
    if n==1:
        return np.sqrt(2.)*x*np.pi**(-0.25)*np.exp(-x**2/2)
    h_i_2=np.ones_like(x)*np.pi**(-0.25)
    h_i_1=np.sqrt(2.)*x*np.pi**(-0.25)
    sum_log_scale=np.zeros_like(x)
    for i in range(2, n+1):
        h_i=np.sqrt(2./i)*x*h_i_1-np.sqrt((i-1.)/i)*h_i_2
        h_i_2, h_i_1=h_i_1, h_i
        abs_h_i=np.abs(h_i)
        log_scale=np.log(abs_h_i,out=np.zeros_like(abs_h_i),where=abs_h_i>1).round()
        scale=np.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*np.exp(-x**2/2+sum_log_scale)

Any other changes are cosmetic. I hope I haven't introduced any new errors. You should check this for errors if you use it or convert it to Matlab. This approach using the recursion relation is most suitable when you want to calculate the $h_n(x)$ for many values of $n$ up to some maximum, in which case this function can be re-written to return all of them. I've done some preliminary tests, here for instance is the result for $n=800$ at high $x$, pretty much copying what was on Bauke's blog.

---

EDIT following OP edit to, and comments on, question.

The Hermite function recursion relations are \begin{align*} h_n'(x) &= \sqrt{\frac{n}{2}} h_{n-1}(x) - \sqrt{\frac{n+1}{2}} h_{n+1}(x) \tag{1} \\ xh_n(x) &= \sqrt{\frac{n}{2}} h_{n-1}(x) + \sqrt{\frac{n+1}{2}} h_{n+1}(x) \tag{2} \end{align*} Subtracting (2) - (1) gives the equation you are using to calculate $h_{n+1}(x)$, with the "raising operator" appearing on the left, and $h_n'(x)$ approximated using FFT. Eqn (2) is used in the algorithm given above in my answer. I have used this to generate $h_n(x)$ functions as a reference, and hence obtained accurate derivatives through eqn (1). Comparing with your Fourier evaluation of $h_n'(x)$, from the same accurate functions $h_n(x)$, shows that Gibbs-type oscillations do indeed start appearing, due to the discontinuities at the edges. Iterating your algorithm would be expected to amplify these effects at each stage. I cannot be sure that this is the complete answer to the problems with the FFT approach, but that's as far as I'm willing to go. If you are continuing the analysis, I suggest that you look at the accuracy of each derivative step in the way I just described, rather than simply examining the output of the whole algorithm, which will show the accumulated errors.

Answer 2

I did some more investigating about the accuracy of the standard recurrence relation method vs. the newer Bunck algorithm. It seems that in fact the Bunck algorithm is generally more accurate for all $n$ and $x$. I made a python script using the mpmath arbitrary precision library to calculate the Hermite functions to lots of decimal places, and then compared that result with the GSL recurrence relation approach and the Bunck method as implemented in LonelyProf's answer.

Below is a plot for "small" $x$ in $[-1,1]$. The blue dots are mpmath-GSL (i.e. the recurrence relation result). The green dots are mpmath-Bunck (labeled NEW). Here I used $n = 10,000$. The red lines show $\pm \epsilon_{mach} = \pm 2.22e-16$. For this range of $x$, there isn't much difference in accuracy, and both algorithms give an absolute error close to the machine epsilon. The Bunck algorithm gives a slightly smaller rms of $2.50e-16$ compared with $2.59e-16$ for GSL.

Now to explore larger $x$ ranges, I have the below figure:

where the horizontal axis is plotted as $1/x$, corresponding to $x \in [-\infty,-1]$ and $x \in [1,\infty]$. You can see that near $1/x \approx 0$, the error in the GSL algorithm based on recurrences blows up substantially. This occurs near $1/x \approx 0.025$ which corresponds to $x = 40$, which as the previous authors note, is where $\exp{-x^2/2}$ underflows in double precision. Reducing the $y$ scale in this figure gives us:

which shows that both methods give similar accuracy until the "large x" region in the middle.

I made similar plots for a wide range of different $n$ values and found basically the same behavior overall. My conclusion is that the new Bunck algorithm is as accurate as recurrences for small $x$, and more accurate for large $x$. If Bunck's claim is true that his method is also faster than recurrences, which he says in his paper (I have not done any benchmarks), then I would conclude his method is superior to the recurrence method in general.