Departamento de Ciência da Computação - IME-USP

Quarto Exercício-Programa
Data de entrega: até 27 de novembro de 1998

Prof. Ronaldo Fumio Hashimoto

1  Introdução

O objetivo deste EP é implementar a resolução de um sistema de equações lineares aplicando o Método de conhecido como Método de Decomposição LU.

Definição: Um sistema de equações lineares (de n equações e n incógnitas) é um sistema do tipo:

onde os a_i,j e os y_i (i,j = 0,1,2,¼, n-1) são dados, e os x_i (i = 0,1,2,¼,n-1) são as incógnitas do sistema. Por exemplo,

é um sistema de 3 equações e 3 incógnitas.

Sejam

Então, podemos escrever o sistema (1) como A·x = y.

Um sistema dado pode não ter solução, ou ter uma ou mais soluções (inclusive infinitas). No exemplo (*****), x₀ = -2, x₁ = 1 e x₂ = 2 é a única solução do sistema. As vantagens de uma representação por matrizes são, além de ser uma notação mais compacta, permite trabalhar em um computador.

2  Casos Especiais: Matrizes Triangulares

Temos dois casos especiais que são de resolução imediata:

1. A matriz A é uma matriz triangular superior, isto é, a_i,j = 0 se i > j. Assim, os valores abaixo da diagonal de A são todos nulos. Por exemplo, seja o sistema:

A resolução é então trivial: da última equação achamos x₂ = 1. Com x₂ e a segunda equação, achamos x₁ = 0. E com x₁ e x₂ achamos, da primeira equação, o valor de x₀ = -5. Este método chama-se de ``back-substitution''.

2. A matriz A é uma matriz triangular inferior, isto é, a_i,j = 0 se i < j. Assim, os valores acima da diagonal de A são todos nulos. Por exemplo, considere dado o seguinte sistema:

Novamente, a resolução fica trivial: de primeira equação, x₀ = 2. Com x₀ e a segunda equação, obtemos x₁ = -1, e finalmente, da última equação, e com os valores de x₀ e x₁, temos que x₂ = 4. Este método chama-se de ``forward-substitution''.

3  Método de Decomposição LU

Dada uma matriz A, vamos tentar achar (se for possível) uma matriz L triangular inferior e uma matriz U triangular superior tais que A = LU. Nesse caso,

como A = LU, temos:

e podemos fazer:

Fazendo Ux = z, temos:

Dessa forma, resolver o sistema Ax = y equivale a resolver dois sistemas lineares (triangulares) simples: primeiramente o sistema (triangular inferior) Lz = y por ``forward-substitution'' e depois o sistema (triangular superior) Ux = z por ``back-substitution''.

3.1  Cálculo da matriz triangular inferior L

Dada uma matriz A = (a_i,j), queremos obter a partir dela uma matriz L com zeros acima da diagonal. Para isso, vamos começar ``zerar'' gradualmente os elementos a_0,1, a_0,2, ¼ , a_0,n-1. Seja l tal que a_0,1-la_0,0 = 0. Então l = a_0,1/a_0,0 (assumir a_0,0 não nulo). Calculamos, então, os novos valores para a segunda coluna da seguinte forma: a primeira coluna menos l vezes a segunda coluna. Ou seja,

Isto é,

Claramente, obtemos a¢_0,1 = 0 (pela escolha do l). Repetimos o procedimento para a terceira coluna, isto é, para ``zerar'' a_0,2, recalculamos a terceira coluna a partir da primeira. Seja l¢ = a_0,1/a_0,0. Então

Ou seja,

Claramente, obtemos a¢_0,2 = 0 (pela escolha do l¢). Repetindo esse processo, vamos ``zerar'' todos elementos a_0,1, a_0,2, ¼ , a_0,n-1, obtendo:

Devemos prosseguir para ``zerar'' a segunda linha. Para ``zerar'' a¢_1,2, recalculamos a terceira coluna a partir da segunda (utilizando a coluna já modificada no passo anterior!). Seja b = a¢_1,2/a¢_1,1 (assumir a¢_1,1 não nulo). Então:

Ou seja,

Assim, obteremos a¢¢_1,2 = 0.

Da mesma forma, para obtermos a¢¢_1,3 = 0, recalculamos a quarta coluna a partir da segunda. Seja b¢ = a¢_1,3/a¢_1,1. Então:

Ou seja,

Seguindo este procedimento, vamos ``zerar'' todos os elementos a¢_1,2, a¢_1,3, ¼ , a¢_1,n-1, obtendo:

Prosseguimos assim para a terceira linha, zerando a¢¢_2,3, a¢¢_2,4, ¼ , a¢¢_2,n-1, sempre usando a terceira coluna. Depois, a quarta linha usando a quarta coluna. Seguindo esse processo para as demais linhas, até obter uma matriz L diagonal inferior.

3.2  Cálculo da matriz triangular superior U

A matriz U é obtida simplesmente armazenando os coeficientes usados em cada passo do procedimento (anterior) para obter a matriz L, da seguinte forma:

Definimos u_i,i = 1 e u_i,j = 0, para i=1, 2, ¼, n-1 e j < i. Assim obteremos a matriz:

