I am trying to approximate the wavefunction of a particle in a delta potential $U(x) = -U_0 \delta(x)$ with $V_0 \gt 0$. I am using the following formula to calculate the wavefunction:

$\psi(x+\Delta x) = 2 \psi(x) - \psi(x-\Delta x) + (\Delta x)^2 (V(x) - \epsilon)\psi(x)$

This can be derived from the one dimensional, time-independent, one-dimensional schrödinger-equation $\psi''(x) = (U(x)-\epsilon)\psi(x).$

The following C++-code sippet shows the calculation:

    double phi_x_prev = 0;
    double phi_x_curr = 1e-100;
    double phi_x_next = 0;


    double stepSq = step * step;

    for (double x = -xmax + step; x <= xmax; x += step)
        phi_x_next =
            (2 * phi_x_curr - phi_x_prev)                           // O(1)
            + stepSq * (system->Potential(x) - eps)*phi_x_curr;     // O(step^2)


        phi_x_prev = phi_x_curr;
        phi_x_curr = phi_x_next;

The solution looks good for $x<0$, but for $x>0$ the wavefunction starts growing exponentially, which is not compatible with the boundary condition $\psi(x) \rightarrow \infty$ for $x \rightarrow \pm \infty$.

Screenshot of the plotted solution

What did I do wrong? Thanks in advance!


closed as off-topic by Norbert Schuch, Kyle Kanos, John Rennie, Jon Custer, Martin Jan 8 at 20:58

  • This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ What value of $\epsilon$ do you use? Delta-potential is known to have a unique bound state for the particular value of energy. $\endgroup$ – Gec Jan 6 at 14:18
  • 3
    $\begingroup$ I'm voting to close this question because questions regarding the writing and debugging are considered off-topic here. They should be modified to fit Stack Overflow or Computational Science or another of the other programming-related SE sites. $\endgroup$ – Kyle Kanos Jan 7 at 10:59

There are two generic sources of errors.

The first is the choice of guess energy for your solution. Quantization of energy is the result of applying boundary conditions, and it is by forcing the boundary conditions that one recovers discrete energy levels. If you do this in the reverse direction, i.e. you start with some energy and solve the Schrodinger equation, there is no guarantee that the boundary condition (here at $\infty$) will be satisfied unless you precisely pick one of the solution energies: for pretty much any other value of energy the solution will eventually diverge. The tolerance on $\epsilon$ is quite small: you can be off by less than $1\%$ and the solution will eventually diverge.

The second is numerical accuracy of the scheme. Even if you start with the exact and correct value of energy, the integration scheme will pick up some small round off errors as you move away from your initial conditions. This produces the same effect as picking a slightly wrong value of energy: the solution will behave for a while but eventually, for some large enough values of the argument, will start to exponentially diverge. There is no solution to this beyond using increasingly accurate integration scheme, either by increasing the order of the scheme, or by making the integration mesh finer.

As examples of these two types of error, consider the numerical solution to the Schrodinger equation for a harmonic potential. On the first figure, I solved the equation using the exact eigenvalue for $n=2$: the solution still eventually diverges; this illustrates the second type of error. On the bottom I took $n=2.001$, i.e. an error in guess eigenvalue of 1 part in 2000. Note how the solution starts to diverge for smaller $\vert x\vert$. This shows how a generic integration scheme is very sensitive to the guess eigenvalue; this illustrates the first source of errors. Both integrations were done using the same NDSolve function in Mathematica so there can be no issue of implementation.

enter image description here enter image description here

In your specific potential the control of numerical errors is critical because the derivative is discontinuous at the $\delta$. I suspect that some small numerical error occur when going across the $\delta$ and this messes up your solution when you continue to integrate on one side.


Not the answer you're looking for? Browse other questions tagged or ask your own question.