# Crank Nicolson method

1. Jun 21, 2011

### WiFO215

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!

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