Determine numerical infinity for Schrodinger equation $−\psi''(x) + x^ 2 \psi(x) = E\psi(x)$

Consider the following Schrodinger equation for the harmonic oscillator with real $x$: $$−ψ''(x) + x^ 2 ψ(x) = Eψ(x).$$

I solve the last equation using shooting method and implicit Runge-Kutta integration. For the eigenvalue $E=1$, the corresponding eigenfunction looks like: where I set $x=4$ as $+∞$ and $x=-4$ as $-∞$ . From the Fig, we see that the eigenfunction exponentially decay at $x=\pm 4$ .

However, when eigenvalue $E=17$, the corresponding eigenfunction is too complicated so that it won't exponentially decays at $x=\pm 4$ as shown on the following Fig. In fact, the eigenfunction decays at $x=\pm 6$ as shown on the following Fig. My Question

The numerical infinity should be large (like $x=6$) for large eigenvalue $E$, and small (like $x=4$) for small eigenvalue. If independent from the eigenvalue, I uniformly set numerical infinity as a very large number, it would waste computational resource. So, is there any formula to determine the exact numerical infinity based on different inputted eigenvalue please?

• What do you call "decaying"? What is the right threshold for it? 1/10, 1/100, 1/1000... The exact numerical infinity is at inifinity, because that is the domain of the problem. On the other hand, how do you know the eigenvalues beforehand? Apr 18 '16 at 17:28
• The eigenvalue can be analytically calculated with results: $E = 1, 3, 5, 7...$. I am sorry that due to low reputation, I only upload two figures. Now from the 3rd Fig, the eigenfunction is decayed in the sense that the derivative of the eigenfunction with respect to $x$ at the numerical infinity is almost zero.
– pig
Apr 18 '16 at 20:54
• The eigenvectors can be analytically calculated as well. Apr 18 '16 at 20:55
• Please don't cross post. Apr 18 '16 at 21:59
• Crossposting is discouraged. See this discussion. I am not a mathematician, but a lot of people in SciComp are... so I think that your claim is quite incorrect. Apr 18 '16 at 22:15

You look at the question the wrong way. First, let us consider the exact problem and its solution: $$-\psi''(x) + x^2 \psi(x) = E\psi(x),$$ with boundary values $\psi(\pm\infty)=0$. Now you want to consider truncating your domain somewhere, i.e., you want to consider a different problem that reads $$-\tilde\psi''(x) + x^2 \tilde\psi(x) = E\tilde\psi(x),$$ with boundary values $\tilde\psi(\pm x_\infty)=0$ where $x_\infty$ is a finite termination point of your domain.

Clearly, these two problems will have different answers, $\psi$ and $\tilde\psi$. The question you are asking is, in essence, how large you have to choose $x_\infty$. But this cannot be answered without also saying how small you want the difference $\psi-\tilde\psi$ to be!

If you choose $x_\infty$ small, then you can solve the problem of course, but $\psi-\tilde\psi$ will be large. On the other hand, if you want the difference to be small, you will likely have to choose $x_\infty$ to be large. In general, if you require that $\|\psi-\tilde\psi\|\le\varepsilon$, then there will be some value $\bar x(\varepsilon)$ so that any $x_\infty\ge \bar x(\varepsilon)$ will be a valid choice of the domain size.

• If I understand your answer correctly, you suggest that I need numerically solve two equations simultaneously - the original equation and the complex conjugate of the original equation - afterwards, I compare the difference $\psi(x) - \bar \psi(x)$ from the two solutions at the chosen numerical infinity $x_\infty$. If the difference is not vanished, then I have to increase the chosen value for $x_\infty$ and numerically solve the two equations simultaneously again. Then do the comparison again. I think that your method may work, but numerically it is expensive. Anyway, I give you +1.
– pig
Apr 23 '16 at 6:25
• No, no, you misunderstand. $\bar\psi$ isn't the complex conjugate (I should have chosen a different notation, maybe $\tilde\psi$). It's just the solution of a different problem, on a finite domain, whereas $\psi$ is the solution on the infinite domain. Of course, numerically, you can't solve for $\psi$ because -- well -- the domain is infinite. Apr 23 '16 at 13:01
• I've edited the answer to use $\tilde\psi$ instead of $\bar\psi$. Apr 23 '16 at 13:03
• Thank you for your clarification! Looks like I need to do the following things: After numerically calculating $\tilde \psi(x)$, then due to exponentially decay (i.e. $\psi'(\infty) = 0$), $\tilde \psi'(x_\infty)$ must be vanished as well. Numerically, I can set $\tilde \psi'(x_\infty) < \epsilon$. If the last inequality is not satisfied, then I must increase $x_\infty$, then solve the ODE again and find the new solution $\tilde \psi(x)$ until the last inequality is satisfied. All looks well, except that I may never find $\tilde \psi(x)$ at first hand due to incorrect selection of $x_\infty$.
– pig
Apr 23 '16 at 22:19
• The Runge-Kutta integration may not guarantee to converge if $x_\infty$ is not appropriately chosen. If it's not converged, then $\tilde \psi(x)$ can not be obtained, consequently there is no way to do the following test $\tilde \psi'(x_\infty) < \epsilon$. But thank you for your answer!
– pig
Apr 23 '16 at 22:26

You should maybe look at the analytical form of these solutions (here's the wikipedia page), this would probably help you to derive a relation between the energy level $n$ and the typical length of your state (I would look at the $H_n$ functions given in the link above).

• It would be much better to determine the length of my state without knowing the analytical form. However, I give you +1, since your answer implies that I should look at asymptotic solution for the eigenfunctions. The asymptotic form must provide information where $\psi'(x)$ and $\psi(x)$ vanishes on the $x$.
– pig
Apr 18 '16 at 22:02