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.

    http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

    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

    #----------------------------------------------------#

    #GIVEN DATA

    h = 0.01
    deltat = 0.001

    J = 500
    K = 500

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

    #BOUNDARY CONDITIONS#

    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))

       return(b)

    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]
       
       return(d)

    #----------------------------------------------------#

    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]

       return

    #----------------------------------------------------#

    for k in range(0,K-1):
       tdms(a,c,k,phi_a,b,d)

    #----------------------------------------------------#

    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')
    anphi_a.close()
    This is written in Python, which is nearly pseudo-code, so most of you should be able to understand it. Please help!
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

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



Similar Discussions: Crank Nicolson method
  1. Bisection Method (Replies: 2)

  2. Newton Method (Replies: 1)

  3. Method of Complements (Replies: 3)

Loading...