# How do you randomly draw samples from the probability density function of the quantum harmonic oscillator in MATLAB?

The Quantum Harmonic Oscillator in the ground state is specified by the following Gaussian PDF in two dimensions:

$$p(x,y)= \frac{M \omega_x}{\pi h}\sqrt{ \frac{\omega_y}{\omega_x}} e^{-\frac{M}{h}(\omega_xx^2+\omega_yy^2)}$$

where M = mass of particle, h = plank's constant, $$\omega_x=\frac{h}{2M\sigma_x^2}$$ and $$\omega_y=\frac{h}{2M\sigma_y^2}$$ are related to the variances $$\sigma_x$$ and $$\sigma_y$$in the x and y direction respectively.

How do I draw a 1000 samples from a normal distribution like this in MATLAB?

I do not think I can simply use randn(1000,1) because I need to somehow take into consideration the constants $$M,h,\pi, \omega_x, \omega_y$$ when drawing samples. Is there way to accomplish this in MATLAB assuming all these constants are predefined?

$$p(x,y)=p_x(x)p_y(y).$$
This means that your $$x$$ and $$y$$ random variables are independent and can be generated from their respective PDFs $$p_x$$ and $$p_y$$. After you get expressions for (normalized!) $$p_x$$ and $$p_y$$ you can simply use inverse transform sampling to generate $$x$$ and $$y$$ respectively.
• @Ruslan- Thanks for your insight. I just factored out $p(x,y)$ like you said above and then plugged in $\omega_x$ into $p_x(x)$ and $\omega_y$ into $p_y(y)$ and wind up getting the normal distribution with zero mean and variances $\sigma_x$ and $\sigma_y$ respectively. Do you know what the typical values of the variance of a QHO is? Does it make sense to set the variance in the x direction to 1 and the variance in the y direction to 10? is that behavior reflective of reality? – user1068636 May 19 at 14:28