Pular para conteúdo

Lista 08 · Jacobi e Gauss-Seidel

Acompanha o capítulo 8. Métodos iterativos para sistemas: chutar, melhorar e parar. Confira a dominância diagonal antes de iterar.

Abrir no Colab

Baixar o caderno

Como esta lista funciona

Abra o caderno no Colab, escreva a sua solução no lugar do # TODO e rode a célula de teste logo abaixo. Ela mesma diz o que passou e o que falhou — você não precisa esperar correção de ninguém.

O código de cada exercício também está aqui embaixo, caso você prefira colar num editor local. Os testes comparam números com tolerância.

Exercícios marcados com ✏️ são à mão, no papel, com 4 casas decimais — como na parte em papel das avaliações.

Exercício 01

✏️ À mão. Considere o sistema

\[ 10x + y = 12, \qquad x + 5y = -4. \]

a) Ele tem dominância diagonal? Justifique.

b) Faça 2 iterações de Jacobi a partir de \((0, 0)\).

c) Faça 2 iterações de Gauss-Seidel a partir de \((0, 0)\).

d) A solução exata é \(x = 64/49 = 1{,}3061\), \(y = -52/49 = -1{,}0612\). Qual dos dois chegou mais perto?

Exercício 02

Escreva diagonal_dominante(A), que devolve True se, em todas as linhas, o valor absoluto do elemento da diagonal é maior que a soma dos valores absolutos dos outros elementos da linha.

def diagonal_dominante(A):
    """True se, em toda linha, |A[i, i]| > soma dos |A[i, j]| com j != i."""
    n = len(A)
    # TODO: para cada linha, some os abs fora da diagonal e compare com abs(A[i, i])
    pass


# --- testes (não altere) ---
import numpy as np


assert diagonal_dominante(np.array([[3.0, -0.1, -0.2], [0.1, 7.0, -0.3], [0.3, -0.2, 10.0]])) is True
assert diagonal_dominante(np.array([[1.0, 3.0], [4.0, 1.0]])) is False
assert diagonal_dominante(np.array([[4.0, 1.0], [1.0, 3.0]])) is True
assert diagonal_dominante(np.array([[2.0, 2.0], [1.0, 5.0]])) is False, "igual não basta: tem de ser maior"
print("Exercício 02: todos os testes passaram!")

Exercício 03

Escreva jacobi(A, b, tol, max_iter), o método de Jacobi para qualquer sistema: a cada iteração, calcule todos os valores novos num array novo, usando só os valores antigos de x; depois compare os dois (maior \(\varepsilon_a\) em %) e troque x = novo.

import numpy as np


def jacobi(A, b, tol, max_iter):
    """Solução de A x = b por Jacobi a partir de zeros; para quando o maior eps_a (%) < tol."""
    n = len(b)
    x = np.zeros(n)
    # TODO: em cada iteração: novo = np.zeros(n); preencha novo[i] com a fórmula
    # (usando x, o antigo); calcule o maior eps_a; x = novo; break se < tol
    return x


# --- testes (não altere) ---
import numpy as np


A = np.array([[3.0, -0.1, -0.2], [0.1, 7.0, -0.3], [0.3, -0.2, 10.0]])
b = np.array([7.85, -19.3, 71.4])
x = jacobi(A, b, 1e-8, 100)
for i in range(3):
    assert abs(x[i] - [3, -2.5, 7][i]) < 1e-8, f"veio {x}"
x = jacobi(A, b, 1e-8, 1)
assert abs(x[1] - (-19.3 / 7)) < 1e-12, "na primeira iteração, x2 = -19.3/7 (Jacobi usa os antigos, que são zero)"
print("Exercício 03: todos os testes passaram!")

Exercício 04

Escreva gauss_seidel(A, b, tol, max_iter): igual ao Jacobi do ex03, mas cada valor novo vai direto para x[i], e já é usado nas linhas seguintes da mesma iteração.

import numpy as np


def gauss_seidel(A, b, tol, max_iter):
    """Solução de A x = b por Gauss-Seidel a partir de zeros."""
    n = len(b)
    x = np.zeros(n)
    # TODO: como o Jacobi, mas grave cada valor novo direto em x[i]
    return x


# --- testes (não altere) ---
import numpy as np


A = np.array([[3.0, -0.1, -0.2], [0.1, 7.0, -0.3], [0.3, -0.2, 10.0]])
b = np.array([7.85, -19.3, 71.4])
x = gauss_seidel(A, b, 1e-8, 100)
for i in range(3):
    assert abs(x[i] - [3, -2.5, 7][i]) < 1e-8, f"veio {x}"
x = gauss_seidel(A, b, 1e-8, 1)
assert abs(x[1] - (-19.3 - 0.1 * 7.85 / 3) / 7) < 1e-12, "na primeira iteração, x2 já usa o x1 novo"
print("Exercício 04: todos os testes passaram!")

Exercício 05

