aResolução de Sistemas Lineares1. O que é um sistema linear1.1. IntroduçãoUma equação linear de uma variável (incógnita), x, e um dado, y, é simplesmente(1)y=𝛼x+𝛽,cuja solução geral é dada porx=(y-𝛽)/𝛼,onde 𝛼 e 𝛽 são constantes do problema. Como estamos resolvendo para x, tanto 𝛼, 𝛽 e y são "dados" do problema, assim normalmente não escrevemos y explicitamente, mas apenas a combinação y-𝛽, assim temos o mesmo sistema como: 𝛼x=b, onde b=(y-𝛽), cuja solução é x=b/𝛼, como mostrado. No entanto podemos ter o caso de que ambos x e y são variáveis a serem resolvidas e os dados são apenas as constantes:a1x+a2y=b,que igualmente representa uma linha reta, significando que não existe euma solução única. Por exemplo(2)x+2y=3,é a reta abaixo em azul (y=-x/2+3/2). Se acrescentarmos uma segunda equação, por exemplo:(3)2x-y=-1,teremos a reta em verde (y=2x+1). As duas equações juntas (Eq. 2 e Eq. 3) formam um sistema:{a

x+2y=-3

2x-y=-1

|,cuja única solução será o ponto de interseção:P=(-1,-1).

Exercício: Cheque esta resposta resolvendo explicitamente o sistema. 1.1.1. Python codeVamos testar a representação gráfica usando Python. import numpy as npimport matplotlib.pyplot as plt x=np.linspace(-4,4,50)y1 = -3/2-x/2y2 = 2*x-1 plt.figure(1)plt.title('Duas Funções')plt.plot(x,y1)plt.plot(x,y2)plt.grid(True)ax=plt.gca()ax.set_xlim([-4.5, 4.5])ax.set_ylim([-3.5, 1.5])ax.set_aspect('equal') 1.2. Multiplas váriáveisEm geral, quando temos mais n variáveis, x1,x2,...,xn, e igualmente n dados b1,b2,...,bn, podemos ter um sistema linear, dado por:(4){a

a11x1+a12x2+…a1nxn	=	b1
a21x1+a22x2+…a2nxn	=	b2
…
an1x1+an2x2+…annxn	=	bn

|,onde x1, x2,..., xn são as n incógnitas, a11, a12, etc..., são os n2 coeficientes das equações, e b1, ..., bn são os n termos constantes. Note: Quando temos duas variáveis (x1 e x2, n=2), cada linha do sistema representa uma linha reta (uma dimensão). Se tivermos três variáveis (n=3), cada linha representará um plano (duas dimensõeS), e a solução será a intersecão dos planos, dentro de um espaço de dimensão 3. Em geral, para n variáveis, cada linha representa um hiperplano de n-1 dimensões, e a solução será a interseção desses hiperplanos. 1.3. Notação MatricialPodemos escrever o mesmo sistema Eq. 2 em termos de matrizes.Definindo a matrix A (quadrada, n×x), os vetores x e b (n×1) como sendo: A=[a

a11	a12	…	a1n
a21	a22	…	a2n
an1	an2	…	ann

] x=[a

x1

x2

…

xn

] b=[a

b1

b2

…

bn

] o sistema pode ser reescrito como: Ax=b. NOTE: Por simplicidade vamos escrever, daqui para frente, xi, aij e bi para indicar qualquer uma dessas quantidades, e usaremos os números apenas quanto nos referirmos a uma componente particular. Exercício: Dada a matriz de coeficientes A e o vetor constante B abaixo, escreva o sistema linear equivalente. Note que em 3 dimensões voce pode usar x, y, e z ao invés de x1, x2, x3.A=[a

2	1	1
3	2	1
2	1	2

]. B=[a

1

2

3

] Escreva as matrizes A e B para o sistema:{a

x+2y=-3

2x-y=-1

|. 1.4. Solução:Dado um sistema linear na forma matricial,AX=B, o vetor de soluções X é dado pela inversão da matriz A: X=A-1B.Isso só será possível caso a matriz A tenha inversa, e uma consição para isso é que seu determinante não seja zero, det(A)≠0. 1.4.1. Programa Python import numpy as np A=np.array([[2,1,1],[3,2,1],[2,1,2]])B=np.array([[1,2,3]]).Tadet = np.linalg.det(A)if(adet!=0): print("Det (A) = ", adet, " OK.")else: print("Det(A)=0 - Não Existe Inversa.") X=np.linalg.inv(A)@Bprint(X) 1.5. Matrizes TriangularesMatrizes triangulares superiores tem a forma U=(a

3	4	1
0	2	3
0	0	4

).Existem também as diagonais inferiores:L=(a

3	0	0
4	2	0
1	3	4

).Note que se U é triangular superior, UT é triangular inferior (como no exemplo). Tudo que será discutido aqui pode ser transformado para se aplicar em matrizes triangulares superiores. 2. Sistemas triangularesConsidere o sistema abaixo, chamado de Triangular Superior: a

2x1-x2+x3	=	2
x2+2x3	=	3
x3	=	1

,onde arrumamos as variáveis x1, x2, x3, nas mesmas colunas. Note que a primeira linha tem todas as 3 variáveis, a segunda linha tem apenas 2, e a última linha tem apenas 1 variável. Esse sistema é facilmente resolvido, substituindo de baixo para cima: x3=1x2=3-2x3→x2=3-2=12x1-x2+x3=2→x1=1. Para termos um algorítmo para solucionar um sistema triangular superior, vamos primeiramente pensar no caso geral com n=3:{a

a11x1	+a12x2	+a13x3	=b1
	a22x2	+a23x3	=b2
		a33x3	=b3

|,e assim, a solução geral é (indo de baixo para cima): {a

x3	=	b3/a33
x2	=	(b2-a23x3)⏠⏣⏣⏣⏡⏣⏣⏣⏢soma/a22
x1	=	(b1-a13x3-a12x2⏠⏣⏣⏣⏣⏡⏣⏣⏣⏣⏢soma)/a11

|, onde, depois de calcular a linha de cima, podemos calcular o valor do termo "soma" para então calcular o valor de x subsequente em baixo. Exercício:a) Usando a solução acima, resolva o sistema:{a

x1+4x2+x3	=	1
6x2+4x3	=	1
x3	=	1

|.b) Reescreva esse sistema como uma equação matricial. Defina, A, B, e X. 2.1. Caso geral No caso geral, de n variáveis e n equações, teremos{a

a11x1	+a12x2	+…	+a1(n-1)xn-1	+a1nxn	=	b1
	a22x2	+…	+a2(n-1)xn-1	+a2nxn	=	b2
		…
			a(n-1)(n-1)xn-1	+an-1,nxn	=	bn-1
				annxn	=	bn

|,onde seguimos o mesmo procedimento, onde a variável xi será igual à (bi-"resto") dividido por aii:{a

xn	=	bn/ann
xn-1	=	(bn-1-an-1,nxn)⏠⏣⏣⏣⏣⏡⏣⏣⏣⏣⏢soma/an-1,n-1
xn-2	=	(bn-2-an-2,n-1xn-1-an-1,nxn)⏠⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏢soma/an-2,n-2
		=(bn-2-n∑i=n-1an-2,ixi)1 an-2,n-2
	…
x2	=	(b2-a2nxn-…-a23x3)/a22
		= (b2-n∑i=3a2ixi)1 a22
x1	=	(b1-a1nxn-…-a13x3-a12x2)/a11
		=(b1-n∑i=2a1ixi)1 a11

|.Confira com o caso n=3:{a

x3	=	b3/a33
x2	=	(b2-a23x3)⏠⏣⏣⏣⏡⏣⏣⏣⏢soma/a22
x1	=	(b1-a13x3-a12x2⏠⏣⏣⏣⏣⏡⏣⏣⏣⏣⏢soma)/a11

|. Assim a fórmula geral, para cada variávels xj, onde j=1...n, é dada por(5)xj=1

ajj(bj-n∑i=j+1ajixi),onde precisamos começar com xn e ir descendo até x1. Note que, para haver solução, a diagonal da matrix não pode ser zero, ou seja aii≠0. Exercício: Confira a Eq. 5 no caso de n=3. 2.2. Criando listas no PythonNo Python, e na biblioteca numpy, existe algumas maneiras de se criar sequências de números. 1. range(a,b,s)(a) Cria uma sequência de números inteiros começando no valor "a" e indo até, mas sem incluir, "b", indo de passos do tamanho de "s". Apenas números inteiros.(a) range(1,10)=range(1,10,1)=[1,2,3,4,5,6,7,8,9](b) range(1,10,2)=[1,3,5,7,9](c) range(2,10,2)=[2,4,6,8](d) range(10,1,-1)=[10,9,8,7,6,5,4,3,2]2. numpy.arange(a,b,s)(a) A mesma coisa, mas cria um vetor do numpy (ndarray), e usa o tipo de número que estiver nos argumentos (inteiro ou float).(a) numpy.arange(1,3,.5)=[1.0, 1.5, 2.0, 2.5]3. numpy.linspace(a,b,n)(a) Cria uma sequencia de "a" até "b", com "n" valores uniformemente espaçados. Sempre inclui o ponto final. Usa o tipo numérico que estiver nos argumentos.(a) numpy.linspace(0,2,5)=[0.0 0.5 1.0 1.5 2.0] Ex. compare numpy.arange(0,0.31,.1) com numpy.linspace(0,0.3,4) 2.2.1. "Concertando" a função rangeSe você estiver confuso com relação a função range, voce pode (e deve) usa-la da seguinte forma:Se voce quiser uma lista de "a" até "b" (incluindo "b"), de passo "s", voce deve fazer: range(a,b±1,s), por exemploaini=1afim=10apas=2for ii in range(aini, afim+1, apas): print(ii) for ii in range(afim, aini-1, -apas): print(ii) Mas quando acessamos ararys, temos que lembra que seus índices começam em zero! Por exemplo, o vetor abaixov=(1,2,3,4,5),v1=1,v2=2,...v5=5,seria escrito matematicamente como vi=i,mas em python (ou em C também), teríamos que fazerv[0]=1v[1]=2...v[4]=5 ou v=numpy.zeros(5)for i in range(0,5): v[i]=i+1print(v) ou for i in range(1,5+1): v[i-1]=iprint(v) Vai do gosto. Voce vai ter que sempre pensar no que voce está fazendo! 2.3. Algorítmo (Sistema Triangular Superior) e programa Python Note que o cerne do algorítmo é a fórmula:(6)xj=1

ajj(bj-n∑i=j+1ajixi),que será dividido em duas partes, a soma (7)s=n∑i=j+1ajixi,e o cálculo de xj (8)xj=1

ajj(bj-s),para cada valor do índice j, das linhas do sistema, de n até 1. Adicionando os cheques de dimensionalidade e a condição de que nenhum elemento diagonal seja zero. O algorítmo segue como: Declara função com dois argumentos de entrada (A,B) e um vetor X de saída:1. Obtem a forma de A (alin, acol). (a) Checa se é quadrada (alin=acol), se não, emite erro e sai(b) Define a dimensionalidade do problems: N=alin.2. Checa se A é triangular superior: (verifica se o triângulo inferior é zero) (a) Para ilin de 1 até N:i. Para jcol de 1 até ilin-1:A. Se A(ilin, jcol) != 0, emite erro e sai.3. Obtem a forma de B (blin, bcol).(a) Se blin≠N ou bcol≠1, emite erro e sai.4. Declara o vetor coluna X com forma (Nx1), com zeros.(a) Iteração principal5. Para jlin de [n, n-1, ..., 1]: (loop nas linhas do sistema triangular)(a) soma=0(b) Para icol=jlin+1 até n: (Eq. 4)i. soma=soma+A(jlin, icol)*x(icol,1)(c) X(jlin)=(B(jlin)-soma)/A(jlin,jlin) (Eq. 5)6. Retorna X import numpy as np def trisup(A,B): # 1. Obtem a forma de A alin,acol=A.shape # Testa se é quadrada if(alin!=acol): print('Matriz A não é quadrada') return None n=alin # 2. Testa se A é triangular superior for ilin in range(0,n): for jcol in range(0,ilin-1): if(A[ilin][jcol]!=0): print('Matriz A não é triangular inferior') print(A) return None # 3. Obtem a forma de B blin,bcol=B.shape if(bcol!=1 or blin!=n): print('Vetor B não é compativel') #4. Declara X X=np.zeros([n,1]) #5. Iteração principal # Calcula linha por linha, da última até a primeira for jlin in range(n,1-1,-1): s=0 for icol in range(jlin+1,n+1): s=s+A[jlin-1][icol-1]*X[icol-1,0] X[jlin-1]=(B[jlin-1,0]-s)/A[jlin-1][jlin-1] return X # Teste da rotina trisup# Define uma matriz triangular superiorA=np.array([[1,1,2],[0,2,1],[0,0,5]])# Define um vetor coluna - note que usamos o método transpose().B=np.array([[1, -3, -5]]).transpose()B=np.array([[1],[2],[3]]) print('Matriz A=')print(A)print('Vetor B=')print(B) X=trisup(A,B)print('Solução X=')print(X) Exercícios:Faça testes para chegar as salvaguardas do programa - ou seja, para checar os testes internos. Resolva o sistema {a

x1	+2x2	+3x3	+4x4	=	30
	5x2	+6x3	+7x4	=	56
		8x3	+9x4	=	60
			10x4	=	40

| Calcule o determinante usando a rotina definida no curso, e calcule manualmente (observe que matrizes triangulares tem uma fórmula super simples para o determinante!). 3. Sua bibliotecaAntes de terminarmos, vamos cria a sua própria biblioteca de rotinas. Coloque a rotína de soma, multiplicação, traço, determinante e agora de solução de sistemas lineares em um único arquivo. Chamado de linlib.py: import numpy as np def msoma(A,B): .... def mtransposto(A): .... def mtraco(A): .... def mmult(A,B): .... def det(A): .... E crie um arquivo com todos os testes, chamado teste.py import numpy as npimport linlib as ll # coloque os testes ....# por exemplo:A=np.array([[1,2,3],[4,5,6],[7,8,9]])print('det(A)=',ll.det(A)) .... 4. Eliminação GaussianaPara resolvermos sistemas lineares em geral (que não sejam triangulares), precisamos notar o seguinte:1. Toda equação permanece válida quando realizamos a mesma operação matemática em ambos os lados:(a) Ex: (3x+2y=5 )×4 ↔ 12x+8y=202. Num sistema, se uma equação for somada (ou subtraída) de outra:(a) Gerará uma terceira equação válida(b) Este irá substituir uma das duas originais {a

3x+2y	=5 1
2x-2y	=10 2

|→{a

5x	=15	1+2
2x-2y	=10	2

|→{a

x	=3	1+2 /5
x-y	=5	2 /2

|→{a

x	=3
y	=-2

| 3. A permuta (troca) de linhas não altera o resultado: {a

3x+2y	=5
2x-2y	=10

|⇄{a

2x-2y	=10
3x+2y	=5

|implica que ambas as equações abaixo geram o mesmo resultado:(a

3	2
2	-2

)(a

x

y

)=(a

5

10

)⇆(a

2	-2
3	2

)(a

x

y

)=(a

10

5

). A Eliminação Gaussiana usa essas duas propriedades para transformar um sistema geral em um sistema triangular. O método é o seguinte.Considere o sistema: {a

a11x1+a12x2+…a1nxn	=	b1	1
a21x1+a22x2+…a2nxn	=	b2	2
a31x1+a32x2+…a3nxn	=	b3	3
…
an1x1+an2x2+…annxn	=	bn	n

|Observe a sequencia de operações:1. Começe com a linha 1:1.1) Divida cada equação pelo seu respectivo coeficiente de x1:{a

x1+a12/a11x2+…a1n/a11xn	=	b1/a11	1 /a11= 1'
x1+a22/a21x2+…a2n/a21xn	=	b2/a21	2 /a21= 2'
x1+a32/a31x2+…a3n/a31xn	=	b3/a31	3 /a31= 3'
…
x1+an2/an1x2+…ann/an1xn	=	bn/an1	n /an1= n'

| Note 1: se o coeficiente a11=0, este não poderá ser usado e esta linha deve ser trocada por outra (pivoteamento). Note 2: se em alguma linha, ai1=0, ela é ignorada neste passo e no seguinte. 1.2) Subtraia a primeira equação de todas as outras:{a

x1+a12/a11x2+…a1n/a11xn	=	b1/a11
(a22/a21-a12/a11)x2+…(a2n/a21-a1n/a11)xn	=	b2/a21-b1/a11
(a32/a31-a12/a11)x2+…(a3n/a31-a1n/a11)xn	=	b3/a31-b1/a11
…
(an2/an1-a12/a11)x2+…(ann/an1=a12/a11)xn	=	bn/an1-b1/a11

| produzindo o sistema{a

x1+a12'x2+…a1n'xn	=	b1'	1'
a22'x2+…a2n'xn	=	b2'	2'
a32'x2+…a3n'xn	=	b3'	3'
…
an2'x2+…ann'xn	=	bn'	n'

| 2. Repita o procedimento agora com a linha 2 como cabeça. O algorítmos para esse procedimento pode ser construído assim: • Para cada linha i, de 1 até n, faça: – Se aii=0 troque a ordem das linhas por uma que não tenha esse problema (algorítmo não mostrado)– Para cada linha j, de i até n (de i para baixo), faça:* Se aji=0, passe para a linha seguinte.* bj'=bj/aji* Para cada coeficiente k, (da linha j), de j até n, faça:· ajk'=ajk/aji.– Para cada linha j, de i+1 até n (abaixo de i), faça:* Se aji=0, passe para a linha seguinte* bj''=bj'-bi'* Para cada coeficiente k, (da linha j), de j até n, faça:· ajk''=ajk'-aik'. A troca de linhas é um algorítmo a parte e será feito diferentemente. Este algorítmo pode ser escrito em python assim:(lembre que: 1. As matrizes tem que ter os índice subtraído de 1aij → a[i-1][j-1] 2. Para gerar uma lista de "a" até "b" usando a função range, fazemos: range(a,b+1) def elimgauss(A,B,n): a=copy(A) # necessário para não estragar o valor de A for i in range(1,n+1): # Pivoteamento (troca de linhas) # A,B=pivoteia(A,B,i,n) for j in range(i,n+1): if(a[j-1][i-1]==0): continue b[j-1]=b[j-1]/a[j-1][i-1] for k in range(j,n+1): a[j-1][k-1]=a[j-1][k-1]/a[j-1][i-1] for j in range(i+1,n+1): if(a[j-1][i-1]==0): continue b[j-1]=b[j-1]-b[i-1] for k in range(j,n+1): a[j-1][k-1]=a[j-1][k-1]-a[i-1][k-1] # Zera diagonal inferior (normalmente não é necessário) for i in range(1,n+1): for j in range(1,i+1): a[i-1][j-1]=0 return a,b # testeimport numpy as npA=np.array([[1,2,3],[4,5,6],[7,8,9]]);B=np.array([[1],[2],[3]]);print(A)print(B) At,Bt=elimgauss(A,B,3) print(At)print(Bt) 4.1. Eliminação gaussiana com pivoteamentoPara evitar problemas com coeficientes zero ou muito próximos de zero, uma estratégia comum é a do pivoteamento: na i-ésima iteração, ao invés de usar a linha i automaticamente, busca-se a linha j com coeficiente aji de maior valor absoluto, e usa-se essa linha j no lugar da linha i. Algorítmo:• ilin=i• am=0• Para cada linha j, de i até n, faça:– Se(abs(aji)>am)* am=aji* ilin=j• Se ilin≠i, troque as linhas ilin e i:– Para cada coluna k, de i até n:* Atemp[1][k]=A[i][k]* A[i][k]=A[ilin][k]* A[ilin][k]=Atemp[1][k]– Btemp=B[i]– B[i]=B[ilin]– B[ilin]=Btemp– Ou em python def pivoteia(A1,B1,i1,n): # A1 - matriz # B1 - vetor coluna # i - linha em questão (de 1 até n) # n - dimensionalidade de A e B import copy A=copy.copy(A1) B=copy.copy(B1) # O programa abaixo está configurado para usar índices de 0 até n-1. É necessário então subtrair i1 de 1. i = i1-1 ilin=i am=0 Atemp=np.zeros(n) for j in range(i,n): if(abs(A[j][i])>am): am=A[j][i] ilin=j if(ilin!=i): for k in range(i,n): Atemp[k]=A[i][k] A[i][k]=A[ilin][k] A[ilin][k]=Atemp[k] Btemp=B[i][0] B[i]=B[ilin][0] B[ilin][0]=Btemp return A,BA=np.array([[1,2,3],[4,5,6],[7,8,9]]);B=np.array([[1],[2],[3]]);print(A)print(B) At,Bt=pivoteia(A,B,0,3) print(At)print(Bt) 4.2. Rotina Final para resolver sistemas linearesVamos juntar todas as três rotinas para resolver sistemas lineares em geral: def linsis(A,B): At,Bt=elimgauss(A,B) X=trisup(At, Bt) return X Resolva o sistema: {a

2x+2y-3t = 2z+1

3x-2y+t-z=0

-4y+x-3z -2 =0

x+t = -1

| 5. A Matriz InversaFinalmente vamos calcular a matriz inversa. Ela é definida pela equação matricial: (9)AA-1=A-1A=I,onde I é a matriz identidade. Todas as matrizes são da forma n×n. Uma condição necessária para a existência da inversa é que detA≠0.Para resolver a Eq. 9 coluna por coluna de A-1. Vamos chamar a coluna i da matriz inversa como Xi, e assimAXi=Bi,onde Bi é o vetor com todas a componentes zero exceto na componente i, por exemplo, se n=3 e i=2, B2=(a

0

1

0

).Então, para resolver para A-1, basta resolver as n equações para X1, X2, ... até Xn. O algorítmo seria:1. Checa se det(A)≠02. Faz loop de 1 até n para calcular todas as colunas de A-1 2.1. Define Bi 2.2. Calcula Xi usando eliminação gaussiana e solução triangular superior 2.3. Coloca o valor de Xi na coluna i da nova matriz A-1. Ou em python def inversa(A): if(det(A)==0): print('det(A)=0') return None for ii in range(1,n+1): Bi=np.zeros(n) Bi[ii-1]=1 Xi=linsis(A,Bi) for jj in range(1,n+1): Ainv[jj-1][ii-1]= Xi[jj-1] return Ainv # TesteA=np.array([[1., 2],[3, 4]])Ai=inversa(A)print(mmult(A,Ai))print(mmult(Ai,A)) Aicorreto=np.array([[4,-2],[-3,1]])*(-.5)print(Ai-Aicorreto) 6. Comandos do Numpy para realizar essas tarefasimport numpy A=np.array([[1., 2],[3, 4]])B=np.array([[2,3],[4,5]]) V1=np.array([[1,2,3]])V2=np.array([[0,1,2]]) # Soma de matrizesprint(A+B)print(A-B) # MultiplicaçãoA2=np.matmul(A,A)print(A2) print(A@B) # InversaAi=np.linalg.inv(A)print(Ai) print(np.matmul(A,Ai))print(A@Ai)print(Ai@A) # Produto escalar (Vetor com vetor)print(V1.dot(V2)) # E assim existem vários outros comandos.... 7. Exemplos de Problemas Lineares7.1. Circuito de ResistoresUm circuiro elétrico de várias malhas com apenas baterias e resistores pode ser escrito como um sistema linear de equações, cuja incógnita são as correntes elétricas:

E1E2E2R1R1R1R1R2i1i1i3i2i2i2i1Na malha esquerda:a3-2i1-4i3-6-2i1=0-3-i1(2+2)-4i3=0-3-4i1-4i3=04i1+4i3=-3.Na malha direita - vamos rodar no sentido anti-horário:a-R2i3-E2-R1i2+E2-R1i2=0-4i3-6-2i2+6-2i2=0-4i3-4i2=0Lei dos nós:i1+i2=i3{a

4i1		+4i3	=-3
	-4i2	-4i3	=0
i1	+i2	-i3	=0

|R1=2ΩR1=4Ωℇ1=3Vℇ1=6V Definindo então as matrized A e B, teremos: import numpyA=np.array([[4,0,4],[0,-4,-4],[1,1,-1]])B=np.array([[-3,0,0]]).TI=np.linalg.inv(A)@Bprint(I)