<a href="https://colab.research.google.com/github/EnzoItaliano/calculoNumericoEmPython/blob/master/Lista_3.ipynb" target="_parent"></a>

Universidade Tecnológica Federal do Paraná  
Professor: Wellington José Corrêa  
Orientando: Enzo Dornelles Italiano  
Cálculo Numérico

# Solução de Sistemas Lineares

Antes de iniciar é necessário rodar uma vez o código do tópico **Códigos**

##Códigos


```python
!pip3 install prettymatrix
import math
import copy
import numpy as np
import prettymatrix
from prettytable import PrettyTable
from sympy import *

def Determinantes(A, ordem):
    k = 1
    tempA = [[]]
    while k <= ordem:
        for i in range(k):
            for j in range(k):
                tempA[i].append(A[i][j])
        if(round(np.linalg.det(tempA), 8) == 0): return (k,round(np.linalg.det(tempA), 8))
        if k < ordem:
            for i in range(k):
                for j in range(k):
                    tempA[i].pop()
            tempA.append([])
        k+=1
    return (0,round(np.linalg.det(tempA), 8))

def Transposta(A, ordem):
    k = ordem
    mT = []
    for c in range(k):
        mT.append([])
    for i in range(k):
        for j in range(k):
            mT[i].append( A[j][i] )
    soma = 0
    for i in range(k):
        for j in range(k):
            if(mT[i][j] != A[i][j]):
                soma += 1
    if soma == 0:
        return 1

def valoresG(A, G, ordem):
    k = ordem
    for i in range(k):
        for j in range(k):
            if(i == 0 and j == 0):
                G[0][0] = A[0][0] ** (1/2)
            elif(i == j):
                soma = 0
                for c in range(0, i):
                    soma += G[i][c]**2
                G[i][j] = (A[i][j] - soma)**(1/2)
            elif(j == 0 and i >= 1):
                G[i][0] = A[i][0]/G[0][0]
            elif(j == 1 and i >= 2):
                G[i][1] = (A[i][1] - G[i][0]*G[1][0])/G[1][1]
                G[i][1] = round(G[i][1], 8)
            elif(1 <= j and j < i):
                soma = 0
                for c in range(i):
                    soma += G[i][c]*G[j][c]
                G[i][j] = (A[i][j] - soma)/G[j][j]
    
    G = np.asarray(G)
    return G

def sistLinear(G, B, ordem):
    y = symbols('y:'+str(ordem))
    mY = []
    for i in range(len(y)):
        mY.append(y[i])
    D = np.linalg.det(G)
    tempG = G.copy()
    for j in range(ordem):
        for i in range(ordem):
            tempG[i][j] = B[i]
        tempD = np.linalg.det(tempG)
        tempG = G.copy()
        mY[j] = round(tempD/D, 8)
    mTemp = []
    for i in range(len(mY)):
        mTemp.append([mY[i]])
    mY = mTemp.copy()
    mY = np.asarray(mY)
    return mY

def valoresLU(L, U, A, ordem):
    for i in range(ordem):
        for j in range(ordem):
            if i == 0:
                U[0][j] = A[0][j]
            elif i == 1 and j >= 1:
                U[1][j] =  A[1][j] - L[1][0]*U[0][j]
            elif i <= j:
                soma = 0
                for c in range(i):
                    soma += L[i][c]*U[c][j]
                U[i][j] = A[i][j] - soma
        
        for j in range(ordem):
            if i == 0 and j >= 1:
                L[j][0] = A[j][0]/U[0][0]
            elif i == 1 and j >= 2:
                L[j][1] = (A[j][1] - L[j][0]*U[0][1])/U[1][1]
            elif j > i:
                soma = 0
                for c in range(i):
                    soma += L[j][c]*U[c][i]
                L[j][i] = (A[j][i] - soma)/U[i][i]

    return L, U

def pivotParcial(A,B,MA,MB,ordem,k):
    if A[k][k] == 0:
        for i in range(k+1,ordem):
            if A[i][k] != 0:
                temp1 = copy.copy(A[i])
                temp2 = copy.copy(B[i])
                temp3 = copy.copy(MA[i])
                temp4 = copy.copy(MB[i])
                A[i] = copy.copy(A[k])
                B[i] = copy.copy(B[k])
                MA[i] = copy.copy(MA[k])
                MB[i] = copy.copy(MB[k])
                A[k] = copy.copy(temp1)
                B[k] = copy.copy(temp2)
                MA[k] = copy.copy(temp3)
                MB[k] = copy.copy(temp4)
                break
    save = 0
    maior = -math.inf
    for i in range(ordem):
        if abs(A[i][k]) > maior and i >= k:
            maior = copy.copy(abs(A[i][k]))
            save = copy.copy(i)
    temp1 = copy.copy(A[save])
    temp2 = copy.copy(B[save])
    temp3 = copy.copy(MA[save])
    temp4 = copy.copy(MB[save])
    A[save] = copy.copy(A[k])
    B[save] = copy.copy(B[k])
    MA[save] = copy.copy(MA[k])
    MB[save] = copy.copy(MB[k])
    A[k] = copy.copy(temp1)
    B[k] = copy.copy(temp2)
    MA[k] = copy.copy(temp3)
    MB[k] = copy.copy(temp4)

def MA(A, B, ordem):
    MA = np.zeros((ordem,ordem))
    MB = np.zeros(ordem)

    for i in range(ordem):
        MB[i] = B[i]
        for j in range(ordem):
            MA[i][j] = A[i][j]

    mAum = np.zeros((ordem,ordem+1))
    for k in range(ordem-1):
        pivotParcial(A,B,MA,MB,ordem,k)
        for i in range(k+1,ordem):
            MB[i] = round(B[i] - (B[k] * A[i][k])/A[k][k], 8)
            for j in range(ordem):
                MA[i][j] = round(A[i][j] - (A[k][j]*A[i][k])/A[k][k], 8)


        for i in range(ordem):
            for j in range(ordem+1):
                if j == ordem:
                    mAum[i][j] = MB[i]
                else:
                    mAum[i][j] = MA[i][j]
        mAum = np.asarray(mAum)
        print(prettymatrix.matrix_to_string(mAum, name='Matriz aumentada = '))

        for i in range(ordem):
            B[i] = MB[i]
            for j in range(ordem):
                A[i][j] = MA[i][j]
    return MA, MB

def Cholesky(A, B):
    ordem = np.shape(A)
    if ordem[0] == 0: return print("Isso não é uma matriz")
    if ordem[0] == ordem[1]:
        det = Determinantes(A, ordem[0])
        if( det[1] <= 0):
            print("Determinante <= 0 para k = " + str(det[0]))
            return
        trans = Transposta(A, ordem[0])
        if( not trans ):
            return print("Matriz não é igual sua transposta")

        mG = np.zeros((ordem[0],ordem[1]))
        mG = valoresG(A, mG, ordem[0])
        print(prettymatrix.matrix_to_string(mG, name='G = '))

        mY = sistLinear(mG, B, ordem[0])
        print(prettymatrix.matrix_to_string(mY, name='Y = '))

        transpG = np.zeros((ordem[0],ordem[1]))
        for i in range(ordem[0]):
            for j in range(ordem[0]):
                transpG[i][j] = mG[j][i]

        transpG = np.asarray(transpG)

        mX = sistLinear(transpG, mY, ordem[0])
        print(prettymatrix.matrix_to_string(mX, name='X = '))

def LU(A,B):
    ordem = np.shape(A)
    if ordem[0] == 0: return print("Isso não é uma matriz")
    if ordem[0] != ordem[1]: return print("Esta matriz não é quadrada")
    det = Determinantes(A, ordem[0]-1)
    if( det[1] == 0):
        return print("Determinante = 0 para k = " + str(det[0]))

    L = np.eye(ordem[0])
    U = np.zeros((ordem[0],ordem[1]))

    mLU = valoresLU(L, U, A, ordem[0])

    mY = sistLinear(L, B, ordem[0])
    mX = sistLinear(U, mY, ordem[0])

    print(prettymatrix.matrices_to_string(L, U, names=['L','U']))
    print(prettymatrix.matrix_to_string(mY, name='Y = '))
    print(prettymatrix.matrix_to_string(mX, name='X = '))

def gaussParcial(A, B):
    ordem = np.shape(A)
    if ordem[0] == 0: return print("Isso não é uma matriz")
    if ordem[0] != ordem[1]: return print("Esta matriz não é quadrada")
    mA, mB = MA(A, B, ordem[0])
    mX = sistLinear(mA, mB, ordem[0])
    print(prettymatrix.matrix_to_string(mX, name='X = '))

def norma_inf(A):
    somaMaior = 0
    for i in range(len(A)):
        if abs(A[i]) > somaMaior:
            somaMaior = abs(A[i])
    
    return somaMaior

def condicaoEDD(A, ordem):
    for i in range(ordem):
        for j in range(ordem):
            if A[i][j] > A[i][i]:
              return false

def condicaoF(A, ordem, tipo):
    F = np.zeros((ordem,ordem))
    for i in range(ordem):
        for j in range(ordem):
            if i != j:
                F[i][j] = -(A[i][j]/A[i][i])
    print(prettymatrix.matrix_to_string(F, name='Matriz F = '))
    print("Cálculo da soma em valor absoluto de cada linha da matriz F:")
    if norma_linha(F) < 1:
        return true
    elif tipo == 1:
        print("Cálculo da soma em valor absoluto de cada coluna da matriz F:")
        if norma_coluna(F) < 1:
            return true
    else:
        return false

def Sanssenfeld(A, ordem):
    B = []
    maior = -math.inf
    for i in range(ordem):
        soma = 0
        for j in range(i):
            if i == 0: continue
            soma += abs(A[i][j]/A[i][i]) * B[j]
        temp = copy.copy(soma)
        soma = 0
        for j in range(i+1,ordem):
            soma += abs(A[i][j]/A[i][i])
        print("\u03B2", i, " = ", round(temp+soma,8))
        B.append(temp+soma)
        if B[len(B)-1] > maior:
            maior = B[len(B)-1]
    if maior < 1:
        return true
    else:
        return false

def norma_linha(A):
    ordem = np.shape(A)
    if ordem[0] == 0 or ordem[1] == 0: return print("Isso não é uma matriz")
    
    somaMaior = 0
    for i in range(ordem[0]):
        soma = 0
        for j in range(ordem[0]):
            soma += abs(A[i][j])
        if soma > somaMaior:
            somaMaior = soma
        print("Soma dos elementos da linha ", i+1, " = ", round(soma, 8))
    
    return somaMaior

def norma_coluna(A):
    ordem = np.shape(A)
    if ordem[0] == 0 or ordem[1] == 0: return print("Isso não é uma matriz")
    
    somaMaior = 0
    for j in range(ordem[0]):
        soma = 0
        for i in range(ordem[0]):
            soma += abs(A[i][j])
        if soma > somaMaior:
            somaMaior = soma
        print("Soma dos elementos da coluna ", j+1, " = ", soma)
    
    return somaMaior

def norma_euclidiana(A):
    ordem = np.shape(A)
    if ordem[0] == 0 or ordem[1] == 0: return print("Isso não é uma matriz")
    soma = 0
    for i in range(ordem[0]):
        for j in range(ordem[0]):
            soma += A[i][j]**2
        
    soma = soma ** (1/2)
    return soma

def GaussJacobi(A, B, X0, e):
    flag = 0
    ordem = np.shape(A)
    if ordem[0] == 0 or ordem[1] == 0: return print("Isso não é uma matriz")
    if condicaoEDD(A,ordem[0]) == false:
        flag += 1
        print("A matriz não é E.D.D.")
    if condicaoF(A, ordem[0], 1) == false:
        flag += 1
        print("A matriz não contém os requisitos quando feita a matriz F")
    if flag != 2:
        xk = []
        x_x = []
        xk.append(X0)
        x_x.append('-')
        controle = 0
        end_condition = 0
        while not end_condition:
            array = []
            for i in range(len(X0)):
                array2 = []
                temp = 0
                for j in range(len(X0)):
                    if i != j:
                        temp += A[i][j]*X0[j]
                array.append(round((1/A[i][i]) * (B[i] - temp),8))

            for i in range(len(array)):
                array2.append(array[i]-xk[controle][i])
            
            xk.append(array)
            x_x.append(norma_inf(array2)/norma_inf(array))
            X0 = array
            controle += 1
            if norma_inf(array2)/norma_inf(array) < e:
                end_condition = 1
        
        Table = PrettyTable(["k", "xk", "||x^(k+1) - x^(k)||/||x^(k+1)||"])
        for k in range(0, len(xk)):
            Table.add_row([k, xk[k], x_x[k]])
        
        print(Table)

def GaussSeidel(A, B, X0, e):
    flag = 0
    ordem = np.shape(A)
    if ordem[0] == 0 or ordem[1] == 0: return print("Isso não é uma matriz")
    if condicaoEDD(A,ordem[0]) == false:
        flag += 1
        print("Matriz não é E.D.D.")
    if condicaoF(A, ordem[0], 2) == false:
        flag += 1
        print("Não atende aos requisitos quando transformada na matriz F")
    print("Critério de Sanssenfeld:")
    if Sanssenfeld(A, ordem[0]) == false:
        flag += 1
        print("A matriz não atende ao critério de Sanssenfeld")
    if flag != 3:
        xk = []
        x_x = []
        xk.append(X0)
        x_x.append('-')
        controle = 0
        end_condition = 0
        while not end_condition:
            array = []
            array2 = []
            for i in range(len(X0)):
                temp = 0
                for j in range(len(X0)):
                    if i != j:
                        if j < i:
                            temp += A[i][j]*array[j]
                        else:
                            temp += A[i][j]*X0[j]
                array.append(round((1/A[i][i]) * (B[i] - temp),8))

            for i in range(len(array)):
                array2.append(array[i]-xk[controle][i])
            
            xk.append(array)
            x_x.append(norma_inf(array2)/norma_inf(array))
            X0 = array
            controle += 1
            if norma_inf(array2)/norma_inf(array) < e:
                end_condition = 1
        
        Table = PrettyTable(["k", "xk", "||x^(k+1) - x^(k)||/||x^(k+1)||"])
        for k in range(0, len(xk)):
            Table.add_row([k, xk[k], x_x[k]])
        
        print(Table)

    else:
        return print("A matriz não atende aos requisitos do método")
```

