Vamos considerar aqui os métodos para solução de sistemas lineares com $n$ equações e $n$ incógnitas da forma $$ \begin{array}{ccccc} a_{11}x_1 & + a_{12}x_2 & + \cdots & + a_{1n}x_n& = b_1 \\ \vdots & & & & \vdots \\ a_{n1}x_1& + a_{n2}x_2& + \cdots & + a_{nn}x_n & = b_n \end{array}$$ ou na forma matricial: $$ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} &\cdots & a_{nn} \end{bmatrix}\begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ \vdots \\ b_n \end{pmatrix}$$
Em primeiro lugar caracterizamos os sistemas fáceis de resolver. Um sistema é fácil de resolver quando a matriz de coeficientes $A$ é triangular superior, isto é quando tivermos: $$ a_{ij}=0 \text{ para } i > j. $$ Neste caso resolvemos o sistema começando da última equação e chegando até a primeira incógnita por substituição. Claro que o sistema só será possivel e determinado quando cada elemento $a_{ii}$ da diagonal for diferente de zero. E podemos aplicar a fórmula: $$ x_k = \frac{1}{a_{kk}}\left(b_k - \sum_{k < j \leq n} a_{kj}x_j\right) $$
import numpy as np
# para a maioria dos exemplos vamos usar o numpy
def BackStep(A,b):
""" Se a matriz for triangular superior resolve Ax=b
A é uma matriz triangular.
b é uma lista de números reais de mesma dimensão de A. """
#i Vamos testar se A é TS e se a a diagonal não tem zero
# Esta parte consome muito tempo e não deve ser usada depois
n = len(A)
for i in range(n):
teste1 = A[i,i]==0
if i==n:
teste2 = False
else: teste2 = False in (A[i+1:n,i]==0)
if teste1 or teste2:
raise ValueError(" Matriz A não pode ser parametro do BackStep ")
# Esta foi só a parte de teste.
x=np.zeros(n)
k=n-1
while k>=0:
x[k] = (1/A[k,k])*(b[k]-sum([A[k,j]*x[j] for j in range(k+1,n)]))
k-=1
return x
O método direto consiste em transformar um sistema $Ax=b$ num sistema $A_1x=b_1$ sem alterar as soluções do sistema e tal que o segundo sistema seja fácil de resolver.
São as operações básicas que não alteram o conjunto solução e que iremos utilizar.
$E_{ij}$ Significa trocar as posições das equações $j$ e $i$. Na matriz de coeficientes e no vetor $b$ corresponde a trocar estas linhas.
$E_{ij,\alpha}$: isto corresponde a fazer a substituição $L_i = L_i + \alpha L_j$, ou seja, somar a equação $i$ com um múltiplo da equação $j$
$E_{i,\beta}$ para $\beta \neq 0$ correspode a multiplicação da $i$-ésima equação por $\beta$ ($=E_{ii,{\beta-1}}$).
Para o $k$ de $1$ até $n-1$ executamos:
Passo $k$: Caso $a_{kk} \neq 0$, faça de $ i=k+1$ até $ i=n $: $ m_{ik} = a_{ik}/a_{kk} $ e aplique ao sistema a operação elementar $ E_{ik,m_{ik}}$.
Ao fim do passo $n-1$ o sistema estará triangularizado. Para evitar a possibilidade de não se triangularizar o sistema, mesmo ele sendo possível e determinado, executamos a pivotação que é o seguinte procedimento: antes da execução do passo $k$ do MEG. Buscamos, na $k-$ésima coluna, a partir da $k-$ésima linha o índice $i$ tal que $|a_{ik}|$ seja o maior possível, e executamos $E_{ik}$. Com isto garantimos a triangularização de qualquer sistema que tenha solução única e evitamos perdas de algarismos significativos.
# A is type array from numpy
# remember: A[0.0] is the first element
# Elementary row operations:
def troca_linha(A,i,j):
'''Troca as linhas i e j da matriz A'''
buffer = A[i].copy()
A[i] = A[j]
A[j] = buffer
return A
def mult_linha(A,i,alfa):
'''Multiplica por alfa a linha i'''
A[i] = alfa*A[i]
return A
def subs_linha(A,i,k,alfa):
'''soma a linha i com um multiplo da linha k'''
A[i] = A[i] + alfa*A[k]
return A
# step k of gauss elimination:
def gauss_step(A,k):
'''O k-esimo passo na eliminacao de gauss'''
L=range(len(A))
for j in L[k+1:]:
m=-A[j,k]/A[k,k]
A=subs_linha(A,j,k,m)
A[j,k]=-m # Aproveitamos a matriz A pra guardar o multiplicador
return A
# triangle form, without care and pivoting
def triangle_form(A):
'''Retorna a matriz escalonada,
com os multiplicadores na parte triangular inferior'''
L=len(A)
for k in range(L-1):
p=pivo(A,k)
troca_linha(A,k,p)
gauss_step(A,k)
return A
def pivo(A,n):
''' retorna o indice de pivotacao da matriz A na coluna n,
a partir da linha n'''
tamanho = shape(A)
l=tamanho[0] # numero de linhas de A
c= tamanho[1] # numero de colunas
if (l<=n) or (c<=n):
raise ValueError("n deve ser menor que dimensao de A")
p=n # inicio o pivo como n, e vou aumentando
for k in range(n,l):
if abs(A[k,n]) > abs(A[p,n]):
p=k
return p
#Teste
# do Backstep
A=np.matrix([[1,2,3],[0,1,2],[0,0,3]])
print(A)
b=[1,2,3]
x=BackStep(A,b)
print(x)
[[1 2 3] [0 1 2] [0 0 3]] [-2. 0. 1.]
# teste da função troca_linha
A1=troca_linha(A,0,1)
A1
matrix([[0, 1, 2], [1, 2, 3], [0, 0, 3]])
A
matrix([[0, 1, 2], [1, 2, 3], [0, 0, 3]])
A2 = troca_linha(A,1,2)
A2
matrix([[0, 1, 2], [0, 0, 3], [1, 2, 3]])
A
matrix([[0, 1, 2], [0, 0, 3], [1, 2, 3]])
A3 = mult_linha(A,1,1/3)
A3
matrix([[0, 1, 2], [0, 0, 1], [1, 2, 3]])
# teste da função de pivotação
pivo(A,1)
2
#teste da função de triangularização
triangle_form(A)
matrix([[1, 2, 3], [0, 1, 2], [0, 0, 1]])