Python: Successive-over-relaxation (Gauss-Seidel with relaxation) code is not converging -

i trying implement quantitative mobility spectrum analysis (qmsa) in python 2.7 + numpy + matplotlib. qmsa tool analysis of magnetic field-dependent hall-effect measurements. in semiconductor, more 1 type of carrier (electron or hole) different concentration , mobility can found. qmsa can seperate number of different types of carriers, density, mobility, , sign. there versions of improved version i-qmsa. however, building standart qmsa 0 hard work do.

the job highly know. there alot of scientific articles subject. however, because of copyrights of each article held publisher can not share you. mostly, can reach them university account. there thesis like:

1.) digital.library.txstate.edu/bitstream/handle/10877/4390/cunningham-thesis.pdf?sequence=1

2.) wrap.warwick.ac.uk/56132/1/wrap_thesis_kiatgamolchai_2000.pdf

3.) etd.lib.nsysu.edu.tw/etd-db/etd-search/getfile?urn=etd-0629104-152644&filename=etd-0629104-152644.pdf (i think in chinese)

4.) fbetezbankasi.gazi.edu.tr/pdf-indir/22233741 (i think in turkish. equations given in qmsa manual listed in thesis @ pages 73-75)

code trying successive-over-relaxation (sor) iterative approach done. firstly, prepare simple code produce artificial experimental data of magnetic field dependent conductivity tensor sigmaxx(b) , sigmaxy(b). experimental input , predefine mobility values, code running...

for in range (0,n,1):     bxx[i] = data[i][1]     bxy[i] = data[i][2]     j in range (0,m,1):         if data[i][0] == 0:             data[i][0] = 0.001         axx[j,i]=1/(1+(mobin[j]**2)*(data[i][0]**2))         axy[j,i]=(mobin[j]*data[i][0])/(1+(mobin[j]**2)*(data[i][0]**2)) 

here, bxx, bxy, mobin , data[i][0] experimental sigmaxx, experimental sigmaxy, predefined mobility values taken text file , experimental magnetic field points, respectively. therefore trying solve 2 equations sor in form of ax=b. xx part of problem a, x , b renamed axx, solxx , bxx. xy part of problem a, x , b renamed axy, solxy , bxy.

for sor, need parameter called omega. found optimum omega values gauss-seidel (here showing xx part of conductivity. same procedures done xy):

print "iterations finding omega..." it_count in range(1,501):     in range(0,n,1):         s1xx = np.dot(axx[i, :i], solxx_new[:i])         s2xx = np.dot(axx[i, + 1:], solxx[i + 1:])         solxx_new[i] = (bxx[i] - s1xx - s2xx) / axx[i, i]     dx = sqrt(np.dot(solxx_new-solxx,solxx_new-solxx))     in range(0,n,1):         solxx[i] = solxx_new[i]     if it_count == k: dx1 = dx     if it_count == k + 1:         dx2 = dx         omegaxx = 2.0/(1.0 + sqrt(abs(1.0 - (dx2/dx1))))         break print "i think best omega value these xx calculations ", omegaxx 

this "finding optimum omega" procedure taken kiusalaas, jaan. numerical methods in engineering python 3. cambridge university press, 2013, page 83.

after finding omega, time same iterations done sor:

for it_count in range(iteration_limit):    in range(0,n,1):        s1xx = np.dot(axx[i, :i], solxx_new[:i])        s2xx = np.dot(axx[i, + 1:], solxx[i + 1:])        solxx_new[i] = (1-omegaxx)*solxx[i-1]+omegaxx*(bxx[i] - s1xx - s2xx) / axx[i, i]    if np.allclose(solxx, solxx_new, rtol=1e-9):        break    print "iteration:",it_count    in range(0,n,1):        solxx[i] = solxx_new[i] 

then, calculated conductivity spectrum values each mobility with:

for in range (0,n,1):     if == 0:         deltamob = 100     else:         deltamob = mobin[i] - mobin[i-1]     sn[i] = abs((solxx[i] - solxy[i]))/(2*deltamob*1.6e-19)     sp[i] = abs((solxx[i] + solxy[i]))/(2*deltamob*1.6e-19)     x[i] = mobin[i]     b[i] = data[i][0] 

then x vs sn , x vs sp must mobility spectrums. thing can single gaussian peak. , without hole carrier, electron , hole spectrums identical. problem solxx , solxy getting larger values after each iteration. problem may caused sor code written python 3. using python 2.7.

i can send files if necessary.

thanks responses.