Pode-se verificar que A = LU, e que a decomposição é única.

Observação: Se em algum momento a¢^¼¢_i,i = 0, o método deve parar.

4  Um exemplo de cálculo de L e U

Primeiro passo, zerar a_0,1 e a_0,2, isto é, 4 e 2. Para isso, primeiro calculamos

l = a_0,1/a_0,0 = 4/2 = 2, e recalculamos a segunda coluna:

a¢_0,1 = 0,

a¢_1,1 = a_1,1-la_1,0 = 7-2·4 = -1,

a¢_2,1 = a_2,1-la_2,0 = -7-2·(-2) = -3.

Para zerar a_0,2, calculamos l¢ = a_0,2/a_0,0 = 2/2 = 1, e assim, recalculamos a terceira coluna (a partir da primeira):

a¢_0,2 = 0,

a¢_1,2 = a_1,2-la_1,0 = 7-1·4 = 3,

a¢_2,2 = a_2,2-la_2,0 = 5-1·(-2) = 7.

Os multiplicadores são l = 2 e l¢ = 1, assim u_0,1 = 2 e u_0,2 = 1 (como U tem 1 na diagonal, u_0,0 = 1). A matriz depois deste passo será:

Segundo passo, zerar o resto da segunda linha, isto é, a¢_1,2 = 3. Para isso, calculamos l = a¢_1,2/a¢_1,1 = 3/(-1) = -3 (agora usamos a segunda coluna!). Assim, recalculamos a terceira coluna:

a¢¢_1,2 = 0,

a¢¢_2,2 = a¢_2,2-la¢_2,1 = 7-(-3)·(-3) = -2.

O multiplicador foi l = -3, e assim u_1,2 = -3. Uma vez que U é triangular superior e tem em sua diagonal, então devemos ter u_1,0 = 0 e u_1,1 = 1.

Agora acabou o processo, u_2,0 = u_2,1 = 0 e u_2,2 = 1, e obtemos:

Dessa forma, basta resolvermos:

E finalmente, resolvendo:

Observe que, como os computadores tem uma precisão limitada, necessariamente vamos ter erros de arredondamento nos cálculos, e então uma solução para o sistema Ax = y não vai ser exata. Seja x a solução encontrada pelo seu programa. O erro pode ser medido achando-se a diferença entre o vetor Ax e y. A diferença é calculada da seguinte forma:

5  O problema a ser resolvido!

Fazer um programa em Linguagem C tal que, lida uma matriz real quadrada A de ordem n e um real vetor y de n elementos, resolva (se for possível) o sistema linear Ax = y, seguindo os seguintes passos:

1. Achar uma decomposição A = LU;

2. Resolver (por forward-substitution) Lz = y;

3. Resolver (por back-substitution) Ux = z;

4. Achar o erro Erro = ||Ax - y|| da solução encontrada x, onde

||v|| = max{|v_i| : 0 £ i £ n-1}.

5. Imprimir A, L e U, z, x e o erro E.

6  O que você DEVE fazer e utilizar

Para isso, você deve definir as seguintes funções:

void copy_mat (int n, double A[][MAX_MAT], double B[][MAX_MAT]);

Copia a matriz quadrada A de ordem n na matriz B quadrada, também de ordem n.

int calc_lu (double A[][MAX_MAT], double LU[][MAX_MAT]);

Esta função obtêm L e U. Inicialmente devemos ter L = A. Não use nenhuma matriz auxiliar nesta função! Na matriz LU deve ficar a matriz L na parte inferior de LU e a matriz U (sem a diagonal de 1) na parte superior de LU. Para o exemplo (*****), teríamos a matriz LU como:

LU = é ê ê ê ê ë 2 2 1 4 -1 -3 -2 -3 -2 ù ú ú ú ú û

Observação: Se na execução de calc_lu achamos algum l_i,i = 0, então a função deve retornar um valor especial para permitir ao programa principal detectar tal situação, imprimindo na tela uma mensagem de erro.

void prod (int n, double A[][MAX_MAT], double x[], double y[]);

Calcula o produto da matriz A_n×n pelo vetor x_n×1, guardando o resultado em y_n×1, ou seja,

yi = n-1 å j = 0  ai,j·xj

void forw_subst (int n, double LU[][MAX_MAT], double y[], double z[]);

Resolve Lz = y por forward-substitution, onde L é a parte inferior de LU.

void back_subst (int n, double LU[][MAX_MAT], double z[], double x[]);

Resolve Ux = z por back-substitution, onde U é a parte superior (sem considerar a diagonal) de LU, sendo que cada elemento de sua diagonal de U deve valer 1.

double norma (int n, double x[]);

Devolve o valor de ||x ||. Você pode usar a função fabs (double a) da biblioteca de C para calcular o valor absoluto de a.

void imprime_vet (int n, double v[]);

Imprime o vetor x com n elementos.

void imprime_mat (int n, double A[][MAX_MAT]);

Imprime a matriz quadrada A de ordem n.

7  Faça um teste!

Imprima a saída do seu programa para a seguinte entrada:

A: uma matriz 10×10 (obs.: .7030 representa 0.7030);

y: um vetor de 10 elementos.