- #1
WiFO215
- 420
- 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
This is written in Python, which is nearly pseudo-code, so most of you should be able to understand it. Please help!
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:
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!