On using Python to solve Time Independent Schrodinger Equation, the eigenfunctions have their values "pushed" to one of the boundaries?

I am having trouble using numerical methods to solve Time Independent Schrodinger Equation. I am considering a quartic potential function: $$V(x) = x^4 -4x^2.$$ $$-\frac{d^2\psi(x)}{dx^2} + V(x) \psi(x) = E \psi(x)$$

I wish to get a few solutions of the eigenproblem (about 150).

Here is the code I have written:

import numpy as np
import matplotlib.pyplot as plt

def solveTISE(xmin,xmax,h,potFunc,params):
xvec = np.arange(xmin,xmax,h)
Nx = len(xvec)
potFunctVec = potFunc(xvec,params)
mat = -((np.tri(Nx,Nx,1)-np.tri(Nx,Nx,-2)-3.*np.eye(Nx))/h**2-np.diag(potFunctVec))
print(mat)
eigenValues,eigenVectors = np.linalg.eig(mat)
idx = eigenValues.argsort()[::1]
eigenValues = eigenValues[idx]
eigenVectors = eigenVectors[:,idx]
return eigenValues,eigenVectors,xvec

def doubleWell(x,params):
return np.piecewise(x, [x < 0, x >= 0], [lambda x : params[2]*(x**4 - params[0]*x**2), lambda x :  params[3]*(x**4 - params[1]*x**2)] )

betaValue = 0.01;
hValue = 0.01;
eigEnergy1i, eigFunc1i, xvec1i = solveTISE(-10,10,hValue,doubleWell,[4,4,1,1]);

plt.plot(xvec1i,eigFunc1i[0]);
plt.show()
plt.plot(xvec1i,eigFunc1i[1]);
plt.show()


However, I get the following output:

When I changed the limits of x to : -20 and 20, the resulting wavefunction (the 0th index one/lowest energy state) looked like this:

I don't know why the wavefunction is being "pushed" to the right boundary.

• Hi, can you provide a little more info? Define "TISE" so we know what problem you're attacking. How did you get data for [-60, +60] if your input was [-20,+20] ? By the way, since your potential function is quadratic in $x^2$ you may be able to simplify your work. May 21, 2021 at 12:14
• I'd recommend testing code first with problems where you know the solution analytically to test it. For example a simple harmonic potential. Does your script work there ? May 21, 2021 at 13:06
• Hi @CarlWitthoft! Thanks for the remarks. I have edited the questions to add the correct image and stated that TISE is the time-independent Schrodinger equation. May 21, 2021 at 13:07
• Hi @HansWurst! I did test the code for Infinite Potential Well and worked fine for it. However, the code had the same issues for harmonic potential. May 21, 2021 at 13:08
• Please note that generally, if you want help with your code, the code should be complete: i.e. runnable without second-guessing what is missing. This time I've completed your code in my edit. May 21, 2021 at 13:59

Your script looks fine on first glance. Your problem is that your eigenfunctions are the columns of your eigFunc1 Matrix and not the rows. Use this to plot them.

plt.plot(xvec1i,eigFunc1i[:,n]);


where n is the integer of the corresponding eigenfunction. I.e. to plot the groundstate use

plt.plot(xvec1i,eigFunc1i[:,0]);


Also note that some_array[::1] doesn't change the array, as you you are stepping over all values with a step size of 1.

EDIT:

Due to popular demand a small plot of the first eigenvectors n=0 and n=1 at the height of their eigenvalue as baseline which can be read off at the left hand y-axis.

The figure was created with

fig, ax = plt.subplots()
ax.plot(xvec1i,eigFunc1i[:,0]+eigEnergy1i[0], c="b", lw=3);
ax.plot(xvec1i,eigFunc1i[:,1]+eigEnergy1i[1], c="orange", lw=3);
ax.axhline(eigEnergy1i[0],c="b", alpha=0.5, lw=3, label="n=0")
ax.axhline(eigEnergy1i[1], c="orange", alpha=.5, lw=3, label="n=1")
ax_twinx = ax.twinx()
ax_twinx.plot(xvec1i,doubleWell(xvec1i,[4,4,1,1]))
ax.legend()
ax.set_xlim(-10, 10)

• Thanks a lot for the clarification! It is working now. May 21, 2021 at 13:47
• @PhysicsKidDyinginSchool Well then, update your question to include the correct plots! I'm curious to see what they look like! :) May 21, 2021 at 13:49
• @Philip Ideally Hanswurst should post the correct plots, rather than including the "answer" material in the question. Given the completenesss of the code & parameters, we all should be able to plot on our own :-) May 21, 2021 at 14:09
• @CarlWitthoft Good point! I thought it was a little forward to ask HansWurst to do it, which is why I asked the OP :) But you are right, of course. May 21, 2021 at 15:08

Your another (small) problem is that your values of $$x$$, called xvec, are not in the range $$\{-10, -9.99, \dots, +9.99, +10\}$$. Instead they lack the $$+10$$, because arange(begin,end,step) returns the range [begin,end). The result for your code is that the states that should have very close energies will instead differ more than expected, due to the broken symmetry. The fix for this is to replace your generation of xvec with

    xvec = np.arange(xmin,xmax+h,h)


This is much better, but it still gives you not ideal values, so the precision of the eigenvalues will be a bit smaller than you may like. Here's the output of print(xvec) after setting np.set_printoptions(precision=17):

[-10.                 -9.99               -9.98              ...
9.979999999999574   9.989999999999574   9.999999999999574]


I suppose it's due to the way arange() is implemented. To fix this, you can try to make sure that the step is integral, so that you don't accumulate too much rounding error. An example (just a dirty hack for demonstration, not intended for production):

    xvec = np.arange(xmin*1000,(xmax+h)*1000,h*1000)/1000


This code does generate

[-10.    -9.99  -9.98 ...   9.98   9.99  10.  ]

• You sure you're not just running into floating point precision there? There's no way Python stores the apparent value "9.98" as a float without the binary precision error In R, for example,  x <- rep(-10,10); for (j in 2:10) x[j] = x[j-1]+0.01 sprintf("%1.20f", x) # [1] "-10.00000000000000000000" "-9.99000000000000021316" [3] "-9.98000000000000042633" "-9.97000000000000063949" [5] "-9.96000000000000085265" "-9.95000000000000106581" [7] "-9.94000000000000127898" "-9.93000000000000149214" [9] "-9.92000000000000170530" "-9.91000000000000191847"  May 21, 2021 at 14:19
• @CarlWitthoft indeed, almost all of these values are not representable in binary64. But this is no excuse for using a numerically poor algorithm of generating a range. If you just accumulate 0.01 like x[n]=x[n-1]+step, you'll get worse precision than if you calculate it as x[n]=xFirst+(xLast-xFirst)*n/(N-1). This is because roundoff error happens on each iteration when you accumulate. Unfortunately, it seems that numpy.arange does it the former, poor, way. May 21, 2021 at 14:41