0
$\begingroup$

I'm doing homework, which I have to build a code that resolves a system of linear equations.

But I'm having a problem. I have to use a matrix and a column matrix in the 3 functions. But, when I do a def main(), and call the 3 functions, these matrix are changing, and I cannot run the code.

For example: I start running the Gauss Elimination first. It runs very well. But after finished, the matrix is changed. And, when starts to run the Jacobi method, it fails, because the matrix is not the original.

How can I deal with it? Can I reset, or do something with the matrix?

from numpy import array, dot, diag
import numpy as np

'''def seidel(A,b):
    A = np.array(s)
    b = np.array(c)
    print ("Metodo de Gauss-Seidel")
    print("Sistema Linear:")
    for i in range(A.shape[0]):
        row = ["{}*I{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
        print(" + ".join(row), "=", b[i])
    print()

    iteracao = 100 #numero de iterações
    A [[0,1]] = A[[1,0]]    #troca da primeira com a segunda linha
    print("Depois de trocar as linhas, temos:")
    for i in range(A.shape[0]):
        row = ["{}*I{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
        print(" + ".join(row), "=", b[i])
    print()
    x = np.zeros_like(b)    #para x receber o numero de zeros iguais ao tamanho de b
    for cont in range(1, iteracao): #laço para fazer os calculos
        xnovo = np.zeros_like(x) #o novo valor de x, recebe o mesmo espaço do vetor b
        print("Iteração {0}: {1}".format(cont, x)) #mostro o valor atual de x
        for i in range(A.shape[0]):  #vai fazer o calculo de acordo com o tamanho da dimensão do array de A
            mult1 = np.dot(A[i, :i], xnovo[:i])  #aqui faz a multiplicação usando o comando dot. Pegamos do ponto i até o
            # iesimo elemento da Matriz A, e multiplica pelo até iesimo elemento do vetor x
            mult2 = np.dot(A[i, i + 1:], x[i + 1:]) #aplicamos o mesmo comando, porém com a alteração de começar
            # a conta a partir do iesimo elemento
            xnovo[i] = (b[i] - mult1 - mult2) / A[i, i] #aqui calculamos o valor do vetor x, devido aos valores anteriores
        if np.allclose(x, xnovo, rtol=1e-3): #se as variaveis x e xnovo tiverem valores pequenos, com uma tolerancia
            #pequena, a condição faz com que pare o laço. (exemplo: 1.0 e 0.999999 são proximos, se a tolerancia for 0.01)
            break
        x = xnovo

    print("Solução: {0}".format(x))
    erro = np.dot(A, x) - b
    print("Erro: {0}".format(erro))
    return 0'''

