Solving Schrodinger Eq. w/ Crank Nicolson Method

[WiFO215]

I have tried discretizing the Schrödinger 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

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!

 

[lofgran]

Hi there,

It's difficult to pinpoint the exact issue without more information, but here are a few suggestions that may help:

1. Double check your boundary conditions: Make sure they are correctly implemented and are consistent with the initial conditions.

2. Check your indexing: Make sure that you are using the correct indices when accessing and assigning values to arrays. A simple mistake in indexing can lead to unexpected results.

3. Consider using a different algorithm: The tridiagonal matrix algorithm is known to have stability issues for certain types of matrices. You may want to try using a different algorithm, such as the LU decomposition method, to solve the resulting matrix.

4. Check for convergence: Make sure that your solution is converging to a steady state. If it is not, then there may be a problem with your discretization or algorithm.

5. Debug your code: Use print statements or a debugger to track the values of your variables and make sure they are behaving as expected.

I hope these suggestions help. Good luck with your research!