HarmonicOscillatorWaveFunction fails for small Fock states

[joshcombes]

Bug Description

Using the HarmonicOscillatorWaveFunction class one can compute and plot the position wave function.

However, for relatively low Fock states the routine used by Qutip begins to fail. This makes it hard to plot e.g. the GKP position wavefunction which requires a lot of Fock states to converge.

The bug starts to appear for Fock state n=43. Below I have only plotted n = 47 because that's the first time it becomes very noticeable.

The (crappy) code snippet I provided seems to fail around n = 151.

Code to Reproduce the Bug

import numpy as np
import matplotlib.pyplot as plt
import qutip as qt
from qutip.distributions import HarmonicOscillatorWaveFunction
from numpy.polynomial.hermite import hermval
import math


def harmonic_oscillator_wavefunction(psi, xvecs, omega0=1.0):
    n_coeff = []
    for n in range(psi.shape[0]):
        n_coeff.append(pow(omega0 / np.pi, 0.25) / \
                np.sqrt(2.0 ** n * math.factorial(n)) 
                *  psi.data[n, 0])
    return hermval(xvecs[0], n_coeff)* np.exp(-xvecs[0] ** 2 / 2.0)

xmax = 12
steps = 3000
x = np.linspace(-xmax,xmax,num=steps)

# choose a fock state
n = 47
state = qt.basis(n+1,n)

wav_fn = HarmonicOscillatorWaveFunction(state, omega=1, extent=[-xmax,xmax], steps=steps)
y_qutip = wav_fn.data
y_good = harmonic_oscillator_wavefunction(state,[x])

plt.plot(x,y_qutip,label="qutip")
plt.plot(x,y_good,':',label="homebrew fn")
plt.plot(x,y_qutip-y_good,label= "qutip - homebrew")
plt.legend()
plt.title('Harmonic Oscillator Wave Function n = '+str(n))
plt.xlabel('position')

Code Output

No response

Expected Behaviour

Your Environment

QuTiP Version:      4.7.2
Numpy Version:      1.25.1
Scipy Version:      1.11.1
Cython Version:     None
Matplotlib Version: 3.7.2
Python Version:     3.11.4
Number of CPUs:     8
BLAS Info:          OPENBLAS
OPENMP Installed:   False
INTEL MKL Ext:      False
Platform Info:      Darwin (arm64)

Additional Context

No response

[Ericgig]

Thank you for reporting.

Everything in qutip.distribution is experimental: not tested and little documentation.
We probably won't have time to look at it soon, but we will keep the issue open until we do.