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

I am trying to implement Quantitative Mobility Spectrum Analysis (QMSA) in Python 2.7 + NumPy + MatPlotLib. QMSA is a tool for analysis of magnetic field-dependent Hall-effect measurements. In a semiconductor, more than one type of carrier (electron or hole) with different concentration and mobility can be found. QMSA can seperate the number of different types of carriers, their density, mobility, and sign. There are some versions of it like improved version i-QMSA. However, even the building the standart QMSA from zero is a hard work to do.

The job is highly know. There are alot of scientific articles about the subject. However, because of copyrights of each article is held by a publisher I can not share with you. Mostly, you can reach them with an university account. There are some thesis about it 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 it is in Chinese)

4.) fbetezbankasi.gazi.edu.tr/pdf-indir/22233741 (I think it is in Turkish. Equations which is given in QMSA manual is listed in thesis at pages 73-75)

Code will I am trying to do Successive-over-relaxation (SOR) iterative approach as originally done. Firstly, I prepare I simple code to produce artificial experimental data of magnetic field dependent conductivity tensor sigmaxx(B) and sigmaxy(B). With this experimental input and a predefine mobility values, code is running…

for i in range (0,n,1):
    bxx[i] = data[i][1]
    bxy[i] = data[i][2]
    for 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 and data[i][0] is experimental sigmaxx, experimental sigmaxy, predefined mobility values taken from a text file and experimental magnetic field points, respectively. Therefore we are trying to solve two equations with SOR in the form of Ax=b. For XX part of the problem A, x and b are renamed as Axx, solxx and bxx. For XY part of the problem A, x and b are renamed as Axy, solxy and bxy.

For the SOR, you need a parameter called omega. I found optimum omega values with GAuss-Seidel (here I am showing this for XX part of the conductivity. Same procedures are also done for XY):

print "Iterations for finding omega..."
for it_count in range(1,501):
    for i in range(0,n,1):
        s1xx = np.dot(Axx[i, :i], solxx_new[:i])
        s2xx = np.dot(Axx[i, 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))
    for i 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 for these XX calculations is ", omegaxx

This “finding optimum omega” procedure is taken from Kiusalaas, Jaan. Numerical methods in engineering with Python 3. Cambridge university press, 2013, page 83.

After finding omega, this time same iterations are done with SOR:

for it_count in range(ITERATION_LIMIT):
   for i in range(0,n,1):
       s1xx = np.dot(Axx[i, :i], solxx_new[:i])
       s2xx = np.dot(Axx[i, 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
   for i in range(0,n,1):
       solxx[i] = solxx_new[i]

Then, I calculated conductivity spectrum values for each mobility with:

for i in range (0,n,1):
    if i == 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 and x vs sp must be your mobility spectrums. Only thing that I can get is a single Gaussian like peak. And even without any hole carrier, electron and hole spectrums are identical. Problem is solxx and solxy is getting larger values after each iteration. Problem may caused by the SOR code is written for Python 3. But I am using Python 2.7.

I can send files if necessary.

Thanks for any responses.


Source: python

Leave a Reply