Fatoração LU

Teoria Completa

Introdução

Nesse tópico vamos começar a ver outro método direto para a resolução de sistemas lineares, a fatoração LU. A ideia aqui é fatorar a matriz A em um produto de matrizes L U , onde L é uma matriz triangular inferior e U é uma matriz triangular superior.

Primeiro vamos ver o método e depois como podemos aplicar ok? As matrizes L e U podem ser encontradas por meio de fórmulas, mas isso pode dificultar nossa vida mais pra frente. Sendo assim eu tenho uma péssima notícia pra quem achou que estava livre da Eliminação de Gauss, porque nós vamos usá-la aqui haha. Vamos ver um exemplo genérico de sistema 3 × 3 :

a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3

A matriz A do sistema A x = b seria dada por:

A = a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33

Agora se a gente aplicar a Eliminação de Gauss como faríamos? Calcularíamos os multiplicadores né? Por exemplo:

m 21 = a 21 a 11

Lembrou?

Depois de aplicar toda a eliminação de Gauss íamos chegar à matriz triangular superior que sempre obtemos no final da eliminação:

U = u 11 u 12 u 13 0 u 22 u 23 0 0 u 33

E vamos a chamar de U . Já aqueles multiplicadores vamos jogar numa matriz triangular inferior assim:

L = 1 0 0 m 21 1 0 m 31 m 32 1

Beleza até aqui? Depois de fazermos isso, teremos duas matrizes triangulares tal que A = L U ! Se não tiver acreditando no exemplo eu te mostro que dá mesmo certo haha. Beleza, já vimos como fazer, mas como isso pode ser útil pra gente?

Aplicações

Vou citar aqui algumas aplicações em que essa fatoração é útil.

Primeira aplicação: Cálculo de determinante

Como A = L U , det ⁡ A = det ⁡ L U = det ⁡ L ∙ det ⁡ U

E o determinante de uma matriz triangular é simplesmente o produto dos elementos da sua diagonal! Como a diagonal de L são sempre vários 1 , teremos det ⁡ L = 1 . Já o det ⁡ U será det ⁡ U = u 11 ∙ u 22 ∙ … ∙ u n n . Então ficamos com:

det ⁡ A = det ⁡ L ∙ det ⁡ U = 1 ∙ u 11 ∙ u 22 ∙ … ∙ u n n = u 11 ∙ u 22 ∙ … ∙ u n n

Então para achar o determinante de A basta multiplicarmos os elementos da diagonal de U .

Segunda aplicação: Inversa de A

Como A = L U , A - 1 = L U - 1 = U - 1 L - 1

E como U e L são triangulares, fica mais fácil inverter elas do que inverter a A .

Mas nenhuma dessas aplicações é o nosso objetivo primário aqui. Nós estaremos mais focados na terceira:

Terceira aplicação: Resolução de sistemas lineares de mesma matriz A

Você provavelmente reparou que, diferente da Eliminação de Gauss, aqui nós não usamos o vetor b pra nada! E isso foi proposital. O poder desse método está em resolver vários problemas A x = b para um mesmo A , mas vários b diferentes! Uma vez que A = L U , temos:

A x = b ⇔ L U x = b ⇔ L U x = b

Se chamarmos agora U x de y . Teremos:

U x = y   s i s t e m a   t r i a n g u l a r   s u p e r i o r

L y = b   s i s t e m a   t r i a n g u l a r   i n f e r i o r

Usando a Eliminação de Gauss, sempre que trocássemos o vetor b teríamos que repetir todo aquele processo chato. Usando a fatoração L U , não importa qual o vetor b , nós sempre teremos que resolver apenas 2 sistemas triangulares, que são muito mais fáceis!

Vamos ver um exemplo! Resolva o sistema abaixo para os vetores b dados:

x 1 + x 2 = b 1 2 x 1 + 3 x 2 = b 2

b a = 5,6 ; b b = 0,1 ; b c = 2,3

Temos o melhor caso possível para o nosso método, uma matriz A inalterada e o vetor b variando. Então vamos aplica-lo:

A = 1 1 2 3

L 2 ← L 2 - m 21 L 1 = L 2 - 2 1 L 1 :

1 1 0 1

Logo,

L = 1 0 2 1 ; U = 1 1 0 1

E se fizermos o produto L U veremos que:

L U = 1 0 2 1 1 1 0 1 = 1 1 2 3 = A

Agora podemos resolver os sistemas pedidos:

b a = 5,6 :

L U x = b a

L y = b a :

1 0 2 1 y 1 y 2 = 5 6

y 1 = 5

10 + y 2 = 6 ⇒ y 2 = - 4

