# Numeric Solution vs Analytic Solution (1D Quantum Harmonic Oscillator)

We probably all know the analytic solution of the 1D Quantum Harmonic Oscillator is given as follows:

$$\Psi_n(x) = \dfrac{1}{\sqrt{2^n\,n!}}\,\left(\dfrac{m\,\omega}{\pi\,\hbar}\right)^{(1/4)}\,e^{\textstyle -m \,\omega\,x^2/(2\,\hbar)}\,H_n\left(\sqrt{\dfrac{m\,\omega}{\hbar}}\,x\right)$$

where $$\Psi$$ is the nth eigenfunction of the oscillator in position space.

You get these result by solving: $$\left(-\dfrac{\hbar}{2\,m}\,{\partial_x}^2 + \dfrac{1}{2}\,m\,\omega^2\,x^2\right)\Psi = E\,\Psi\quad \text{where} \quad E = (n + 1/2)\,\hbar\,\omega$$

Now to the numeric approach (Finite Difference Method)

By approximating the derivative $${\partial_x}^2\,\Psi_i \ = \dfrac{\Psi_{i+1} - 2\,\Psi_i + \Psi_{i-1}}{{\Delta x}^2}\quad \text{(with discrete points)}$$

the Schrodinger equation can be written:

\begin{align} &-\dfrac{\hbar}{2\,m}\left(\dfrac{\Psi_{i+1} - 2\,\Psi_i + \Psi_{i-1}}{{\Delta x}^2}\right) + V_i\,\Psi_i = E\,\Psi \\[12pt] &-\dfrac{\hbar}{2\,m\,{\Delta x}^2} \,\Psi_{i+1} + \left(\dfrac{\hbar}{m\,{\Delta x}^2} + V_i\right)\,\Psi_i - \dfrac{\hbar}{2\,m\,{\Delta x}^2}\,\Psi_{i-1} = E\,\Psi_i \end{align}

$$\text{where} \quad V_i = \dfrac{1}{2}\,m\,\omega^2{x_i}^2$$

Writing this as a Matrix-equation brings:

$$\left[\begin{array}{c} \frac{\hbar}{m\,{\Delta x}^2} + V_1 & -\frac{\hbar}{2\,m\,{\Delta x}^2} & 0 & 0 & 0 \\ -\frac{\hbar}{2\,m\,{\Delta x}^2} & \frac{\hbar}{m\,{\Delta x}^2} + V_2 & -\frac{\hbar}{2\,m\,{\Delta x}^2} & 0 & 0 \\ 0 & \ldots & \ldots & \ldots & 0 \\ 0 & 0& -\frac{\hbar}{2\,m\,{\Delta x}^2} &\frac{\hbar}{m\,{\Delta x}^2} + V_{N-2} &-\frac{\hbar}{2\,m\,{\Delta x}^2} \\ 0 & 0 & 0 & -\frac{\hbar}{2\,m\,{\Delta x}^2} & \frac{\hbar}{m\,{\Delta x}^2} + V_{N-1} \end{array}\right] \: \left[\begin{array}{c}\Psi_1 \\[6pt] \Psi_2 \\[6pt] \ldots \\[6pt] \Psi_{N-2} \\[6pt] \Psi_{N-1}\end{array}\right] \: = E\:\left[\begin{array}{c}\Psi_1 \\[6pt] \Psi_2 \\[6pt] \ldots \\[6pt] \Psi_{N-2} \\[6pt] \Psi_{N-1}\end{array}\right]$$

this automatically accounts for the boundary conditions $$\Psi_0 = 0$$ and $$\Psi_N = 0$$

Now all that has to be done is to solve for the Eigenfunctions $$\Psi$$ of this big Matrix numerically.

However what I found and what brings me to this question is that the analytic and numeric $$\Psi$$ are not identical

(This can be seen in the graphics I attached below)

If you were to zoom in close, the numeric solution indeed reflects the wave function-shape, but it's never at the same height as the analytic one.

I just realized this might be a better question for Programmers, but maybe someone of you is experienced.

• I find your boundary condition a bit strange so there might be an issue. However, it seems that your problem is mostly due to normalization. What is the normalization condition for the numeric solution? Commented Jun 4, 2023 at 15:02
• Yea, I didn't set a normalization myself. I just took the Eigenfunctions as they are. (I figured norm(v) = 1, but this ain't probably the same as in Quantum)
– Leon
Commented Jun 4, 2023 at 22:11

To elaborate on @Javi's answer, yes this is a normalization issue, but there's some detail that might be helpful. An eigenvector can be scaled by any complex number and is still an eigenvector (with the same eigenvalue). So having a wavevector which is correct except for scaling is a problem with a pretty straightforward solution.

