7. Eliminação de Gauss¶

Até aqui, cada problema tinha uma incógnita. Muitos têm várias, ligadas entre si: as correntes nas malhas de um circuito, as quantidades de cada alimento numa dieta, a posição de um celular a partir das torres, as forças numa treliça. Quando as ligações são lineares (somas de incógnitas multiplicadas por números), o problema é um sistema linear, \(A\,x = b\).

Você resolveu sistemas de 2 ou 3 equações na escola, isolando e substituindo. Esse jeito não escala: um simulador de circuito resolve sistemas de milhares de equações, e uma previsão do tempo, de milhões. O método que escala é a eliminação de Gauss: o mesmo escalonamento da escola, organizado em laços. Esta unidade começa pelos métodos que resolvem (este capítulo e o 8) e termina, no capítulo 9, com a parte que nenhuma biblioteca faz por você: montar o sistema a partir de um problema.

Sistemas lineares e matrizes¶

🌍 Contexto: as malhas de um circuito

Num circuito com várias malhas (os "anéis" do desenho), cada malha tem uma corrente. A lei de Kirchhoff diz que, em cada malha, a soma das tensões nos resistores é igual à tensão das fontes. Cada malha dá uma equação, e cada equação envolve a corrente dela e a das vizinhas (os resistores que elas dividem). Três malhas dão um sistema de três equações e três incógnitas. Como se chega a estas equações a partir do desenho é assunto do capítulo 9.

Um sistema como

\[ \begin{aligned} 15\,i_1 - 5\,i_2 \phantom{{}- 10\,i_3} &= 10\\ -5\,i_1 + 20\,i_2 - 10\,i_3 &= 0\\ -10\,i_2 + 25\,i_3 &= 5 \end{aligned} \]

se escreve como \(A\,x = b\): a matriz \(A\) guarda os coeficientes, \(x\) as incógnitas e \(b\) o lado direito. Com o numpy, a matriz é um array de duas dimensões.

🧰 Matriz: o array de duas dimensões

np.array recebe uma lista de linhas, e cada linha é uma lista de números. A[i, j] é o elemento da linha i e coluna j (contando do zero), A[i] é a linha i inteira e len(A) é o número de linhas:

import numpy as np

A = np.array([[2.0, 1.0, -1.0],
              [4.0, 5.0, 3.0],
              [1.0, 0.0, 2.0]])   # uma lista de linhas
print(A)
print(A[1, 2])      # linha 1, coluna 2 (contando do zero)
print(A[0])         # a linha 0 inteira
print(len(A))       # quantas linhas

[[ 2.  1. -1.]
 [ 4.  5.  3.]
 [ 1.  0.  2.]]
3.0
[ 2.  1. -1.]
3

🧰 Comando novo: .copy()

A eliminação de Gauss modifica a matriz. Para não estragar a original, o método trabalha numa cópia. Mas atenção: b = a não copia um array; só dá a ele um segundo nome. Para copiar de verdade, use a.copy():

import numpy as np

a = np.array([1.0, 2.0, 3.0])
b = a               # NAO copia: b e outro nome para o MESMO array
b[0] = 99
print(a)            # a mudou junto!

c = a.copy()        # copia de verdade
c[1] = -5
print(a, c)         # agora a nao muda

[99.  2.  3.]
[99.  2.  3.] [99. -5.  3.]

No quadro: a eliminação de Gauss¶

🧑‍🏫 No quadro — eliminação e substituição

1. Eliminação. Use a equação 0 para zerar a primeira incógnita em todas as equações de baixo: da linha \(i\), subtraia a linha 0 multiplicada por \(\text{fator} = A_{i0}/A_{00}\). O número \(A_{00}\), que divide, é o pivô.

2. Repita com a equação 1 para zerar a segunda incógnita nas linhas de baixo, e assim por diante. No fim, a matriz é triangular superior: só há números na diagonal e acima dela.

3. Substituição. A última equação tem uma incógnita só: \(x_{n-1} = b_{n-1}/A_{n-1,n-1}\). A penúltima tem duas, uma já conhecida. Subindo, cada equação revela uma incógnita:

\[ x_i = \frac{b_i - \sum_{j > i} A_{ij}\,x_j}{A_{ii}}. \]

4. É o escalonamento da escola. A diferença é que, organizado assim, ele vira três laços for que funcionam para 3 ou para 3 000 equações.

A substituição anda de trás para frente: começa na última linha e sobe.

🧰 Python: range de trás para frente

range(inicio, fim, passo) com passo negativo conta para trás. Como sempre, o fim fica de fora: range(n - 1, -1, -1) dá n-1, n-2, ..., 1, 0. Com n = 3, dá 2, 1, 0.

# Circuito com tres malhas: as leis de Kirchhoff dao um sistema linear nas
# correntes de malha i1, i2, i3 (em ampere).
#   15 i1 -  5 i2 -  0 i3 = 10
#   -5 i1 + 20 i2 - 10 i3 =  0
#    0 i1 - 10 i2 + 25 i3 =  5
import numpy as np