U x = y :

1 1 0 1 x 1 x 2 = 5 - 4

x 2 = - 4

x 1 - 4 = 5 ⇒ x 1 = 9

Logo,

x 1 , x 2 = 9 , - 4

b b = 0,1 :

L U x = b b

L y = b b :

1 0 2 1 y 1 y 2 = 0 1

y 1 = 0

y 2 = 1

U x = y :

1 1 0 1 x 1 x 2 = 0 1

x 2 = 1

x 1 + 1 = 0 ⇒ x 1 = - 1

Logo,

x 1 , x 2 = - 1,1

Por fim,

b c = 2,3 :

L U x = b c

L y = b c :

1 0 2 1 y 1 y 2 = 2 3

y 1 = 2

4 + y 2 = 3 ⇒ y 2 = - 1

U x = y :

1 1 0 1 x 1 x 2 = 2 - 1

x 2 = - 1

x 1 - 1 = 2 ⇒ x 1 = 3

Logo,

x 1 , x 2 = 3 , - 1

Tranquilo né? Vamos ver o algoritmo desse negócio agora!

Algoritmo

Bom, podemos dividir esse algoritmo em três fases bem distintas. Fase 1 será a eliminação de Gauss. A diferença é que dessa vez você só vai receber a matriz A e não a A e a b . Além disso, essa fase deve te dar duas matrizes ao invés de uma só! A Fase 2 receberia a matriz L vinda da Fase 1 e o vetor b para resolver o sistema L y = b , que seria resolvido não por retrosubstituição, mas por substituição simples. A Fase 3 receberia o resultado da fase 2 , o vetor y , e a matriz U da Fase 1, e aplicaríamos nosso algoritmo da retrosubstituição. Vamos lá então:

Algoritmo da Fase 1:

Passo 1. Dados iniciais (o que seu programa vai receber de você):

a. Matriz quadrada A n × n

Passo 2. Definimos uma matriz com apenas 1 na diagonal e atribuímos à variável L .

Passo 3. Definimos U = A

Passo 4. Percorremos as linhas de U com uma variável i (de 1 a n ).

Passo 5. Escolher pivô da i -ésima linha de U   e armazenar na variável p . O pivô será sempre o i -ésimo número da i -ésima linha.

Passo 6. Para cada ‘ i ’, percorremos as linhas de U de i + 1 até n com uma variável k .

Passo 7. Definimos o multiplicador m i k = a k i p .

Passo 8. Guardar m i k na posição linha i e coluna k de L .

Passo 9. Para cada ‘ k ’, percorremos as colunas de U i até n com uma variável j .

Passo 10. a k j = a k j - m i k ∙ a i j

Pronto, tá feito o nosso algoritmo. Daí vamos tirar nossas matrizes U e L . Agora precisamos criar um algoritmo para a substituição simples.

Algoritmo da Fase 2:

Passo 1. Dados iniciais (o que seu programa vai receber de você):

a. Matriz quadrada triangular L n × n (que você vai receber do algoritmo acima).

b. Vetor solução b 1 × n

Passo 2. y 1 = b 1 a 11

Passo 3. i = n + 1

Passo 4. y i = b i

Passo 5. s o m a = 0

Passo 6. j = 1

Passo 7. s o m a = s o m a + a i j ∙ x j

Passo 8. j = j + 1

Passo 9. Enquanto j < i , repetimos os Passos 7 e 8. Quando j ≥ i prosseguimos.

Passo 10. y i = y i - s o m a a i i

Passo 11. i = i + 1

Passo 12. Enquanto i ≠ n + 1 , repetimos os Passos 4 ao 11. Quando i = n + 1 , temos a nossa resposta que pode ser dada em forma de uma lista y 1 , y 2 , … , y n .

Para a Fase 3 teremos exatamente o mesmo algoritmo da retrosubstituição. Só iremos alterar as entradas.

Algoritmo da Fase 3:

Passo 1. Dados iniciais (o que seu programa vai receber de você):

a. Matriz quadrada triangular U n × n (que você vai receber do algoritmo da Fase 1).

b. Vetor solução y 1 × n (recebido do algoritmo da Fase 2).

Passo 2. x n = b n a n n (lembrando que n é a dimensão da sua matriz)

Passo 3. i = n - 1

Passo 4. x i = b i

Passo 5. s o m a = 0

Passo 6. j = n

Passo 7. s o m a = s o m a + a i j ∙ x j

Passo 8. j = j - 1

Passo 9. Enquanto j > i , repetimos os Passos 7 e 8. Quando j ≤ i prosseguimos.

Passo 10. x i = x i - s o m a a i i