As @Javi pointed out, normalization is necessary even in analytical solutions. In an analytical solution, you normalize according to an integral condition: $$\int_{-\infty}^\infty\Psi^*(x)\Psi(x)\,dx = 1$$ so that the probability density function has the right properties. If you can't do an integral, there are two ways of normalizing the wavevector, depending on what value you want to hold consistent: $$\sum_{i}\Psi_i^*\Psi_i = 1\quad\text{(Euclidean norm)}$$ or $$\frac{1}{\Delta x}\sum_{i}\Psi_i^*\Psi_i = 1\quad\text{(approximation of an integral)}$$

A computer calculating an eigenvalue will use the Euclidean norm by default, scaling the vector by $$C = \left(\sum_{i}\Psi_i^*\Psi_i\right)^{-\frac{1}{2}}.$$ Note that the sum will diverge as $$\Delta x$$ approaches zero, so it will scale your wavefunction smaller and smaller as you use a finer and finer approximation. In order to avoid this, it's better to use the other form, and scale by $$C = \left(\frac{1}{\Delta x}\sum_{i}\Psi_i^*\Psi_i\right)^{-\frac{1}{2}}.$$

This is because the computer interprets the meaning of the values in the approximation differently as you change $$\Delta x$$. Suppose you're approximating the interval from -10 to 10. If $$\Delta x = 0.1$$, then the computer interprets $$\Psi_0^*\Psi_0$$ as the probability of measuring the particle between -10 and -9.9. If $$\Delta x = 0.001$$, then $$\Psi_0^*\Psi_0$$ is interpreted as the probability of measuring the particle between -10 and -9.999. The computer scales the wavevector in accordance to this interpretation. When $$\Delta x = 0.1$$, the probabilities should be around 100 times larger than the probabilities when $$\Delta x = 0.001$$, and the wavevector components around 10 times larger.

But you don't want $$\Psi(-10)$$ to depend on the value of $$\Delta x$$ that you used, because you aren't looking for an actual discrete wavevector. You're looking for an approximation of a continuous wavevector. You want to know what is the infinitesimal probability density at a certain point. For this, you have to use the other normalization condition, where you scale the wavevector by $$C = \left(\frac{1}{\Delta x}\sum_{i}\Psi_i^*\Psi_i\right)^{-\frac{1}{2}}.$$

By contrast, if you were considering electron spins, you could use the Euclidean norm: $$\Psi = \frac{1}{2}\Psi_{up} + \frac{\sqrt{3}}{2}\Psi_{down}$$ $$\left(\frac{1}{2}\Psi_{up}\right)^*\left(\frac{1}{2}\Psi_{up}\right) + \left(\frac{\sqrt{3}}{2}\Psi_{up}\right)^*\left(\frac{\sqrt{3}}{2}\Psi_{up}\right) = \frac{1}{4} + \frac{3}{4} = 1,$$ but this is usually only the case when you're dealing with variables that will be measured discretely.

An excellent treatment of the change between discrete and continuous wavevectors can be found in R. Shankar Principles of Quantum Mechanics section 1.10.

• Thank you very much for these detailed explanation and for the great book source. I was seeking such a good guidance for a long time. Most likely I will find my remaining queries right up there, but one thing to add on already are the major differences that are built into the boundary conditions both of Analytic (0 at inf) and Numeric (0 at endpoints) which will probably cause them to never line up
– Leon
Commented Jun 6, 2023 at 17:02

I think that the numerical scheme is incomplete because it lacks the proper normalization of the Eigenfunction.

Say that $$\Psi_n$$ for $$n\in[0,N]$$ is the numerical solution that you find applying finite differences and the boundary conditions. Then the proper Eigenfunction should be $$\Psi'_n=C\Psi_n$$ with $$\frac{1}{C}=\left(\sum_{n=0}^N\Psi_n (\Psi_n)^*\right)^{\frac{1}{2}}.$$

Note that the normalization step is also usually required when you solve the Schrödinger equation analytically.

Also, the errors introduced by the discretization will decrease as $$N$$ becomes bigger.

• This is correct. Eigenvectors are not unique, they can be scaled by any complex number. When the eigenvector is determined numerically, the scale isn't chosen with quantum normalization in mind, so it must be rescaled. Commented Jun 4, 2023 at 16:34
• Just for Completeness: are you saying I just need to multiply my numeric state by a scalar and eventually it will cover up with the analytic one? Well this is not the case, both also have different stretching. However I assume it is to difficult to derive where this is coming from since it's all heavily based on the interior of the programming code.
– Leon
Commented Jun 4, 2023 at 22:17
• @Leon The spreading is due to the finite difference method being an approximation. If you use a smaller $\Delta x$ or a wider interval, the errors will decrease. Commented Jun 4, 2023 at 23:17