Add Gauss-Seidel (with relaxation) and Jacobi for solving Ax=b

[EmmanuelMess]

Chapter Request

Description

With Gauss-Seidel and Jacobi one can find solutions for Ax=b in an iterative way, it is very useful.

Preconditions for convergence to a solution should also be added, (checking if T=I−(N^−1)A, is dialgonally dominant, if any norm is < 1, if largest eigenvalues is <1).

Code

SciLab:

//Gauss-Seidel with relaxation
function x = SOR(A, b, x0, vv, eps)
    sz = size(A,1)
    x = x0
    xant = x
    suma = 0
    it = 1
    cont = 0
    
    while(abs(norm(x-xant)) > eps | cont == 0) 
    xant = x
        for i = 1:sz
            suma = 0
            for j = 1:i-1 
                suma = suma + A(i,j)*x(j)
            end
            
            for j = i+1:sz
                suma = suma + A(i,j)*x(j)
            end
            x(i) = (1 - vv)*xant(i) + vv/(A(i,i))*(b(i)-suma)
        end
     cont = cont + 1
    end
    
    mprintf("Cantidad de iteraciones: %d\n",cont);
    
endfunction

function x = jacobisolver(A,b,x0,eps) 
    
    sz = size(A,1)
    x = x0
    xant = x
    suma = 0
    //it = 1
    cont = 0
    
    while(abs(norm(x-xant)) > eps | cont == 0) 
    xant = x
        for i = 1:sz
            suma = 0
            for j = 1:sz 
                if (i <> j)
                    suma = suma + A(i,j)*xant(j)
                end
            end
            x(i) = 1/(A(i,i))*(b(i)-suma)
        end
     cont = cont + 1
    end
    
    mprintf("Cantidad de iteraciones: %d\n",cont);
    
endfunction

Python:

def gaussseidel(x0, A, b):#A esta en fc
    x1 = x0
    size = len(A)

    for k in range(1000):
        suma = 0
        i = 0
        for j in range(1, size):
            suma += A[i][j] * x0[j]

        x1[i] = 1 / A[i][i] * (b[i] - suma)

        for i in range(1, size-1):
            suma = 0
            for j in range(i):
                suma += A[i][j]*x1[j]

            for j in range(i+1, size):
                suma += A[i][j]*x0[j]

            x1[i] = 1 / A[i][i] * (b[i] - suma)

        suma = 0
        i = size-1
        for j in range(size-1):
            suma += A[i][j] * x1[j]

        x1[i] = 1 / A[i][i] * (b[i] - suma)

        x0 = x1


    return x1

For Algorithm Archive Developers

• [ ] This chapter can be added to the Master Overview (if it cannot be, explain why in a comment below -- lack of sufficient technical sources, not novel enough, too ambitious, etc.)
• [ ] There is a timeline for when this chapter can be written
• [ ] The chapter has been added to the Master Overview
• [ ] The chapter has been written (Please link the PR)