Numerical approximation of the wavefunction in a delta-potential [closed]

Question

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;

    m_value.push_back(phi_x_curr);

    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)

        m_value.push_back(phi_x_next);

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

What did I do wrong? Thanks in advance!

Answer 1

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.

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.