A = np.array([[15.0, -5.0, 0.0],
              [-5.0, 20.0, -10.0],
              [0.0, -10.0, 25.0]])
b = np.array([10.0, 0.0, 5.0])

n = len(b)
A = A.copy()        # a eliminacao muda A e b: trabalhamos em copias
b = b.copy()

# 1) eliminacao: zera tudo abaixo da diagonal, coluna por coluna
for k in range(n - 1):
    for i in range(k + 1, n):
        fator = A[i, k] / A[k, k]
        for j in range(k, n):
            A[i, j] = A[i, j] - fator * A[k, j]
        b[i] = b[i] - fator * b[k]
print("matriz triangular:")
print(np.round(A, 4))

# 2) substituicao de tras para frente
x = np.zeros(n)
for i in range(n - 1, -1, -1):
    soma = b[i]
    for j in range(i + 1, n):
        soma = soma - A[i, j] * x[j]
    x[i] = soma / A[i, i]
print("correntes (A):", np.round(x, 6))

matriz triangular:
[[ 15.      -5.       0.    ]
 [  0.      18.3333 -10.    ]
 [  0.       0.      19.5455]]
correntes (A): [0.790698 0.372093 0.348837]

A matriz triangular mostra a eliminação trabalhando: abaixo da diagonal, só zeros. As correntes saem na substituição.

Pivotamento¶

A eliminação divide pelo pivô \(A_{kk}\). Se ele for zero, a conta quebra. Se for muito pequeno, ela não quebra, mas erra.

# Quebre o metodo: um pivo (A[k, k]) igual a zero, ou muito pequeno.
import numpy as np


def gauss_sem_pivo(A, b):
    A = A.copy()
    b = b.copy()
    n = len(b)
    for k in range(n - 1):
        for i in range(k + 1, n):
            fator = A[i, k] / A[k, k]
            for j in range(k, n):
                A[i, j] = A[i, j] - fator * A[k, j]
            b[i] = b[i] - fator * b[k]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        soma = b[i]
        for j in range(i + 1, n):
            soma = soma - A[i, j] * x[j]
        x[i] = soma / A[i, i]
    return x


# x + y = 2 e 1e-17 x + y = 1: a resposta e x = 1, y = 1 (quase exatos)
A = np.array([[1e-17, 1.0],
              [1.0, 1.0]])
b = np.array([1.0, 2.0])
print("pivo minusculo em cima:  ", gauss_sem_pivo(A, b))

# as mesmas equacoes, com as linhas trocadas
A_trocada = np.array([[1.0, 1.0],
                      [1e-17, 1.0]])
b_trocado = np.array([2.0, 1.0])
print("linhas trocadas (pivo 1):", gauss_sem_pivo(A_trocada, b_trocado))

pivo minusculo em cima:   [0. 1.]
linhas trocadas (pivo 1): [1. 1.]

💥 Quebre o método

Com o pivô \(10^{-17}\) em cima, o fator é \(10^{17}\), e a linha de baixo vira \(1 - 10^{17} \approx -10^{17}\): o "1" original some no arredondamento (capítulo 1). A resposta sai \(x = 0\), errada, sem nenhum aviso. Trocando as linhas de lugar, o pivô passa a ser 1 e a resposta volta a \(x = 1\).

A defesa se chama pivotamento parcial: antes de eliminar cada coluna, troque a linha do pivô pela linha que tem o maior valor absoluto naquela coluna. Todo programa sério de sistemas lineares faz isso.

Confira com a biblioteca¶

🧰 Comando novo: np.linalg.solve e o produto @

np.linalg.solve(A, b) resolve \(A\,x = b\) com eliminação de Gauss com pivotamento, escrita em código otimizado. E A @ x é o produto da matriz pela coluna: a conta que monta o lado esquerdo do sistema. Substituir a resposta e conferir que A @ x dá b é o jeito de verificar qualquer sistema.

# Confira com a biblioteca, e confira de novo substituindo a resposta.
import numpy as np

A = np.array([[15.0, -5.0, 0.0],
              [-5.0, 20.0, -10.0],
              [0.0, -10.0, 25.0]])
b = np.array([10.0, 0.0, 5.0])

x = np.linalg.solve(A, b)
print("np.linalg.solve:", x)
print("A @ x =", A @ x)            # tem de dar b de volta
print("residuo:", b - A @ x)

np.linalg.solve: [0.79069767 0.37209302 0.34883721]
A @ x = [ 1.00000000e+01 -1.11022302e-16  5.00000000e+00]
residuo: [0.00000000e+00 1.11022302e-16 0.00000000e+00]

O resíduo \(b - A\,x\) é da ordem de \(10^{-16}\): arredondamento, e nada mais.

Sistemas mal condicionados¶

Às vezes o problema não está no método, mas no sistema. Duas equações de retas quase paralelas se cruzam num ponto que muda de lugar com qualquer tremida:

# Sistema mal condicionado: duas retas quase paralelas.
#   x +      y = 2
#   x + 1.0001 y = 2.0001        -> x = 1, y = 1
import numpy as np

