Code Session 5: Sistemas Lineares#

linalg.solve#

A função solve é o método mais simples disponibilizado pelos módulos numpy e scipy para resolver sistemas matriciais lineares. Como a função pertence ao escopo da Álgebra Linear, ela está localizada em submódulo chamado linalg. solve calculará a solução exata do sistema como um método direto se a matriz do sistema for determinada (quadrada e sem colunas linearmente dependentes). Se a matriz for singular, o método retorna um erro. Se for de posto deficiente, o método resolve o sistema linear usando um algoritmo de mínimos quadrados.

Os argumentos de entrada obrigatórios desta função são:

  1. a matriz A dos coeficientes

  2. o vetor independente b

O argumento de saída é:

  • x: o vetor-solução do sistema.

from numpy.linalg import solve

Problema 1#

Uma rede elétrica contém 3 loops. Ao aplicar a lei de Kirchoff a cada loop, o cientista Hely Johnson obteve as seguintes equações para as correntes \(i_1\), \(i_2\) e \(i_3\) em cada loop:

\[\begin{split}\begin{cases} (50 + R)i_1 - Ri_2 - 30i_3 &= 0 \\ - Ri_1 + (65 + R)i_2 - 15i_3 &= 0 \\ -30i_1 - 15i_2 + 45i_3 &= 120 \end{cases}\end{split}\]

Ajude Hely Johnson em seus experimentos e calcule as correntes para os valores de resistência \(R = \{ 5 \Omega, 10 \Omega, 20 \Omega \}\).

Resolução#

Podemos começar definindo uma lista para armazenar os valores das resistências de teste:

R = [5.,10.,15.]

Em seguida, escreveremos a matriz dos coeficientes e o vetor independente (lado direito). Note que precisamos de uma “lista de listas”, ou melhor, um “array de arrays”, onde cada array está associado à uma linha da matriz.

Todavia, vamos definir uma função para montar o sistema em função do valor de \(R\) e resolvê-lo.

# montagem do sistema
def resolve_sistema(R):
    A = np.array([ [50+R,-R,-30],[-R,65+R,-15],[-30,-15,45] ])
    b = np.array([0,0,120])
    x = solve(A,b)
    return x

Além disso, usaremos um laço for para calcularmos todas as respostas necessárias de uma só vez e organizaremos os resultados em dicionário (objeto dict):

# salva soluções agrupadas em um dicionário
sols = {}
for r in R:
    sols[r] = resolve_sistema(r)
    print('Para o valor de resistência R = {0:g} ohms: i1 = {1:g} A, i2 = {2:g} A, i3 = {3:g} A'.format(r, sols[r][0], sols[r][1], sols[r][2]))

Para o valor de resistência R = 5 ohms: i1 = 2.82927 A, i2 = 1.26829 A, i3 = 4.97561 A
Para o valor de resistência R = 10 ohms: i1 = 2.66667 A, i2 = 1.33333 A, i3 = 4.88889 A
Para o valor de resistência R = 15 ohms: i1 = 2.54545 A, i2 = 1.38182 A, i3 = 4.82424 A

Tarefa#

Retorne ao notebook da Aula 09 e use o conteúdo da aula para fazer uma verificação das soluções encontradas para este problema em cada caso. Use a função linalg.cond para calcular o número de condicionamento da matriz do sistema em cada caso. (Sugestão: vide numpy.allclose)

linalg.cholesky#

Assim como solve, a função cholesky está disponível tanto via numpy como scipy para determinar a decomposição de Cholesky de uma matriz simétrica e positiva definida.

Na prática, não é recomendável verificar se uma matriz é positiva definida através dos critérios teóricos. A função cholesky não realiza essa checagem. Portanto, é importante que, pelo menos, se saiba que a matriz é simétrica. Para testar se ela atende a propriedade de definição positiva, a abordagem mais direta é usar cholesky e verificar se ela retorna uma decomposição de Cholesky. Se não for o caso, um erro será lançado.

O argumento de entrada desta função é:

  • A: matriz dos coeficientes

O argumento de saída é:

  • L: o fator de Cholesky

Como importá-la?

from numpy.linalg import cholesky

from numpy.linalg import cholesky

Problema 2#

Calcule o fator de Cholesky para a matriz \(A\) do Problema 1 para \(R = 5\).

R = 5.
B = np.array([ [50+R,-R,-30],[-R,65+R,-15],[-30,-15,45] ])

L = cholesky(B)
L

array([[ 7.41619849,  0.        ,  0.        ],
       [-0.67419986,  8.33939174,  0.        ],
       [-4.04519917, -2.12572731,  4.91097211]])

Podemos verificar a decomposição multiplicando a matriz triangular (fator de Cholesky) pela sua transposta.

np.matmul(L,L.T)

array([[ 55.,  -5., -30.],
       [ -5.,  70., -15.],
       [-30., -15.,  45.]])

B

array([[ 55.,  -5., -30.],
       [ -5.,  70., -15.],
       [-30., -15.,  45.]])

Entretanto, não é verdade que

# por que a igualdade falha 
# para algumas entradas?
B == np.matmul(L,L.T)

array([[ True, False,  True],
       [False,  True, False],
       [ True, False,  True]])

Problema 3#

Verificar se uma matriz simétrica é positiva definida.

# cria matriz simétrica 
C = np.tril(B) - 60
C = np.tril(C) + np.tril(C,-1).T
C

array([[ -5., -65., -90.],
       [-65.,  10., -75.],
       [-90., -75., -15.]])

# erro! 
# matriz não é PD
#cholesky(C)

# cria outra matriz simétrica
D = np.tril(B) + 1
D = np.tril(D) + np.tril(D,-1).T
D

array([[ 56.,  -4., -29.],
       [ -4.,  71., -14.],
       [-29., -14.,  46.]])

# matriz é PD
cholesky(D)

array([[ 7.48331477,  0.        ,  0.        ],
       [-0.53452248,  8.40917866,  0.        ],
       [-3.87528801, -1.91117697,  5.22776678]])

Problema 4#

Resolva o Problema 1 para \(R = 10\) usando a fatoração de Cholesky.

R = 10.
D = np.array([ [50+R,-R,-30],[-R,65+R,-15],[-30,-15,45] ])
b = np.array([0,0,120])

# fator 
L = cholesky(D)

# passo 1
# Ly = b
y = solve(L,b)

# passo 2
# L^T x = y
x = solve(L.T,y)

# solução
x

array([2.66666667, 1.33333333, 4.88888889])

Compare a solução via Cholesky com a do Problema 1:

x, sols[10]

(array([2.66666667, 1.33333333, 4.88888889]),
 array([2.66666667, 1.33333333, 4.88888889]))

Mais uma vez, note que:

# A comparação falha.
# Por quê?
x == sols[10]

array([False,  True,  True])