gauss.py

import numpy as np

matrix01 = np.matrix([
    [3, -0.1, -0.2],
    [0.1, 7, -0.3],
    [0.3, -0.2, 10],
])

solution01 = np.matrix([
    [7.85, -19.3, 71.4]
])

initial01 = np.matrix([
    [0, 0, 0]
])

error01 = 0.001

def abs(matrix, index):
    return matrix.item(np.abs(index))

def converge(m, verbose=False):
    cond1calc = abs(m, 1) + abs(m, 2)
    cond1 = abs(m, 0) >= cond1calc

    cond2calc = abs(m, 3) + abs(m, 5)
    cond2 = abs(m, 4) >= cond2calc

    cond3calc = abs(m, 6) + abs(m, 7)
    cond3 = abs(m, 8) >= cond3calc

    if verbose:
        print("")
        print(abs(m, 0), ">=", abs(m, 1), "+", abs(m, 2))
        print(abs(m, 0), ">=",round(cond1calc, 4), "(", cond1, ")")
        print("")
        print(abs(m, 4), ">=", abs(m, 3), "+", abs(m, 5))
        print(abs(m, 4), ">=",round(cond2calc, 4), "(", cond2, ")")
        print("")
        print(abs(m, 8), ">=", abs(m, 6), "+", abs(m, 7))
        print(abs(m, 8), ">=",round(cond3calc, 4), "(", cond3, ")")
        print("")

    if cond1 and cond2 and cond3:
        return True
    return False

def jacobi(m, s, i, e, verbose=False):

    if verbose:
        print("")
        print("JACOBI:")
        print("")

    if not converge(m, verbose):
        print("Matrix does not have dominant diagonal")
        return

    ok = False
    x = i.item(0)
    y = i.item(1)
    z = i.item(2)
    s0 = s.item(0)
    s1 = s.item(1)
    s2 = s.item(2)

    counter = 1

    while not ok:
        _x = (s0 - m.item(1) * y - m.item(2) * z) / m.item(0)
        _y = (s1 - m.item(3) * x - m.item(5) * z) / m.item(4)
        _z = (s2 - m.item(6) * x - m.item(7) * y) / m.item(8)

        cond1 = np.absolute(_x - x) <= e
        cond2 = np.absolute(_y - y) <= e
        cond3 = np.absolute(_z - z) <= e

        if verbose:
            print("------------------------------------------")
            print("#", counter)
            print("x = (", s0, "-",round(m.item(1) * y, 4), "-", round(m.item(2) * z, 4), ") / ", round(m.item(0), 4)," = ", round(_x,4))
            print("y = (", s1, "-",round(m.item(3) * x, 4), "-", round(m.item(5) * z, 4), ") / ", round(m.item(4), 4)," = ", round(_y,4))
            print("z = (", s2, "-",round(m.item(6) * x, 4), "-", round(m.item(7) * y, 4), ") / ", round(m.item(8), 4)," = ", round(_z,4))
            print("")
            print('|',round(_x,4),"-",round(x,4),'| <=',e)
            print(round(np.absolute(_x - x),4), '<=',e, "(",cond1,")")
            print("")
            print('|',round(_y,4),"-",round(y,4), '| <=',e)
            print(round(np.absolute(_y - y),4), '<=',e, "(",cond2,")")
            print("")
            print('|',round(_z,4),"-",round(z,4), '| <=',e)
            print(round(np.absolute(_z - z),4), '<=',e, "(",cond3,")")
            print("")

        if cond1 and cond2 and cond3:
            ok = True
            if verbose:
                print("------------------RESULT------------------")
                print("x =", _x)
                print("y =", _y)
                print("z =", _z)
                print("------------------------------------------")
            return [_x, _y, _z]
        else:
            x = _x
            y = _y
            z = _z
            counter = counter + 1



def seidel(m, s, i, e, verbose=False):
    if verbose:
        print("")
        print("SEIDEL:")
        print("")

    if not converge(m, verbose):
        print("Matrix does not have dominant diagonal")
        return

    ok = False
    x = i.item(0)
    y = i.item(1)
    z = i.item(2)
    s0 = s.item(0)
    s1 = s.item(1)
    s2 = s.item(2)

    counter = 1

    while not ok:

        _x = (s0 - m.item(1) * y - m.item(2) * z) / m.item(0)
        newX = _x
        _y = (s1 - m.item(3) * newX - m.item(5) * z) / m.item(4)
        newY = _y
        _z = (s2 - m.item(6) * newX - m.item(7) * newY) / m.item(8)

        cond1 = np.absolute(_x - x) <= e
        cond2 = np.absolute(_y - y) <= e
        cond3 = np.absolute(_z - z) <= e

        if verbose:
            print("------------------------------------------")
            print("#", counter)
            print("x = (", s0, "-",round(m.item(1) * y, 4), "-", round(m.item(2) * z, 4), ") / ", round(m.item(0), 4)," = ", round(_x,4))
            print("y = (", s1, "-",round(m.item(3) * newX, 4), "-", round(m.item(5) * z, 4), ") / ", round(m.item(4), 4)," = ", round(_y,4))
            print("z = (", s2, "-",round(m.item(6) * newX, 4), "-", round(m.item(7) * newY, 4), ") / ", round(m.item(8), 4)," = ", round(_z,4))
            print("")
            print('|',round(_x,4),"-",round(x,4),'| <=',e)
            print(round(np.absolute(_x - x),4), '<=',e, "(",cond1,")")
            print("")
            print('|',round(_y,4),"-",round(y,4), '| <=',e)
            print(round(np.absolute(_y - y),4), '<=',e, "(",cond2,")")
            print("")
            print('|',round(_z,4),"-",round(z,4), '| <=',e)
            print(round(np.absolute(_z - z),4), '<=',e, "(",cond3,")")
            print("")

        if cond1 and cond2 and cond3:
            ok = True
            if verbose:
                print("------------------RESULT------------------")
                print("x =", _x)
                print("y =", _y)
                print("z =", _z)
                print("------------------------------------------")
            return [_x, _y, _z]
        else:
            x = _x
            y = _y
            z = _z
            counter = counter + 1



jacobi(matrix01, solution01, initial01, error01, True)

seidel(matrix01, solution01, initial01, error01, True)