A = np.array([[1.0, 1.0],
              [1.0, 1.0001]])
b = np.array([2.0, 2.0001])
print("resposta:            ", np.linalg.solve(A, b))

b_medido = np.array([2.0, 2.0002])   # um erro de 0,005 % num dado
print("com b um pouco errado:", np.linalg.solve(A, b_medido))

resposta:             [1. 1.]
com b um pouco errado: [0. 2.]

Um erro de 0,005 % num único dado mudou a resposta de \((1, 1)\) para \((0, 2)\). O sistema é mal condicionado: pequenas mudanças nos dados viram grandes mudanças na resposta. Nenhum método conserta isso, porque a resposta muda de verdade. O remédio é perceber (perturbe os dados de propósito e veja o quanto a resposta muda) e, se der, reformular o problema ou medir melhor.

💡 Pense assim: achar um cruzamento no mapa

Duas ruas perpendiculares se cruzam num ponto bem definido: errar um pouco o traçado de uma delas mexe pouco no cruzamento. Duas ruas quase paralelas se cruzam muito longe, e um erro pequeno no ângulo de uma delas leva o cruzamento para outro bairro.

Mesmo método, outra área¶

🌍 Contexto: planejar uma refeição

Nutricionistas montam refeições com metas de proteína, carboidrato e gordura. Cada alimento tem uma composição conhecida (gramas de cada nutriente por grama de alimento). Com três alimentos e três metas, a pergunta "quantos gramas de cada?" é um sistema 3 × 3.

# Mesmo metodo, outra area: nutricao. Quantos gramas de arroz, feijao e frango
# dao exatamente 50,6 g de proteina, 77 g de carboidrato e 4,15 g de gordura?
# Composicao por grama (proteina, carboidrato, gordura):
#   arroz  0.025  0.28   0.002
#   feijao 0.048  0.14   0.005
#   frango 0.32   0.0    0.025
import numpy as np

A = np.array([[0.025, 0.048, 0.32],     # proteina de cada alimento
              [0.28, 0.14, 0.0],        # carboidrato
              [0.002, 0.005, 0.025]])   # gordura
b = np.array([50.6, 77.0, 4.15])

gramas = np.linalg.solve(A, b)
print("arroz, feijao, frango (g):", np.round(gramas, 1))
print("confere (proteina, carb., gordura):", np.round(A @ gramas, 3))

arroz, feijao, frango (g): [200. 150. 120.]
confere (proteina, carb., gordura): [50.6  77.    4.15]

As linhas da matriz são os nutrientes, as colunas são os alimentos, e a resposta é a refeição: 200 g de arroz, 150 g de feijão e 120 g de frango. Com metas escolhidas ao acaso, o sistema pode dar quantidades negativas: matematicamente certas, fisicamente impossíveis. Todo resultado precisa passar pela pergunta "isso faz sentido?".

Erros comuns deste capítulo¶

Sintoma	Causa provável
a matriz original apareceu modificada depois da função	faltou o `.copy()`: a eliminação mexeu na matriz de quem chamou
`ZeroDivisionError` ou `inf` na eliminação	pivô zero: troque as linhas (pivotamento)
resposta errada, sem erro nenhum	pivô minúsculo, ou coeficientes na ordem errada: confira `A @ x` contra `b`
`IndexError` na substituição	o `range` de trás para frente: `range(n - 1, -1, -1)`
`LinAlgError: Singular matrix`	o sistema não tem solução única: duas equações dizem a mesma coisa
a resposta muda muito com dados quase iguais	sistema mal condicionado: o problema é dos dados, não do método

Onde isso é usado de verdade¶

Simuladores de circuitos (SPICE e seus descendentes): cada instante de simulação é um sistema linear com uma equação por nó do circuito.
Engenharia de estruturas: o método dos elementos finitos transforma uma ponte ou um prédio num sistema com milhões de equações.
GPS e localização: a posição do receptor sai de sistemas lineares montados a partir das distâncias aos satélites (ou às antenas de celular).
Computação gráfica: cada transformação de uma cena 3D (rotação, projeção) é uma multiplicação por matriz, e desfazê-las é resolver sistemas.
Economia: o modelo de insumo-produto de Leontief (que rendeu um prêmio Nobel) é um sistema linear com um setor da economia por equação.

Resumo¶

Várias incógnitas ligadas linearmente formam um sistema \(A\,x = b\). No numpy, \(A\) é um array de duas dimensões, e A[i, j] é um elemento.
Eliminação de Gauss: zerar abaixo da diagonal, coluna por coluna; depois, substituição de trás para frente.
Pivô zero ou minúsculo estraga a eliminação: troque linhas (pivotamento).
np.linalg.solve confere; A @ x deve devolver b.
Sistema mal condicionado amplifica erros dos dados; não há método que conserte.
Toda resposta precisa fazer sentido no problema (nada de gramas negativos).

Para praticar¶

A Lista 07 começa com uma eliminação à mão (✏️) e segue com um circuito, o balanceamento de uma reação química, uma treliça e a localização de um celular.