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

Crank Nicolson method

  1. Jun 21, 2011 #1
    I have tried discretizing the Schrodinger equation using the crank nicolson method and then have tried to solve the resulting tri-diagonal matrix using the algorithm on Wikipedia.


    However, I am not getting the required result and cannot put my finger on it. Please help! I have attached the code below.

    h is my space step and deltat is my time step. J = total number of space steps and K = total number of time steps

    Code (Text):
    import numpy as np



    h = 0.01
    deltat = 0.001

    J = 500
    K = 500

    phi_a = np.zeros(shape=(K,J), dtype = np.complex)


    for p in range(1,J-1):
       phi_a[0][p] = np.sin(np.pi*(p+1)/J)


    a = np.ones(J, dtype=np.complex)/(4*h**2)
    a[0], a[J-1] = 0,0

    c = np.ones(J, dtype=np.complex)/(4*h**2)
    c[0], c[J-1] = 0,0


    def b(k):
       b = np.ones(J, dtype=np.complex)
       for p in range(1,J-1):
          b[p] = (1j/(deltat)) - (1/(2*h**2))


    def d(k):
       d = np.ones(J, dtype=np.complex)

       for p in range(1,J-2):
          d[p] = -(phi_a[k][p+1]/(4*h**2)) + (phi_a[k][p]*((1j/deltat) + 1/(2*h**2))) - (phi_a[k][p-1]/(4*h**2))

       d[0] = phi_a[k+1][0]
       d[J-1] = phi_a[k+1][J-1]


    def tdms(a,c,k,s,b,d):
       b,d = b(k),d(k)

       for p in range(1,J-1):
          c[p] = c[p]/(b[p] - c[p-1]*a[p])

       for p in range(J):
          d[p] = (d[p]-d[p-1]*a[p])/(b[p] - c[p-1]*a[p])

       for p in range(J-2,0,-1):
          s[k+1][p] = d[p] - c[p]*s[k+1][p+1]



    for k in range(0,K-1):


    anphi_a = open('test.txt','w')
    for W in range(K):
       for Q in range(J-1):
          anphi_a.write(str(phi_a[W][Q]) + ', ')
       anphi_a.write(str(phi_a[W][J-1]) + '\n')
    This is written in Python, which is nearly pseudo-code, so most of you should be able to understand it. Please help!
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted