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

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

$n=10$: wavefunction for n=10

$n=11$: wavefunction for n=11

  • $\begingroup$ Why not using finite differences for the whole thing? $\endgroup$
    – nicoguaro
    Jan 17, 2019 at 17:21
  • 1
    $\begingroup$ I'd suggest using an analytic expression for the derivatives and evaluating it if you need a high-order function. You can easily get different expressions for $\partial_n e^{-x^2}$ from, e.g., WolframAlpha which are ostensibly simpler to evaluate in a stable way. $\endgroup$
    – Nox
    Jan 17, 2019 at 20:16
  • 2
    $\begingroup$ Have you come across the solution described here numbercrunch.de/blog/2014/08/calculating-the-hermite-functions ? They provide a Python implementation, and plot a sample result for $n=800$. $\endgroup$
    – user28077
    Jan 17, 2019 at 21:36
  • $\begingroup$ @nicoguaro I think it'll create trouble for n~500. That's why I'm following fourier transform's way. Thanks. $\endgroup$
    – ph-jd
    Jan 18, 2019 at 12:27
  • $\begingroup$ @Nox Actually I thought about that, but unable to find a way to run and store the values of the wavefunction in mathematica using the approach I presented. I have to look if there is a way to make matlab and mathematica talk and work together in a program. Thanks for the help. Further suggestions are welcome. $\endgroup$
    – ph-jd
    Jan 18, 2019 at 12:28

2 Answers 2


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)
    for i in range(2, n+1):
        h_i_2, h_i_1=h_i_1, h_i
    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.

enter image description here

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.

  • $\begingroup$ This is interesting, but it appears that this algorithm is not as accurate as the standard recurrence relation method for some values of $x$ and $n$. For example, using the inputs $n = 100028$, $x = 0.75$, the mpmath arbitrary precision module reports the value is: $-0.02903467369856961147236598085605007873807791777299928240653$. However your code above reports: $-0.02903467369856992$ for which the last 2 digits (significant digits 14 and 15) are incorrect. $\endgroup$
    – vibe
    Jan 28, 2019 at 6:23
  • $\begingroup$ The paper by Bunck claims his method is accurate to within the machine epsilon, so I wonder if there is a way to modify your code to achieve this? For reference, the GSL implementation (based on recurrence relations) for these input arguments does produce the correct 14th and 15th digits. $\endgroup$
    – vibe
    Jan 28, 2019 at 6:23
  • $\begingroup$ The OP's question was primarily about a stable way of generating these functions for $n$ up to $\sim 2000$, and about the instability in the FFT-based algorithm. I tried to address these points in my answer. Bearing in mind that I am not the originator of the code in my answer, I'm not personally inclined to pursue this aspect of high accuracy, although I am sure it is possible to do so. $\endgroup$
    – user28077
    Jan 28, 2019 at 12:10

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]$. enter image description here 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: enter image description here

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: enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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