def jacobi(s,c):
    print ("Método de Jacobi")
    A = np.array(s)
    print(s)
    b = np.array(c)
    print("Sistema Linear:")
    for i in range(A.shape[0]):
        row = ["{}*I{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
        print(" + ".join(row), "=", b[i])
    print()

    iteracao = 100 #numero de iterações
    A [[0,1]] = A[[1,0]]    #troca da primeira com a segunda linha
    print("Depois de trocar as linhas, temos:")
    for i in range(A.shape[0]):
        row = ["{}*I{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
        print(" + ".join(row), "=", b[i])
    print()
    x = np.zeros_like(b)    #para x receber o numero de zeros iguais ao tamanho de b
    for count in range(iteracao): #laço para fazer os calculos
        print("Solução:", x)    #mostro o valor atual de x
        xnovo = np.zeros_like(x)    #o novo valor de x, recebe o mesmo espaço do vetor b

        for i in range(A.shape[0]): #vai fazer o calculo de acordo com o tamanho da dimensão do array de A
            mult1 = np.dot(A[i, :i], x[:i])    #aqui faz a multiplicação usando o comando dot. Pegamos do ponto i até o
            # iesimo elemento da Matriz A, e multiplica pelo até iesimo elemento do vetor x
            mult2 = np.dot(A[i, i + 1:], x[i + 1:])    #aplicamos o mesmo comando, porém com a alteração de começar
            # a conta a partir do iesimo elemento
            xnovo[i] = (b[i] - mult1 - mult2) / A[i, i] #aqui calculamos o valor do vetor x, devido aos valores anteriores

        if np.allclose(x, xnovo, atol=1e-3, rtol=0.): #se as variaveis x e xnovo tiverem valores pequenos, com uma tolerancia
            #pequena, a condição faz com que pare o laço. (exemplo: 1.0 e 0.999999 são proximos, se a tolerancia for 0.01)
            break

        x = xnovo
    print()
    print("Solução:", x)
    erro = np.dot(A, x) - b
    print("Erro:", erro)



def eliminGauss (s,c):
    print("Eliminação de Gauss")
    A1 = (s)
    b1 = (c)
    n = len(A1)  # vejo o tamanho da matriz
    M = A1  # coloco a matriz A em M
    i = 0  # contador para fazer a solucaoo do sistema de equacao

    for x in M:
        x.append(b1[i])  # aqui eu adiciono uma coluna, onde coloco os valores de b na ultima coluna
        i += 1

    for row in M:  # aqui eu printo, para mostrar como ficou a matriz com os valores de b
        print(row)
    print()

    for k in range(n):  # aqui eu começo a fazer as contas.
        print("Iteração", k)  # vou fazer k iteracoes, onde ocorrera as trocas de linhas e subtracao

        for i in range(k, n):  # nesse laço, eu vejo se o valor de [i][j] é maior que o [i][i] (pivotamento)
            if abs(M[i][k]) > abs(M[k][k]):  # se o valor for maior, eu faco a troca das linhas
                M[k], M[i] = M[i], M[k]
            else:
                pass  # se não for, eu continuo procurando

        for row in M:  # aqui eu printo, para mostrar as trocar que ocorrerao
            print(row)
        print()

        for j in range(k + 1, n):  # nesse laco, eu faco a eliminacao, e deixo zerado os valores abaixo do pivo
            c = float(M[j][k]) / M[k][k]  # aqui eu calculo o valor para zerar as colunas abaixo do pivo
            for m in range(k, n + 1):
                M[j][m] -= c * M[k][m]  # aqui eu faco a eliminacao, usando a subtracao do valor menos o pivo

        for row in M:  # aqui eu printo para mostrar os valores que estao sendo calculados
            print(row)
        print()
    x = [0 for i in range(n)]  # aqui eu monto um vetor com n casas, onde serão calculados os valores das correntes

    print("n = ", n)  # printo o tamanho n, que é o tamanho da Matriz A
    print("x = ", x)  # printo o tamanho do vetor que contem as correntes
    print()
    for row in M:  # aqui eu mostro como a matriz ficou, após a eliminacao
        print(row)
    print()

    x[n - 1] = float(M[n - 1][n]) / M[n - 1][n - 1]  # vou calcular as correntes com esse valor
    for i in range(n - 1, -1, -1):  # aqui, eu calculo o valor de cada corrente, onde eu tenho que colocar o valor da
        # corrente i3, na segunda coluna, para achar a i2, e colocar os valores da corrente i2 e i3 na primeira linha
        # para entao achar o valor da i1 (basicamente ir substituindo os valores encontrados)
        z = 0
        for j in range(i + 1, n):
            z = z + float(M[i][j]) * x[j]
        x[i] = float(M[i][n] - z) / M[i][i]
    print("Erro: ", x)
    print()
    i = 0
    while i <= 2:
        print("O valor da corrente i", i + 1, " é =", round(x[i], 4), "em A")  # printando os valores em separado
        i += 1

def main():
    s = ([[0.0, 5.0, -1.0], [11.0, 0.0, -1.0], [1.0, -1.0, -1.0]])  # o valor da Matriz A
    c = ([5.0, 14.0, 0.0])  # valor da matriz b
    eliminGauss(s,c)
    print()
    jacobi(s,c)
    print()
    #seidel(s,c)
main()

When I run the EliminGauss, I got this answer:

[11.0, 0.0, -1.0, 14.0], [0.0, 5.0, -1.0, 5.0], [0.0, 0.0, -1.1090909090909091, -0.2727272727272727]

After this, the code starts to run the function Jacobi. And it has to use the matrix S, not the answer from the EliminGauss.

This is the problem. I really don't know what happens here

$\endgroup$

put on hold as off-topic by Brian Borchers, Mauro Vanzetto, Christian Clason, Kirill, Biswajit Banerjee yesterday

This question appears to be off-topic. The users who voted to close gave this specific reason:

If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ It's impossible for us to debug your code based on your description. Please include a sample code to reproduce the thing that you are asking. $\endgroup$ – Alone Programmer Sep 20 at 23:59
  • $\begingroup$ Ok, I'll put it. But all the coments are in portuguese. But I'll show the things that happen $\endgroup$ – Daniel M M Sep 21 at 0:16
  • $\begingroup$ Some parts of this code are in a comment, because I was trying to see the problems $\endgroup$ – Daniel M M Sep 21 at 0:23
  • 1
    $\begingroup$ You can remove the parts of the code that are not relevant to the question. Indeed, a Minimal Working Example is preferred. Regarding the comments in Portuguese, you can rewrite them in English. $\endgroup$ – nicoguaro Sep 21 at 13:36
  • $\begingroup$ In EliminGauss, you write A1=(s). This creates another reference to the data stored in s, not a new copy with the same data. The same is true for M=A1. So any changes that you make to M are going to change the original s. You should do a proper fix of this, but you can see this is the issue by swapping the order of Jacobi and EliminGauss. In Jacobi, you call np.array which makes a whole new object so you don't alter your original list and both calls should give the same result (assuming the functions are implemented correctly) $\endgroup$ – Tyberius Sep 21 at 18:56

Browse other questions tagged or ask your own question.