Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python Simulating an infinite well

  1. Jun 1, 2017 #1
    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:
    Code (Python):
    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)2.

    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!
     
  2. jcsd
  3. Jun 1, 2017 #2

    Mark44

    Staff: Mentor

    Code (Python):
    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.
     
  4. Jun 1, 2017 #3
    The code should be

    Code (Python):
        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 :)
     
  5. Jun 1, 2017 #4

    Mark44

    Staff: Mentor

    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.
     
  6. Jun 1, 2017 #5

    ChrisVer

    User Avatar
    Gold Member

    can you explain me what this part does?
    Code (Python):

        for i in range(1,n):
                M[i,i-1] = 1
                M[i-1,i] = 1
     
    I mean what is the array you produce?
     
  7. Jun 14, 2017 #6

    hilbert2

    User Avatar
    Science Advisor
    Gold Member

    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).
     
  8. Jun 14, 2017 #7

    StoneTemplePython

    User Avatar
    Gold Member

    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


    Code (Python):
    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Simulating an infinite well
  1. Quantum Simulation (Replies: 5)

  2. Physics simulation (Replies: 5)

  3. Simulating the world (Replies: 14)

  4. Simulation MCNP6 (Replies: 1)

Loading...