# Probability Function in a 1D Quantum Harmonic Oscillator [closed]

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)$.

## closed as off-topic by Jens, rhermans, MarcoB, rcollyer, eyorbleJul 29 '18 at 20:58

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Jens, rhermans, MarcoB, rcollyer, eyorble
If this question can be reworded to fit the rules in the help center, please edit the question.

• 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

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