## 1. Métodos Diretos

### 1.1 Método de Cholesky

O procedimento a ser usado é Cholesky(A,B).

Exemplo: Resolva o sistema linear $A*X=B$ dado pelas matrizes abaixo:


```python
A = [[4,2,-4],
     [2,10,4],
     [-4,4,9]]
B = [0, 
     6, 
     5]
Cholesky(A,B)
```

Se tivéssemos escolhido a matriz A como sendo,


```python
A = [[4,2,-4],
     [-2,10,4],
     [-4,4,9]]
B = [0, 
     6, 
     5]
Cholesky(A,B)
```

Agora, se considerarmos,


```python
A = [[1,0,-3],
     [0,1,8],
     [-3,8,4]]
B = [0, 
     6, 
     5]
Cholesky(A,B)
```

### 1.2 Decomposição em LU

Aqui o procedimento é LU(A,B)

Exemplo: Usando o método de decomposição LU, resolva o sistema linear $A*X=B$ onde A e B são matrizes dadas a seguir:


```python
A = [[5,2,1],[-1,4,2],[2,-3,10]]
B = [-12, 20, 3]
LU(A,B)
```

Agora, se escolhermos uma matriz com algum menor principal nulo, por exemplo,


```python
A = [[5,2,3,4],
     [1,1,1,2],
     [4,1,2,0],
     [4,3,2,1]]
B = [1,
     2,
     3,
     4]
LU(A,B)
```

