How to Correctly Simulate an Infinite Well in Python?

[svletana]

I'm trying to get the eigenfunctions and eigenvalues (energies) of an infinite well in Python, but I have a few things I can't seem to fix or don't understand...
Here's the code I have:

from numpy import *
from numpy.linalg import eigh
import matplotlib.pyplot as plt
from __future__ import division

def Laplacian(x):
    h = x[1]-x[0] # assume uniformly spaced points
    n = len(x)
    M = -2*identity(n,'d')
    for i in range(1,n):
           M[i,i-1] = M[i-1,i] = 1
    return M/h**2

L = 6
x = linspace(0, L, 100)
T = -0.5*Laplacian(x)
H = T
E,U = eigh(H)

n = 1 # change to see another state
plt.plot(x,U[:,n-1],'o')
plt.plot(x, sqrt(2/L)*sin(n*pi*x/L),'o')
plt.show()

First of all, the eigenfunctions U[:,i] appear to be normalized if I check, but they don't match the analytic solution. Also, the eigenvalues don't match the values of energy for the infinite well 0.5(πn/L)².

What could I be missing? The math seems to be fine.
Would it make it better to use scipy.linalg instead of numpy's?

Thanks!

 

[Mark44]

def Laplacian(x):
    h = x[1]-x[0] # assume uniformly spaced points
    n = len(x)
    M = -2*identity(n,'d')
    for i in range(1,n):
           M[i,i-1] = M[i-1,i] = 1
    return M/h**2

Is that a typo in the next to last line? Should it instead be M[i, i-1] = M[i-1, i] + 1?

BTW, do you have a typo in your user name? svetlana is a name I've seen often, not svletana. If you've made a mistake, I can change it.

 

[svletana]

def Laplacian(x):
    h = x[1]-x[0] # assume uniformly spaced points
    n = len(x)
    M = -2*identity(n,'d')
    for i in range(1,n):
           M[i,i-1] = M[i-1,i] = 1
    return M/h**2

Is that a typo in the next to last line? Should it instead be M[i, i-1] = M[i-1, i] + 1?

BTW, do you have a typo in your user name? svetlana is a name I've seen often, not svletana. If you've made a mistake, I can change it.

The code should be

    for i in range(1,n):
            M[i,i-1] = 1
            M[i-1,i] = 1
    return M/h**2

And my username is not a typo, I normally write it like this because svetlana is most likely taken :)

 

[Mark44]

Now that I understand that part of your code, your two for loop versions do the same thing. What the loop seems to be doing is putting 1's on the diagonals next to the main diagonal of your identity matrix. The only thing that comes to mind that could be causing a problem is a possible off-by-one error. Keep in mind that range(1, n) consists of the numbers 1, 2, 3, ..., n-1, but doesn't include n. Make sure that your loop is setting all of the entries in the matrix that you want to have set.

Other than that, I don't see anything.

 

[ChrisVer]

can you explain me what this part does?

    for i in range(1,n):
            M[i,i-1] = 1
            M[i-1,i] = 1

I mean what is the array you produce?

 

[hilbert2]

I remember encountering a similar problem when doing this... Try to multiply the kinetic energy matrix T with some relatively large number, like 100, and see if the solutions then match the analytical ones (something like this worked, if I remember correctly).

 

[StoneTemplePython]

A couple things:

it's preferable to use "import numpy as np" -- importing * is generally frowned on. (For instance if you use "sum" after importing *, are you using np.sum() or python's built in sum)?

Second, the (finite) graph Laplacian is well understood to be symmetric positive semi-definite with a minimum eigenvalue $\lambda_n = 0$. That is $\mathbf L \mathbf 1 = \lambda_n \mathbf 1 = \mathbf 0$, and since real symmetric, $\mathbf 1^T \mathbf L = \lambda_n \mathbf 1^T = \mathbf 0^T$.

Just before your return statement in your Laplacian function, try adding the following two lines of code

ones_v = ones(M.shape[0])
print dot(ones_v,M)
# if importing numpy as np:
# ones_v = np.ones(M.shape[0])
# print np.dot(ones_v, M)

but this doesn't return the zero vector, which is very strong hint as to what's wrong here

--- alternatively, you could insert similar code just below your assignment to T. Either way the ones vector is not in the nullspace of your Laplacian which means it is not a Laplacian