Solve the stationary Schrödinger equation with the finite difference method [closed]

Question

I am trying to solve the stationary Schrödinger equation for a double-level well potential through the finite difference method.

Here is the shape of the potential I would like to solve it for

where the second walls are infinite.

The Schrödinger equation for this problem is the following

$$ -\frac{\hbar^2}{2m} \partial_{xx}\psi(x) = [E-V(x)]\psi(x)$$

which can be turned into

$$ \partial_{xx}\psi(x) = \frac{2m}{\hbar^2}[V(x)-E]\psi(x) $$

we can discretise the wave function and approximate the second derivative of $\psi_i$ as

$$\partial_{xx}\psi_i \sim \frac{\psi_{i-1} - 2\psi_i + \psi_{i+1}}{h^2}$$

Then the Schrödinger equation becomes

$$ \frac{\psi_{i-1} - 2\psi_i + \psi_{i+1}}{h^2} = \frac{2m}{\hbar^2}[V(x_i)-E]\psi(x) $$

which can be put as an eigenvalue problem as this

$$ \begin{pmatrix} -2 & 1 & 0 & ... & 0\\ 1 & -2 & 1 & ... & 0\\ 0 & .&.&. & 0\\ 0 & ... & 1 & -2 & 1 \end{pmatrix} \cdot \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix} = \frac{2m}{\hbar^2} \begin{pmatrix} V_1 - E & 0 & ... & 0\\ 0 & V_2 - E & ... & 0\\ 0 & 0 & ... & 0\\ 0 & 0 & ... & V_N - E \end{pmatrix} \cdot \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix} $$

Here is my question:

If I didn't have $V_i$ I would be using

scipy.linalg.eigh()

to solve this problem for $E$, however $E$ is now inside the second matrix. Is there a way to solve this equation finding the values of $E$ and relative sets of eigenvectors with Python?

Follow up

I implemented it as suggested in the comments as

$$ \begin{pmatrix} -2-\frac{2m}{\hbar^2}V_1 & 1 & 0 & ... & 0\\ 1 & -2-\frac{2m}{\hbar^2}V_2 & 1 & ... & 0\\ 0 & .&.&. & 0\\ 0 & 0 &... & 1 & -2-\frac{2m}{\hbar^2}V_N \end{pmatrix} \cdot \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix} = -\frac{2m}{\hbar^2}E \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix} $$

and used scipy.linalg.eigh().

Here is my code

import numpy as np
from matplotlib import pyplot as plt
from scipy.linalg import eigh

# number of discrete bins in which the domain is divided
N = 200

# create discretised arrays for the x axis and the value of the wave funciton
xaxis = np.linspace(-1, 1, N)
Vpot = np.linspace(-1, 1, N)

# factor with energy hbar and mass
m = 1
hbar = 1

# define the potential
def V(x):

    if np.abs(x)<0.5:
        val = 0.
    else:
        val = 4.

    return val

# fill an array for the potential
for i in range(len(xaxis)):
    Vpot[i] = V(xaxis[i])

# initialize NxN matrix
Hmat = [[0 for x in range(N)] for y in range(N)]

# fill Hmat like the following
#   -2-2m/h^2*V1 1 0 0 ...0
#   1 -2-2m/h^2*V2 1 0 ...0
#   0 1 -2-2m/h^2*V3 1 ...0
#   .....................0
#   .........1 -2-2m/h^2*VN
for row in range(N):
    for elem in range(N):
        if row == elem: Hmat[row][elem] = -2 -(2*m)/(hbar**2)*Vpot[elem]
        if np.abs(row-elem) == 1: Hmat[row][elem] = 1

# now I get the eigenvalue with w[i] and corresponding eigenvector v[i]
energies, psi = eigh(Hmat, b=None, eigvals_only=False, turbo=True)

# plot potential
#plt.plot(xaxis, Vpot)

# plot ground state
plt.scatter(xaxis, psi[0], s = 1)
plt.show()

Two questions

1 - If I just put the potential (Vpot) to zero I get the classical particle in a hole solution. The solution clearly resembles the analytical one, however I get the signs mixed somehow

I am not sure why this is happening.

2 - When I introduce the double level in the potential, I start getting something wrong, even here I don't understand what's happening.

Any ideas?

Answer 1

You could be doing the same mistake as here, On using Python to solve Time Independent Schrodinger Equation, the eigenfunctions have their values "pushed" to one of the boundaries?

With Python routines typically the columns are the eigenvectors and not the rows. That means you have to dereference your eigenvectormatrix like this [:,0] to obtain the first column/eigenvector. Derefencing by [0] gives you the first row and not the column.

The sign of the eigenvectors also doesn't matter. The sign/phase is just a convention but any eigenvector times a constant remains an eigenvector. The python routines do not always return the eigenvectors with the same sign. Real normed eigenvectors still have an underdetermined sign of $\pm1$. You can fix it by imposing a convention, for example requiring that the first coefficient is positive or picking the coefficient with highest magnitude as positive.

Answer 2

First of all, I suspect one can solve this problem analytically.

As for scipy.linalg.eigh() (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh.html), you just include V in matrix a.