Probability Function in a 1D Quantum Harmonic Oscillator

Question

Let me start off by saying that I am a complete newbie to Mathematica, so I don't really know what I'm doing.

For my assignment I have to find the numerical probability of a particle in a harmonic oscillator potential in between quantum numbers $n=0$ and $n=5$. For simplicity's sake, I am only trying to find the probability where $n=0$.

The wave function of n=0 in a harmonic oscillator is:

$\Psi(x)$ = $N_0H_0 e^{(-x^2/2)}$

So the probability of finding a particle with the given wave function is:

$\int \Psi^2(x) dx$

The classically bound region is defined as $y= \frac 12{kx^2}$

ClearAll["Global`*"];
norm[n_] := (1/(Sqrt[π] 2^n n!))^(1/2)

u[x_, n_] := 1/Sqrt[a] norm[n] HermiteH[n, x/a] exp^[-(a*x^2/2)]

b[x_]  := 0.5 kx^2

φ1[x_, n_] :=  1/Sqrt[2 π] NIntegrate[u[x, n]^2, {x, -1, 1}, {n, 0, 5}]

However, this does not result in any output. I am wondering how you could format this to result in a valid output, and how to get the numerical probability within the bounds of $y= \frac 12{kx^2}$ and $E_v = \hbar \omega (v + 1)$.

This problem just happens to have been solved in the Neat examples under Filling. Type ?Filling, click on the >> and scroll to the bottom. — David G. Stork, Jul 16 '18 at 19:22
Thanks! How would I be able to find the integral of the area underneath the squared functions in between the parabola? — Dylan Hendrickson, Jul 16 '18 at 20:09
There are several syntax errors in your code. Replace exp by Exp, separate k and x^2, define a. The definition φ1[x_, n_] := doesn't make sense. Use φ1[ n_] := instead. Then evaluate φ1[1]. The function is also missing a calculation of the classical turning points as integration limits. Related: Find eigen energies of time-independent Schrödinger equation. — Jens, Jul 16 '18 at 20:14
Ok if I have my function as:
`f[n_, x_] := Abs[((1/Pi)^(1/4) HermiteH[n, x])/(E^(x^2/2) Sqrt[2^n n!])]^2

xtp := Sqrt[(2*n + 1)/a]

a := 1

NIntegrate[f[0, x], {x, -t, t}]` <\br> — Dylan Hendrickson, Jul 16 '18 at 20:23
I don't know if your above function definition will (eventually) get you want you want but there are three things: (1) NIntegrate does not do symbolic integration, i.e., you need to assign values to t, (2) You can use Integrate, and (3) you should put such modifications in the question rather than the comments. That's what the edit button is for. — JimB, Jul 16 '18 at 20:32
You should read about the difference between Set and SetDelayed. You won't get output from SetDelayed (:=), but from =. Basically, := defines functions that will be evaluated when you invoke their name again later (usually with arguments). — Jens, Jul 16 '18 at 21:11

score 2 · Answer 1 · answered Jul 17 '18 at 00:04

If you use NIntegrate, you need to assign all the values. Correcting your syntax:

ClearAll["Global`*"];
norm[n_] := (1/(Sqrt[Pi] 2^n n!))^(1/2)

u[x_, n_] := 1/Sqrt[a] norm[n] HermiteH[n, x/a] Exp[-(a*x^2/2)]

Assign a value to a, which I don't see in your code:

a = 1

You are integrating over x, but you want a table for n, so tell that to MMa.

Table[1/Sqrt[2 Pi] NIntegrate[u[x, n]^2, {x, -1, 1}], {n, 0, 5}]
(*{0.336189, 0.170585, 0.0877831, 0.115384, 0.0808829, 0.0739828}*)

These values give you the probabilities for finding the particle for x from -1 to 1, but that is not the range of x values you want if you want the limits of the parabolic potential. To find those values, you need to set E == V, the potential and energy, and solve for x from

(1/2)*k*x^2 == (n + 1/2)*h*nu

which requires assigning all the values for that equation also.
In this case you can also use Integrate requiring no prior value assignments.

a =.

Table[1/Sqrt[2 Pi] Integrate[u[x, n]^2, {x, -xm, xm}], {n, 0, 5}] // Simplify

Most of the values get too long to print here, but for example for n = 0 we get:

Erf[Sqrt[a]*xm]/(Sqrt[2*Pi]*a^(3/2))

where xm are the x limits found from setting E = V

Probability Function in a 1D Quantum Harmonic Oscillator

1 Answers1