O objetivo deste documento é fazer uma breve revisão sobre da resolução de sistemas lineares utilizando as linguagens R e Python.
Um sistema linear cuja matriz dos coeficientes é uma matriz triangular é denominado um sistema triangular e sua resolução é particularmente mais fácil e de baixo custo. Segue abaixo a resolução em R de um sistema linear por meio de uma matriz triangular.
# 1. Sistema triangular superior:
# Primeiro vamos definir o tamanho da matriz: 3x3
n = 3
A <- matrix(c(2,1,2, 0,-1,4, 0,0,5), nrow = 3, ncol = 3, byrow = T)
A # matriz de coeficientes## [,1] [,2] [,3]
## [1,] 2 1 2
## [2,] 0 -1 4
## [3,] 0 0 5B <- matrix(c(5,9,10), byrow = T)
B # matriz de constantes## [,1]
## [1,] 5
## [2,] 9
## [3,] 10# Vamos criar o vetor x para resolver o sistema
x = numeric(n)
x[n] <- B[n]/A[n,n]
# Vamos fazer um laço de repetição para percorrer a matriz:
for (k in (n-1):1) {
s = 0
for (j in (k+1):n) {
s = s + A[k,j]*x[j]
}
x[k] = (B[k] - s)/A[k,k]
}
cat('A solução do sistema é:', x)## A solução do sistema é: 1 -1 2A solução do sistema foi \(x^* = (1, -1, 2)\)
Agora segue a solução em Python.
import numpy as np # pacote para a construção de matrizes
from scipy.linalg import solve_triangular # para resolução do sistema
A = np.array([[2,1,2], [0,-1,4], [0,0,5]])
A## array([[ 2, 1, 2],
## [ 0, -1, 4],
## [ 0, 0, 5]])B = np.array([5,9,10])
B## array([ 5, 9, 10])x = solve_triangular(A, B, lower=False)
print (x)## [ 1. -1. 2.]Como observado, a solução foi a mesma em ambas as linguagens.
A eliminação gaussiana, também conhecida como escalonamento, é um método para resolver sistemas lineares. Este método consiste em manipular o sistema através de determinadas operações elementares, transformando a matriz estendida do sistema em uma matriz triangular (chamada de matriz escalonada do sistema). Uma vez, triangularizado o sistema, a solução pode ser obtida via substituição regressiva. Naturalmente estas operações elementares devem preservar a solução do sistema e consistem em:
Segue abaixo a solução em R:
## 2. Eliminação de Gauss
n = 3
A <- matrix(c(3,2,4, 1,1,2, 4,3,-2), nrow = 3, ncol = 3, byrow = T)
A # matriz de coeficientes## [,1] [,2] [,3]
## [1,] 3 2 4
## [2,] 1 1 2
## [3,] 4 3 -2B <- matrix(c(1,2,3), byrow = T)
B # matriz de constantes## [,1]
## [1,] 1
## [2,] 2
## [3,] 3for (k in 1:(n-1)) {
for (i in (k+1):n) {
m = A[i,k]/A[k,k]
A[i,k] = 0
for (j in (k+1):n) {
A[i,j] = A[i,j] - m%*%A[k,j]
}
B[i] = B[i] - m%*%B[k]
}
}
# Agora vamos utilizar o código da aula anterior:
x = numeric(n)
x[n] <- B[n]/A[n,n]
for (k in (n-1):1) {
s = 0
for (j in (k+1):n) {
s = s + A[k,j]*x[j]
}
x[k] = (B[k] - s)/A[k,k]
}
print(round(x), 2)## [1] -3 5 0cat('A solução do sistema é:', round(x, 2))## A solução do sistema é: -3 5 0A solução do sistema linear foi de \(x = (-3, 5, 0)\).
Segue abaixo a solução em Python.
import numpy as np # pacote para a construção de matrizes
A = np.array([[3,2,4], [1,1,2], [4, 3, -2]])
A## array([[ 3, 2, 4],
## [ 1, 1, 2],
## [ 4, 3, -2]])B = np.array([[1], [2], [3]])
B
# Vamos criar a inversa da matriz A## array([[1],
## [2],
## [3]])invA = np.linalg.inv(A)
invA
# Agora vamos encontrar as soluções para o sistema:## array([[ 1. , -2. , 0. ],
## [-1.25 , 2.75 , 0.25 ],
## [ 0.125, 0.125, -0.125]])solucao = np.matmul(invA, B)
print('A solução para o sistema foi \n{}'.format(solucao))## A solução para o sistema foi
## [[-3.]
## [ 5.]
## [ 0.]]Como esperado, a solução para o sistema foi \(x = (-3, 5, 0)\) assim como na linguagem R.
A fatoração LU é um dos processos de fatoração mais empregados. Nesta fatoração a matriz L é triangular inferior com diagonal unitária e a matriz U é triangular superior.
Segue abaixo a solução em R.
A0 = matrix(c(3,-4,1, 1,2,2, 4,0,-3), nrow = 3, ncol = 3, byrow = TRUE)
A0## [,1] [,2] [,3]
## [1,] 3 -4 1
## [2,] 1 2 2
## [3,] 4 0 -3p0 <- matrix(c(0,0,1, 0,1,0, 1,0,0), nrow = 3, ncol = 3, byrow = TRUE)
p0## [,1] [,2] [,3]
## [1,] 0 0 1
## [2,] 0 1 0
## [3,] 1 0 0A0 <- p0 %*% A0
A0## [,1] [,2] [,3]
## [1,] 4 0 -3
## [2,] 1 2 2
## [3,] 3 -4 1m0 <- matrix(c(1,0,0, 0,1,0, 0,0,1), nrow = 3, ncol = 3, byrow = TRUE)
m0[2,1] <- -A0[2,1]/A0[1,1]
m0[3,1] <- -A0[3,1]/A0[1,1]
m0## [,1] [,2] [,3]
## [1,] 1.00 0 0
## [2,] -0.25 1 0
## [3,] -0.75 0 1A1 <- m0 %*% A0
A1## [,1] [,2] [,3]
## [1,] 4 0 -3.00
## [2,] 0 2 2.75
## [3,] 0 -4 3.25p1 <- matrix(c(1,0,0, 0,0,1, 0,1,0), nrow = 3, ncol = 3, byrow = TRUE)
p1## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 0 1
## [3,] 0 1 0A1 <- p1 %*% A1
A1## [,1] [,2] [,3]
## [1,] 4 0 -3.00
## [2,] 0 -4 3.25
## [3,] 0 2 2.75L0 <- matrix(c(0,0,0, -m0[2,1], 0,0,-m0[3,1],0,0), nrow = 3, ncol = 3, byrow = T)
L0## [,1] [,2] [,3]
## [1,] 0.00 0 0
## [2,] 0.25 0 0
## [3,] 0.75 0 0L1 <- p1 %*% L0
L1[1,1] <- L1[2,2] <- L1[3,3] <- 1
L1## [,1] [,2] [,3]
## [1,] 1.00 0 0
## [2,] 0.75 1 0
## [3,] 0.25 0 1m1 <- diag(nrow = 3, ncol = 3)
m1[3,2] <- -A1[3,2]/A1[2,2]
m1## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] 0 1.0 0
## [3,] 0 0.5 1A2 <- m1 %*% A1
round(A2, 2)## [,1] [,2] [,3]
## [1,] 4 0 -3.00
## [2,] 0 -4 3.25
## [3,] 0 0 4.38L1[3,2] <- -m1[3,2]
L <- L1
L## [,1] [,2] [,3]
## [1,] 1.00 0.0 0
## [2,] 0.75 1.0 0
## [3,] 0.25 -0.5 1p <- p1%*%p0
p## [,1] [,2] [,3]
## [1,] 0 0 1
## [2,] 1 0 0
## [3,] 0 1 0b <- c(9,3,-2)
m0_ <- m0
aux <- m0_[2,1]
m0_[2,1] <- m0_[3,1]
m0_[3,1] <- aux
m0_## [,1] [,2] [,3]
## [1,] 1.00 0 0
## [2,] -0.75 1 0
## [3,] -0.25 0 1y <- m1%*%m0_%*%p%*%b
y## [,1]
## [1,] -2.00
## [2,] 10.50
## [3,] 8.75A solução do sistema foi \(x = -2, 10.5 e 8.75\).
Segue abaixo a solução em Python.
import pprint
import scipy
import scipy.linalg
A = scipy.array([ [3,-4,1], [1,2,2], [4,0,-3] ])
A## array([[ 3, -4, 1],
## [ 1, 2, 2],
## [ 4, 0, -3]])B = scipy.array([ [9], [3], [-2] ])
B## array([[ 9],
## [ 3],
## [-2]])P, L, U = scipy.linalg.lu(A)
pprint.pprint(P)## array([[0., 1., 0.],
## [0., 0., 1.],
## [1., 0., 0.]])pprint.pprint(L)## array([[ 1. , 0. , 0. ],
## [ 0.75, 1. , 0. ],
## [ 0.25, -0.5 , 1. ]])pprint.pprint(U)## array([[ 4. , 0. , -3. ],
## [ 0. , -4. , 3.25 ],
## [ 0. , 0. , 4.375]])