### 1.3 Método Eliminação de Gauss

O procedimento é gaussParcial(A,B)

Exemplo: Resolva o seguinte sistema abaixo onde as matrizes A e B são dadas a seguir com o método:


```python
A = [[1,2,3],
     [3,1,0],
     [0,3,4]]
B = [3, 4, 3]
gaussParcial(A,B)
```

## 2. Métodos Iterativos

### 2.1 Norma de Matriz

Os procedimentos abaixo nos dão respectivamente, a norma linha, norma coluna e a norma euclidiana de uma matriz A:

norma_linha(A)
norma_coluna(A)
norma_euclidiana(A)

Exemplo: Dada a matriz A, obtenha a norma linha, norma coluna e a norma euclidiana desta matriz.


```python
A = [[3,-5,7],
     [1,-2,4],
     [-8,1,-7]]
```


```python
print("||A||_inf = ", norma_linha(A))
```


```python
print("||A||_1 = ", norma_coluna(A))
```


```python
print("||A||_E = ", norma_euclidiana(A))
```

### 2.2 Método de Gauss Jacobi

Usaremos os seguintes procedimentos que

1. verificam se a matriz A é estritamente diagonalmente dominante e se os critérios da linha e da coluna são satisfeitos para a matriz F, onde $F_{ij}=0$ se $i=j$ e $F_{ij}=-A_{ij}*[A_{ii}]^{-1}$:

