Um sistema linear é um conjunto de duas ou mais equações lineares envolvendo as mesmas variáveis.
Por exemplo:
$${\begin{alignedat}{7}3x&&\;+\;&&2y&&\;-\;&&z&&\;=\;&&1&\\2x&&\;-\;&&2y&&\;+\;&&4z&&\;=\;&&-2&\\-x&&\;+\;&&{\tfrac {1}{2}}y&&\;-\;&&z&&\;=\;&&0&\end{alignedat}}$$A solução para este sistema é:
$$ {\begin{alignedat}{2}x&\,=\,&1\\y&\,=\,&-2\\z&\,=\,&-2\end{alignedat}} $$que satisfaz a todas as equações simultaneamente. A algebra linear é o campo da matemática que estuda estes sistemas.
Na sua forma geral temos:
$$ {\displaystyle {\begin{aligned}a_{11}x_{1}&+a_{12}x_{2}+\cdots +a_{1n}x_{n}=b_{1}\\a_{21}x_{1}&+a_{22}x_{2}+\cdots +a_{2n}x_{n}=b_{2}\\\vdots &\\a_{m1}x_{1}&+a_{m2}x_{2}+\cdots +a_{mn}x_{n}=b_{m},\end{aligned}}} $$Usando notação vetorial também pode-se escrever:
$$x_{1}{\begin{bmatrix}a_{11}\\a_{21}\\\vdots \\a_{m1}\end{bmatrix}}+x_{2}{\begin{bmatrix}a_{12}\\a_{22}\\\vdots \\a_{m2}\end{bmatrix}}+\cdots +x_{n}{\begin{bmatrix}a_{1n}\\a_{2n}\\\vdots \\a_{mn}\end{bmatrix}}={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{m}\end{bmatrix}}$$ou em notação matricial:
$$A{\mathbf {x}} = {\mathbf {b}}$$$$A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{m1}&a_{m2}&\cdots &a_{mn}\end{bmatrix}},\quad {\mathbf {x}}={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\quad {\mathbf {b}}={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{m}\end{bmatrix}}$$A notação matricial é util e nos mostra que a solução pode ser obtida de:
$${\mathbf {x}} = A^{-1}{\mathbf {b}}$$ou seja, através da obtenção da matriz inversa de $A$. No entanto, computacionalmente, inversão de matrizes é uma operação em geral custosa e ineficiente.
Para uma discussão dos tipos de sistemas assim como alguns métodos de solução veja:
https://en.wikipedia.org/wiki/System_of_linear_equations
Pode-se mostrar que em geral os sistemas lineares tem 3 tipos de comportamente:
Abaixo vemos um exemplo de sistema 3x3 e sua interpretação geometrica:
Na figura acima vemos a interpretação geométrica das possibilidades para o caso de sistemas 2x2.
os sistemas podem ser quadrados (nxn) ou não (mxn)
De um modo geral, os sistemas podem ser resolvidos usando dois tipos de técnicas:
regra de Cramer
Eliminação de Gauss
LU decomposition
Gauss-Jacobi
Gauss-Seidel
Consiste basicamente na obtenção da matriz escalonada de A a partir de operações elementares.
A matriz escalonada tem a forma:
$$\left[{\begin{array}{rrrr}2&-3&2&1\\0&1&-4&8\\0&0&0&35\end{array}}\right]$$As operações básicas que podem ser aplicadas a matriz, sem alterar a solução são:
Para uma aplicação passo a passo do método veja:
https://en.wikipedia.org/wiki/Gaussian_elimination#Example_of_the_algorithm
O método da LU Decomposition é uma variante do método de eliminação de Gauss. Nele obtem-se através de algorítmos uma decomposição do sistema na forma:
$${\displaystyle {\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}}={\begin{bmatrix}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\end{bmatrix}}{\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{bmatrix}}.} $$Após esta decomposição a solução é facilmente obtida através da substituição. Neste caso ao contrário do metodo de eliminação de Gauss, precisamos fazer duas sequencias de substituições. Este método é em geral mais eficiente do que a eliminação de Gauss para sistemas grandes.
O Numpy e o Scipy tem pacotes específicos para a solução de sistemas lineares. O pacote Scipy também contém um pacote de funções para realizar operações com algebra linear.
AUma das funções mais usadas é a solve (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.linalg.solve.html) que se encontra no pacote de funções de algebra linear do Numpy. O pacote de algebra linear do numpy contém inumeras funções para realização de operações com vetores e matrizes. Para uma lista e exemplos veja:
https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html
Um tutorial com exemplos de aplicação da decomposição LU podem ser encontradas aqui:
https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/linalg.html
In [5]:
import numpy as np
A = np.array([[2,1,4,1],
[3,4,-1,-1],
[1,-4,1,5],
[2,-2,1,3]],float)
v = np.array([-4,3,9,7],float)
N = len(v)
# Eliminação de Gauss
for m in range(N):
#divide pelo elemento da diagonal
div = A[m,m]
A[m,:] = A[m,:]/div
v[m] = v[m]/div
# subtrair das linhas seguintes
for i in range(m+1,N):
mult = A[i,m]
A[i,:] = A[i,:] - mult*A[m,:]
v[i] = v[i] - mult*v[m]
print '---------------------------'
print A
print v
print ''
# recorrencia para obter solução
x = np.zeros(N,float)
for m in range(N-1,-1,-1):
x[m] = v[m]
for i in range(m+1,N):
x[m] = x[m] - A[m,i]*x[i]
print "a solução é: ", x
# usando os modulos do numpy
print np.linalg.solve(A,v)
print ''
print "a solução com modulo solve é: ", x
In [8]:
#Solução de sistemas lineares usando Scipy e decomposição LU
import scipy.linalg as linalg
import numpy as np
#define A
A = np.array([[2., 1., 1.], [1., 3., 2.], [1., 0., 0.]])
#define B
B = np.array([4., 5., 6.])
#call the lu_factor function
LU = linalg.lu_factor(A)
#solve given LU and B
x = linalg.lu_solve(LU, B)
print "Solutions:\n",x
print ' '
# vamos ver o que foi feito com a matriz original
P, L, U = linalg.lu(A)
print P
print ' '
print L
print ' '
print U
In [7]:
print "%10.20f"% A[2,0]
A outra classe de métodos é a dos iterativos. Neste tipo de solução a meta é obter aproximações sucesivas da solução do sistema.
Para os sistemas lineares que estamos estudando o método pode ser definido pela seguinte recorrencia para Gauss-Jacobi:
$$ x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n. $$Essencialmente o que estamos fazendo aqui é isolar cada variável em função das outras do sistema. Tendo reorganizado o sistema dessa forma adotamos um chute inicial para a solução. Esse chute é substituido no sistema e a partir dele obtemos uma nova solução. Este processo é repetido até que as soluções sucessivas tenham uma diferença menor que um valr pré-estabelecido.
No caso de Gauss-Seidel, fazemos uma pequena melhoria usando as estimativas já obtidas no processo de recursão:
$$ {\displaystyle x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.} $$
In [ ]: