Métodos numéricos

2. Sistemas de equações lineares

  1. Sistemas de equações lineares
    1. Representação matricial
    2. Método de eliminação de Gauss
    3. Método de Gauss-Jordan

2.1. Representação matricial

Um sistema linear de ordem n é um sistema com n variáveis x 1 , x 2 , …, x n , que satisfazem n equações lineares:

(2.1)
a 11 x 1 + a 12 x 2 +· · · + a 1 n x n = b 1 a 21 x 1 + a 21 x 2 +· · · + a 2 n x n = b 2 .. . a n 1 x 1 + a n 2 x 2 +· · · + a n n x n = b n

Se as equações são linearmente independentes, existe sempre uma única solução. Claro está que basta dar a lista de equações e variáveis aos comandos solve ou linsolve do Maxima para obter a solução, mas o objetivo deste capítulo não é obter a resposta mas sim mostrar, passo a passo, como funcionam alguns métodos de resolução.

Chama-se matriz aumentada do sistema à matriz:

(2.2)
a 11 a 12 ·· · a 1 n b 1 a 21 a 22 ·· · a 2 n b 2 .. . . . . . . . . . . .. . a n 1 a n 2 ·· · a n n b n

2.2. Método de eliminação de Gauss

A qualquer equação no sistema pode somar-se outra das equações, multiplicada por uma constante k , e a solução do sistema não se altera. Na matriz aumentada, essa operação corresponde a substituir uma linha L i por L i + k L j . Essa propriedade pode ser usada para obter uma matriz triangular superior, em que a i j é zero, se i < j :

(2.3)
a 11 a 12 ·· · a 1 n b 1 0 a 22 ·· · a 2 n b 2 .. . . . . . . . . . . .. . 0 0 ·· · a n n b n

que corresponde a outro sistema de equações com a mesma solução do sistema original (as constantes a i j e b j não são as mesmas a i j e b j da matriz inicial).

A última equação nesse sistema permite obter facilmente o valor da variável x n

(2.4)
a n n x n = b n x n = b n a n n

Este resultado pode ser substituído na penúltima equação para obter o valor de x n 1 ,

(2.5)
a n 1 , n 1 x n 1 + a n 1 , n x n = b n 1 x n = b n 1 a n 1 , n 1 a n 1 , n b n a n 1 , n 1 a n n

e assim sucessivamente, até se obter o valor de todas as variáveis.

Exemplo 2.1

Resolva o sistema de equações:

8 x + 4 . 5 z 12 . 4 = 3 . 2 y 2 x 0 . 8 y = 7 . 6 + 6 . 1 z x + 3 y + 10 . 2 z = 8 . 4

usando o método de eliminação de Gauss.

Resolução. No Maxima, convém guardar as 3 equações em 3 variáveis

(%i1) eq1: 8*x + 4.5*z - 12.4 = 3.2*y$
(%i2) eq2: 2*x - 0.8*y = 7.6 + 6.1*z$
(%i3) eq3: x + 3*y + 10.2*z = 8.4$

Este sistema pode ser escrito da mesma forma do sistema geral 2.1 para obter a matriz aumentada. A função augcoefmatrix realiza esse procedimento dando como resultado a matriz aumentada. É necessário dar dois argumentos a augcoefmatrix, a lista das equações e a lista das variáveis, na ordem que se quiser. Neste caso, a matriz aumentada pode ser guardada na variável m com o seguinte comando

(%i4) m: float (augcoefmatrix ([eq1, eq2, eq3], [x, y, z]));
(%o4)    8 . 0 3 . 24 . 5 12 . 42 . 0 0 . 8 6 . 1 7 . 61 . 03 . 010 . 2 8 . 4

optou-se por usar a função float para evitar que a matriz seja armazenada em números racionais.

Há que ter em conta que os números na última coluna não são b i mas sim b i (ou seja, no sistema 2.1 todas as equações foram igualadas a zero).

Antes de começar com o processo de eliminação pode fixar-se um número reduzido de algarismos significativos, por exemplo 4, para os resultados que serão apresentados a seguir, de modo a evitar que as matrizes tenham muitos algarismos:

(%i5) fpprintprec: 4$

Como está a ser usado o formato de vírgula flutuante de dupla precisão (8 bytes), internamente os cálculos continuarão a ser feitos com 16 algarismos significativos, mas todos os resultados apresentados no ecrã serão arredondados para 4 algarismos significativos.

Para reduzir a matriz à forma triangular, começa-se por substituir a segunda linha, L 2 , por L 2 - (1/4)  L 1 , onde L 1 representa toda a primeira linha da matriz. O objetivo é fazer com que a 21 fique igual a zero (note-se que 1/4 foi calculado dividindo a 21 por a 11 ).

No Maxima, a função rowop( m , i , j , k ) produz uma nova matriz, em que a linha L i da matriz m é substituída por L i - k L j e as restantes linhas ficam iguais. Assim sendo, a substituição L 2 - (1/4)  L 1 L 2 pode ser feita com o seguinte comando:

(%i6) m: rowop (m, 2, 1, 1/4);
(%o6)    8 . 0 3 . 24 . 5 12 . 40 . 00 . 0 7 . 225 4 . 51 . 03 . 010 . 2 8 . 4

A seguir, para eliminar a 31 = 1.0 usando novamente a 11 = 8.0, faz-se a substituição L 3 - k L 1 L 3 , onde k = a 11 / a 31 = 1/8. Ou seja,

(%i7) m: rowop (m, 3, 1, 1/8);
(%o7)    8 . 0 3 . 24 . 5 12 . 40 . 00 . 0 7 . 225 4 . 50 . 03 . 49 . 637 6 . 85

Neste momento deveria usar-se a 22 para eliminar a 23 . No entanto, isso não é possível porque a 22 é nula ( k = a 23 / a 22 não existe). Em casos como este, troca-se a linha em questão (neste caso a segunda) com uma linha posterior (neste caso apenas pode ser com a terceira). Essa troca de linhas equivale a alterar a ordem das equações, sem que por isso se altere a solução. Usa-se a função rowswap para trocar a ordem de duas linhas:

(%i8) m: rowswap (m, 2, 3);
(%o8)    8 . 0 3 . 24 . 5 12 . 40 . 03 . 49 . 637 6 . 850 . 00 . 0 7 . 225 4 . 5

ficando assim uma matriz diagonal superior, que permite passar à segunda fase do método.

Para resolver as equações associadas a cada linha da matriz é conveniente definir uma lista com as variáveis, na mesma ordem em que foram usadas para encontrar a matriz aumentada, seguidas por 1:

(%i9) vars: [x, y, z, 1]$

Multiplicando a terceira linha da matriz m pela lista de variáveis, obtém-se a última equação, na qual só aparece a variável z ; resolve-se essa equação para encontrar o valor de z , que será guardado para ser usado nos passos seguintes

(%i10) sol: float (solve (m[3].vars));
(%o10)    [ z = -0.6228 ]

Multiplicando a segunda linha da matriz pela lista de variáveis e substituindo o valor que já foi encontrado para z , obtém-se uma equação que depende unicamente de y . A resolução dessa equação conduz ao valor de y e é conveniente acrescentar esse resultado na mesma lista em que foi armazenado o valor de z

(%i11) sol: append (float (solve ( subst (sol, m[2].vars))), sol);
(%o11)    [ y = 3.78, z = -0.6228 ]

Finalmente, repete-se o comando anterior, mas agora multiplicando pela primeira linha da matriz, para determinar o valor de x

(%i12) sol: append (float (solve ( subst (sol, m[1].vars))), sol);
(%o12)    [ x = 3.412, y = 3.78, z = -0.6228 ]

2.3. Método de Gauss-Jordan

O método de eliminação de Gauss tem a vantagem de não alterar o valor do determinante da matriz com coeficientes a i j e, como tal, o determinante pode ser calculado facilmente multiplicando os valores a i i na diagonal da matriz triangular obtida.

O método de Gauss-Jordan consiste em transformar a matriz aumentada até os coeficientes a i j corresponder à matriz identidade, ou seja, iguais a 1 se i = j ou zero caso contrário. A grande vantagem é que os valores finais de b i são os valores das variáveis, sem ser preciso realizar mais contas. Em contrapartida, o determinante da matriz já não é o produto dos elementos na diagonal da matriz reduzida.

Para conseguir transformar a matriz na matriz identidade, começa-se por dividir a primeira equação por a 11 , para que a 11 fique igual a 1, e a seguir usa-se essa primeira linha para eliminar o primeiro elemento em todas as outras linhas. A seguir divide-se a segunda linha por a 22 , para que este fique igual a 1 e usa-se essa segunda linha para eliminar o segundo elemento em todas as outras linhas, incluindo a primeira.

Não existe um comando específico do Maxima para dividir uma linha de uma matriz por uma constante; no entanto, a divisão por k pode ser feita com a função rowop( m , i , i , 1 1/ k ), que substitui a linha L i da matriz m , por L i - (1-1/ k )  L i = L i / k .

Exemplo 2.2

Determine a solução do seguinte sistema, pelo método de Gauss-Jordan

2 x 1 1 . 2 x 2 4 x 3 = 5 3 x 1 + 3 x 2 + 5 x 3 = 122 . 5 x 1 + 4 x 2 3 x 3 = 6

Resolução. Para mostrar como é possível obter soluções exatas, na resolução deste exemplo só será usada a função float no fim e os dois números de vírgula flutuante 1.2 e 2.5 na primeira e terceira equações serão escritos como as fracções 12/10 e 25/10.

Em vez de se usar a função augcoefmatrix, a matriz aumentada pode também ser inserida diretamente usando a função matrix

(%i13) m: matrix ([2,-12/10,-4,5],[-3,3,5,12],[25/10,4,-3,-6]);
(%o13)    2 6 5 45 33 5 12 5 24 3 6

Começa-se por dividir a primeira linha pelo seu primeiro elemento (ou seja 2):

(%i14) m: rowop (m, 1, 1, 1 - 1/m[1][1]);
(%o14)    1 3 5 25 2 33 5 12 5 24 3 6

A seguir, subtrai-se a primeira linha, multiplicada pelo primeiro elemento da segunda linha, para eliminar o primeiro elemento de segunda linha:

(%i15) m: rowop (m, 2, 1, m[2][1]);
(%o15)    1 3 5 25 20 6 5 139 25 24 3 6

e faz-se o mesmo com a primeira e a terceira linhas, para eliminar o primeiro elemento da terceira linha:

(%i16) m: rowop (m, 3, 1, m[3][1]);
(%o16)    1 3 5 25 20 6 5 139 20 11 22 49 4

Repete-se o mesmo procedimento para a segunda coluna; ou seja, divide-se a segunda linha pelo seu segundo elemento,

(%i17) m: rowop (m, 2, 2, 1 - 1/m[2][2]);
(%o17)    1 3 5 25 20 1 5 665 40 11 22 49 4

e usa-se essa nova linha para eliminar o segundo elemento na primeira e na terceira linhas:

(%i18) m: rowop (rowop (m, 1, 2, m[1][2]), 3, 2, m[3][2]);
(%o18)    10 5 249 40 1 5 665 40 0 79 12 813 8

Finalmente, repete-se o mesmo procedimento para a terceira coluna:

(%i19) m: rowop (m, 3, 3, 1 - 1/m[3][3]);
(%o19)    10 5 249 40 1 5 665 40 0 1 2439 158
(%i20) m: rowop (rowop (m, 1, 3, m[1][3]), 2, 3, m[2][3]);
(%o20)    10 0 2081 790 1 0 535 1580 0 1 2439 158

Os três números na última coluna são os valores exatos das 3 variáveis. Para aproximar esses números por números de vírgula flutuante e colocá-los numa lista, extrai-se a quarta linha da matriz transposta e aplica-se a função float

(%i21) float (transpose(m)[4]);
(%o21)    [ -26.34, 3.386, -15.44 ]

Note-se que se tivesse sido usada a função augcoefmatrix para construir a matriz, seria necessário mudar os sinais desses 3 valores obtidos.

Este algoritmo pode ser automatizado facilmente. O mesmo resultado obtido com os comandos (%i13) até (%i20) podia ter-se obtido com o seguinte ciclo:

(%i22) m: matrix ([2,-12/10,-4,5],[-3,3,5,12],[25/10,4,-3,-6]);
(%i23) for i:1 thru 3 do (
    m: rowop (m, i, i, 1 - 1/m[ i ][ i ] ),
    for j:1 thru 3 do (
        if (i # j) then m: rowop (m, j, i, m[ j ][ i ] )));
(%o23)    done

onde # é o operador "diferente de".

Se em alguma das iterações m[ i ][ i ] for nulo, 1/m[ i ][ i ] produz um erro. Para evitar esse tipo de erro seria necessário acrescentar instruções para conferir que m[ i ][ i ] não seja nulo e, caso contrário, trocar a linha L i com alguma das linhas L j com j > i .