Code Session 4: solve
Contents
Code Session 4: solve#
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:
a matriz
Ados coeficienteso vetor independente
b
O argumento de saída é:
x: o vetor-solução do sistema.
from numpy.linalg import solve
import numpy as np
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:
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.,20.]
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 np.linspace(1,50):
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 = 1 ohms: i1 = 3.00565 A, i2 = 1.19774 A, i3 = 5.06968 A
Para o valor de resistência R = 2 ohms: i1 = 2.95652 A, i2 = 1.21739 A, i3 = 5.04348 A
Para o valor de resistência R = 3 ohms: i1 = 2.91099 A, i2 = 1.2356 A, i3 = 5.0192 A
Para o valor de resistência R = 4 ohms: i1 = 2.86869 A, i2 = 1.25253 A, i3 = 4.99663 A
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 = 6 ohms: i1 = 2.79245 A, i2 = 1.28302 A, i3 = 4.95597 A
Para o valor de resistência R = 7 ohms: i1 = 2.75799 A, i2 = 1.2968 A, i3 = 4.9376 A
Para o valor de resistência R = 8 ohms: i1 = 2.72566 A, i2 = 1.30973 A, i3 = 4.92035 A
Para o valor de resistência R = 9 ohms: i1 = 2.69528 A, i2 = 1.32189 A, i3 = 4.90415 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 = 11 ohms: i1 = 2.63968 A, i2 = 1.34413 A, i3 = 4.87449 A
Para o valor de resistência R = 12 ohms: i1 = 2.61417 A, i2 = 1.35433 A, i3 = 4.86089 A
Para o valor de resistência R = 13 ohms: i1 = 2.59004 A, i2 = 1.36398 A, i3 = 4.84802 A
Para o valor de resistência R = 14 ohms: i1 = 2.56716 A, i2 = 1.37313 A, i3 = 4.83582 A
Para o valor de resistência R = 15 ohms: i1 = 2.54545 A, i2 = 1.38182 A, i3 = 4.82424 A
Para o valor de resistência R = 16 ohms: i1 = 2.52482 A, i2 = 1.39007 A, i3 = 4.81324 A
Para o valor de resistência R = 17 ohms: i1 = 2.50519 A, i2 = 1.39792 A, i3 = 4.80277 A
Para o valor de resistência R = 18 ohms: i1 = 2.48649 A, i2 = 1.40541 A, i3 = 4.79279 A
Para o valor de resistência R = 19 ohms: i1 = 2.46865 A, i2 = 1.41254 A, i3 = 4.78328 A
Para o valor de resistência R = 20 ohms: i1 = 2.45161 A, i2 = 1.41935 A, i3 = 4.77419 A
Para o valor de resistência R = 21 ohms: i1 = 2.43533 A, i2 = 1.42587 A, i3 = 4.76551 A
Para o valor de resistência R = 22 ohms: i1 = 2.41975 A, i2 = 1.4321 A, i3 = 4.7572 A
Para o valor de resistência R = 23 ohms: i1 = 2.40483 A, i2 = 1.43807 A, i3 = 4.74924 A
Para o valor de resistência R = 24 ohms: i1 = 2.39053 A, i2 = 1.44379 A, i3 = 4.74162 A
Para o valor de resistência R = 25 ohms: i1 = 2.37681 A, i2 = 1.44928 A, i3 = 4.7343 A
Para o valor de resistência R = 26 ohms: i1 = 2.36364 A, i2 = 1.45455 A, i3 = 4.72727 A
Para o valor de resistência R = 27 ohms: i1 = 2.35097 A, i2 = 1.45961 A, i3 = 4.72052 A
Para o valor de resistência R = 28 ohms: i1 = 2.3388 A, i2 = 1.46448 A, i3 = 4.71403 A
Para o valor de resistência R = 29 ohms: i1 = 2.32708 A, i2 = 1.46917 A, i3 = 4.70777 A
Para o valor de resistência R = 30 ohms: i1 = 2.31579 A, i2 = 1.47368 A, i3 = 4.70175 A
Para o valor de resistência R = 31 ohms: i1 = 2.30491 A, i2 = 1.47804 A, i3 = 4.69595 A
Para o valor de resistência R = 32 ohms: i1 = 2.29442 A, i2 = 1.48223 A, i3 = 4.69036 A
Para o valor de resistência R = 33 ohms: i1 = 2.28429 A, i2 = 1.48628 A, i3 = 4.68495 A
Para o valor de resistência R = 34 ohms: i1 = 2.27451 A, i2 = 1.4902 A, i3 = 4.67974 A
Para o valor de resistência R = 35 ohms: i1 = 2.26506 A, i2 = 1.49398 A, i3 = 4.6747 A
Para o valor de resistência R = 36 ohms: i1 = 2.25592 A, i2 = 1.49763 A, i3 = 4.66983 A
Para o valor de resistência R = 37 ohms: i1 = 2.24709 A, i2 = 1.50117 A, i3 = 4.66511 A
Para o valor de resistência R = 38 ohms: i1 = 2.23853 A, i2 = 1.50459 A, i3 = 4.66055 A
Para o valor de resistência R = 39 ohms: i1 = 2.23025 A, i2 = 1.5079 A, i3 = 4.65613 A
Para o valor de resistência R = 40 ohms: i1 = 2.22222 A, i2 = 1.51111 A, i3 = 4.65185 A
Para o valor de resistência R = 41 ohms: i1 = 2.21444 A, i2 = 1.51422 A, i3 = 4.6477 A
Para o valor de resistência R = 42 ohms: i1 = 2.2069 A, i2 = 1.51724 A, i3 = 4.64368 A
Para o valor de resistência R = 43 ohms: i1 = 2.19958 A, i2 = 1.52017 A, i3 = 4.63977 A
Para o valor de resistência R = 44 ohms: i1 = 2.19247 A, i2 = 1.52301 A, i3 = 4.63598 A
Para o valor de resistência R = 45 ohms: i1 = 2.18557 A, i2 = 1.52577 A, i3 = 4.6323 A
Para o valor de resistência R = 46 ohms: i1 = 2.17886 A, i2 = 1.52846 A, i3 = 4.62873 A
Para o valor de resistência R = 47 ohms: i1 = 2.17234 A, i2 = 1.53106 A, i3 = 4.62525 A
Para o valor de resistência R = 48 ohms: i1 = 2.16601 A, i2 = 1.5336 A, i3 = 4.62187 A
Para o valor de resistência R = 49 ohms: i1 = 2.15984 A, i2 = 1.53606 A, i3 = 4.61858 A
Para o valor de resistência R = 50 ohms: i1 = 2.15385 A, i2 = 1.53846 A, i3 = 4.61538 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, False]])
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([ True, False, False])
A = np.array([[2, 6, 0],[3, -7, 1],[5, 2, 7]])
b = np.array([0, -1, 2])
x = solve(A,b)
x
array([-0.27272727, 0.09090909, 0.45454545])
e = np.matmul(A,x) - b
np.sqrt(np.dot(e,e))
4.440892098500626e-16