2. Uso do método de Gauss-Jacobi para resolver o sistema $A*X=B$ dados $X_0$ e tolerância
$\epsilon$:


GaussJacobi(A,B,X0,epsilon)

Exemplo: Utilize o método de Gauss - Jacobi para resolver o sistema $A*X=B$ com $X_0=[0,0,0]^t$ e
$\epsilon=10^{-2}$.

Solução: Vejamos se o algoritmo descrito nos permitirá usar  o método iterativo de Gauss-Seidel:


```python
A = [[10,2,1],[1,5,1],[2,3,10]]
B = [14, 11, 8]
X0 = [0,0,0]
GaussJacobi(A,B,X0,0.01)
```

### 2.3 Método de Gauss Seidel

Como na seção anterior, usaremos o procedimento com o intuito:

1. Verificar se a matriz A é estritamente diagonalmente dominante e se os critérios de Sanssenfeld, da linha para matriz F e da coluna para a matriz F são satisfeitos, onde a matriz F, onde $F_{ij}=0$ se i=j e $F_{ij}=-A_{ij}*[A_{ii}]^{-1}$:


2. Uso do método de Gauss-Seidel para resolver o sistema $A*X=B$ dados $X_0$ e tolerância
$\epsilon$:

GaussSeidel(A,B,X0,epsilon)

Exemplo: Utilize o método de Gauss - Jacobi para resolver o sistema $A*X=B$ com $X_0=[0,0,0]^t$ e
$\epsilon=10^{-2}$.

Solução: Vamos checar se é possível empregar o método iterativo de Gauss - Seidel:


```python
A = [[5,1,1],[3,4,1],[3,3,6]]
B = [5, 6, 0]
X0 = [0,0,0]
GaussSeidel(A, B, X0, 0.01)
```

1bbfc31f7e9db55afba6343342ba3266ac1f0dfe
ipynb
Jupyter Notebook
Lista_3.ipynb
EnzoItaliano/calculoNumericoEmPython
be3161b823955620be71e0f94a3421288fd28ef0
MIT
2019-12-28T21:23:00.000Z
2019-12-28T21:23:00.000Z
Lista_3.ipynb
EnzoItaliano/calculoNumericoEmPython
be3161b823955620be71e0f94a3421288fd28ef0
MIT
Lista_3.ipynb
EnzoItaliano/calculoNumericoEmPython
be3161b823955620be71e0f94a3421288fd28ef0
MIT
Qwen/Qwen-72B
1. YES
2. YES