Consertando a ordem. Escreva reordena(A, b), que monta uma nova ordem das equações: para cada coluna \(j\), a linha que tem o maior valor absoluto nessa coluna vai para a posição \(j\). Devolva a tupla (A_novo, b_novo), sem modificar A e b.

(Supõe-se que cada coluna tem o seu maior valor numa linha diferente.)

A = np.array([[1.0, 3.0], [4.0, 1.0]])
b = np.array([7.0, 6.0])
reordena(A, b)    # -> (array([[4., 1.], [1., 3.]]), array([6., 7.]))

def reordena(A, b):
    """Troca as linhas para que, em cada coluna j, a linha com o maior |A[i, j]|
    fique na posição j. Devolve (A_novo, b_novo)."""
    n = len(b)
    A_novo = A.copy()
    b_novo = b.copy()
    # TODO: para cada coluna j, ache a linha com o maior abs(A[i, j]) e copie
    # essa linha (de A e de b) para a posição j de A_novo e b_novo
    return (A_novo, b_novo)


# --- testes (não altere) ---
import numpy as np


A = np.array([[1.0, 3.0], [4.0, 1.0]])
b = np.array([7.0, 6.0])
A2, b2 = reordena(A, b)
assert A2[0, 0] == 4.0 and A2[1, 1] == 3.0 and b2[0] == 6.0 and b2[1] == 7.0, f"veio {A2}, {b2}"
assert A[0, 0] == 1.0, "não modifique A"
A = np.array([[1.0, 1.0, 8.0], [9.0, 1.0, 2.0], [1.0, 7.0, 1.0]])
b = np.array([10.0, 12.0, 9.0])
A2, b2 = reordena(A, b)
assert A2[0, 0] == 9.0 and A2[1, 1] == 7.0 and A2[2, 2] == 8.0
assert b2[0] == 12.0 and b2[1] == 9.0 and b2[2] == 10.0
print("Exercício 05: todos os testes passaram!")

Exercício 06

Engenharia térmica — a placa aquecida. Escreva temperatura_centro(topo, base, esquerda, direita), que monta uma placa \(7 \times 7\) (np.zeros((7, 7))), fixa as quatro bordas nas temperaturas dadas, aplica Gauss-Seidel nos pontos de dentro (cada um vira a média dos quatro vizinhos) até a maior mudança ficar abaixo de \(10^{-9}\), e devolve a temperatura do ponto central, T[3, 3].

(Nos cantos, onde duas bordas se encontram, vale a que for escrita por último — eles não entram em nenhuma conta.)

Teste com as bordas (100, 0, 0, 0): o centro deve dar a média das quatro bordas, 25. Isso vale para quaisquer temperaturas de borda?

import numpy as np


def temperatura_centro(topo, base, esquerda, direita):
    """Temperatura no centro de uma placa 7 x 7 com as quatro bordas fixas."""
    n = 7
    T = np.zeros((n, n))
    # TODO: fixe as bordas (linha 0, linha n-1, coluna 0, coluna n-1); depois
    # Gauss-Seidel nos pontos de dentro até a maior mudança < 1e-9
    pass


# --- testes (não altere) ---
assert abs(temperatura_centro(100, 0, 0, 0) - 25) < 1e-6
assert abs(temperatura_centro(80, 20, 40, 60) - 50) < 1e-6
assert abs(temperatura_centro(30, 30, 30, 30) - 30) < 1e-6, "bordas iguais: a placa toda fica igual"
print("Exercício 06: todos os testes passaram!")

Exercício 07

Dominância forte, convergência rápida. Escreva iteracoes_gauss_seidel(A, b, tol), que roda Gauss-Seidel a partir de zeros e devolve quantas iterações foram necessárias para o maior \(\varepsilon_a\) (%) ficar abaixo de tol (no máximo 10 000).

Compare os sistemas com matrizes \(\begin{bmatrix}10 & 1\\ 1 & 10\end{bmatrix}\) e \(\begin{bmatrix}2 & 1{,}9\\ 1{,}9 & 2\end{bmatrix}\), os dois com \(b = (1, 2)\) e tol = 1e-6. Os dois têm dominância diagonal. Qual converge mais depressa, e por quê?

import numpy as np


def iteracoes_gauss_seidel(A, b, tol):
    """Número de iterações de Gauss-Seidel (a partir de zeros) até o maior
    eps_a (%) ficar abaixo de tol (no máximo 10000)."""
    n = len(b)
    x = np.zeros(n)
    # TODO: for it in range(1, 10001): uma iteração de Gauss-Seidel; se o maior
    # eps_a < tol, return it. Depois do laço, return 10000
    pass


# --- testes (não altere) ---
import numpy as np


b = np.array([1.0, 2.0])
forte = iteracoes_gauss_seidel(np.array([[10.0, 1.0], [1.0, 10.0]]), b, 1e-6)
fraca = iteracoes_gauss_seidel(np.array([[2.0, 1.9], [1.9, 2.0]]), b, 1e-6)
assert forte == 6, f"dominância forte: veio {forte}"
assert fraca == 160, f"dominância fraca: veio {fraca}"
print("Exercício 07: todos os testes passaram!")