CNum_ALinComp_06_-_Met._Iter._para__Sist.Lin_e_matrizesesparsas.pdf

Transcript

1. LNCC UFRJ Métodos iterativos para sistemas lineares e matrizes esparsas Prof. Paulo Bordoni
2. Além dos métodos diretos, métodos iterativos, constituem a outra vertente para a resolução de sistemas lineares. Apresento uma visão geral dos especialistas:
3. Eles são fundamentais! Vejam os métodos iterativos no livro do Golub:
4. E no livro do Trefethen/Bau:
5. Outros dois livros fundamentais sobre o assunto, além do Golub e do Trefethen:
6. O “Templates” é um resumo mais focado na parte computacional do tema .
7. Os pré-condicionadores e sua importância.
8. Agora, passaremos a examinar os métodos iterativos mais simples para resolver sistemas lineares.
9. Construiremos três métodos iterativos, de ponto-fixo, para resolver o sistema linear = , conhecidos na literatura como métodos estacionários. Como desejamos achar uma “raiz” da equação vetorial () = 0, em que = − , procuraremos por pontos-fixo de alguma função linear = + onde a matriz é obtida através de manipulação algébrica de ().
10. Para tudo funcionar, precisaremos garantir que a () seja uma contração. E como podemos garantir que essa : ℝ → ℝ é uma contração, Mestra?
11. Como antes, precisamos garantir que − () < − com 0 < < 1. Ora, se = + é linear, temos − = ( − ) e assim ( − ) ≤ ∙ − . Portanto se 0 < < 1 teremos que é uma contração.
12. Para qualquer 0 ∈ ℝ seja a sequência de vetores de ℝ definida por +1 = + , ≥ 0. Então: 1. Se < 1, teremos que () é convergente, i.é, ∃ ∈ ℝ tal que lim →∞ = . 2. O limite satisfaz = + .
13. Então são válidos os seguintes limitantes superiores para o erro: 1. − ≤ − 0 , 2. − ≤ 1− 1 − 0 . Assuma que a sequência de vetores definida por +1 = + , ≥ 0, com 0 ∈ ℝ, é convergente.
14. Para construir a matriz T que define as funções vetoriais : ℝ → ℝ para os métodos de Jacobi, Gauss-Seidel e Gauss-Seidel com relaxação, usaremos a seguinte partição da matriz A: = + + onde: 1. é a diagonal da A, 2. L é a parte triangular inferior de A, abaixo da diagonal, 3. U é a parte triangular superior de A, acima da diagonal.
15. Por exemplo, se = 4 3 −1 1 5 2 2 3 −7 , então: = 4 5 −7 = 1 2 3 = 3 −1 2
16. No método de Jacobi a matriz T é definida por = −−1( + ) e o vetor c por = −1. No método de Gauss-Seidel a matriz T é definida por = − −1 e o vetor c por = ( − )−1.
17. Para a matriz A do exemplo, as matrizes e são: = − 4 0 0 0 5 0 0 0 −7 −1 0 3 −1 1 0 2 2 3 0 = − 4 0 0 1 5 0 2 3 −7 −1 0 3 −1 0 0 2 0 0 0 = 4 3 −1 1 5 2 2 3 −7 = −−1( + ) = −( + )−1
18. Afirmamos que: • ponto-fixo de ⟺ solução de = • ponto-fixo de ⟺ solução de = Agora, seja : ℝ → ℝ a função definida por = + , onde = para o método de Jacobi e = para o método de Gauss- Seidel.
19. Assuma que é ponto-fixo da . Então: = = + = = −−1 + + −1 Portanto = − + + e + + = isto é = .
20. Para o método de Gauss-Seidel as contas são muito parecidas e também vale a equivalência. Como esses passos podem ser revertidos, vale a volta.
21. Agora, seja a sequência de vetores de ℝ definida por +1 = () para = 0,1,2, … com 0 um vetor qualquer de ℝ. Tanto no caso de Jacobi como no de Gauss-Seidel. Se a : ℝ → ℝ for uma contração essa sequência convergirá. Então, sendo = lim →∞ , teremos que = .
22. Sim Professor. Vou esclarecer as ideias na próxima transparência. O método de Gauss-Seidel com fator de relaxação é uma combinação linear dos métodos de Jacobi e de Gauss-Seidel definida por: = + 1 − , para ≥ 0.
23. + (1 − ) A ideia combinação linear. Quando: • = 0 temos o método de Jacobi, • = 1 temos o método de Gauss-Seidel. Para: • 0 < < 1 fala-se em sub-relaxação, • > 1 fala-se em sobre-relaxação. Jacobi Gauss-Seidel
24. = 4 0 0 1 5 0 2 3 −7 −1 0 3 −1 0 0 2 0 0 0 + (1 − ) 4 0 0 0 5 0 0 0 7 −1 0 3 −1 1 0 2 2 3 0 Confirmando na prática: Sejam a matriz do método de Jacobi e a matriz do método de Gauss-Seidel já obtidas. A matriz do método de relaxação é, simplesmente = + 1 − , com ≥ 0. = 4 3 −1 1 5 2 2 3 −7
25. Clamo por exemplos concretos! Atenda ao Surfista, Mestre.
26. Para construir os métodos de Jacobi, Gauss-Seidel e Gauss-Seidel com relaxação, precisaremos de:
27. Detalhes da rotina diag( ):
28. Detalhes da rotina diagflat( ):
29. Detalhes das rotinas tril( ) e triu( ):
30. Apresentarei agora programas para os métodos de Jacobi, Gauss-Seidel e relaxação e exemplos de execução.
31. Eis aqui um programa para o método de Jacobi. Construção da matriz do método de Jacobi.
32. A continuação:
33. Uma execução com convergência.
34. E uma outra execução sem convergência. Veja o raio espectral da matriz .
35. O código para Gauss-Seidel: Marquei onde ele difere do Jacobi.
36. O mesmo problema, e uma convergência bem mais rápida, ~65% mais.
37. O mesmo problema, no qual o método de Jacobi não convergiu. Também sem convergência para GS; veja o raio espectral da matriz .
38. Agora o programa para o método de Gauss-Seidel com relaxação. Aos meus pés a parte de relaxação:
39. A continuação dele:
40. Uma execução do programa, com sobre- relaxação.
41. O resultado da sobre-relaxação. Note que o ganho com relação ao de Gauss-Seidel foi 50%. Já com relação ao de Jacobi foi ~82%.
42. Os resultados da sub-relaxação. O ganho com relação a Gauss- Seidel foi ~ 27%. Já com relação ao Jacobi foi ~75%.
43. Mestre, e como você conseguiu esses valores de ? Ganhou na loteria? Não Loirinha! Usei o programa que eu e a Mestra passaremos a analisar
44. Na primeira parte, construo as matrizes de Jacobi e de Gauss- Seidel, como antes. Depois o Mestre constrói a matriz = + (1 − ) . Para ∈ [0,1] teremos sub- relaxação e para > 1 a sobre- relaxação.
45. Há vários comentários a efetuar sobre este gráfico. Para comparação dos resultados, ampliamos cinquenta vezes o raio espectral, assim a faixa pontilhada horizontal em magenta corresponde a = 1.
46. Em preto e vermelho temos o gráfico da função : 0, +∞ → ℝ, ω ↦ (). Para ∈ [0,1] temos a sub-relaxação e nela: • 0 , é o raio espectral de Jacobi (máximo), • 1 , é o raio espectral de Gauss-Seidel (mínimo), marcados com um quadradinho magenta.
47. Em azul e verde temos o gráfico da função : 0, +∞ → ℕ, ω ↦ (). Para ∈ [0,1] temos a sub-relaxação e: • 0 é o nº de iterações de Jacobi, • 1 é o nº de iterações de Gauss-Seidel. Também marcados com um quadradinho magenta. O ponto de mínimo na sub-relaxação é 0,63 = 13, nas bolinhas azuis.
48. Em vermelho temos o gráfico da função ω ↦ , ω > 1, que corresponde à sobre-relaxação. Note que 1.301 … = 0.352 … é mínimo absoluto; depois cresce ilimitadamente.
49. O gráfico de ω ↦ , ω > 1 está em verde. O mínimo absoluto é 1.261 … = 9 iterações. Depois cresce ilimitadamente e lim →2.748… = +∞, caindo a zero para > 2.748 … . O valor limite é atingido quando = 1 ( linha pontilhada magenta horizontal).
50. Agora o programa. Analise-o com atenção Surfista.
51. A sub-relaxação.
52. A sobre-relaxação.
53. A parte gráfica.
54. Faremos mais um exemplo para ver o que se repete. Mestres, os resultados do exemplo anterior se repetem sempre?
55. Trabalharemos com o sistema linear abaixo, cuja solução é = [ 1 1 0 − 1 1 ]. 7 −4 2 0 0 1 6 3 1 0 0 1 0 2 0 1 8 3 0 4 6 3 1 1 5 1 2 3 4 5 = 3 6 −1 −4 3
56. O método de Jacobi convergiu com 36 iterações.
57. E o método de Gauss-Seidel convergiu com 11 iterações. Bem mais rápido.
58. Eis os resultados gráficos obtidos:
59. Vejam que o número mínimo de iterações da sobre-relaxação e Gauss-Seidel são iguais!
60. Com a sub-relaxação, o ganho relativamente ao Gauss- Seidel foi de apenas 1 iteração!
61. Agora vou apresentar uma outra forma de implementar os três métodos, trabalhando diretamente com as equações.
62. Na forma de equações o sistema = fica: 11 1 + 21 1 + 12 2 + a22 2 + 13 3 + 23 3 + ⋯ + 1 = 1 + 2 = 2 ⋮ 1 1 + 2 2 + 2 2 + ⋯ + = Lembro que no método de: • Jacobi a matriz T é definida por = −−1( + ) e o vetor c por = −1. • Gauss-Seidel a matriz T é definida por = − −1 e o vetor c por = ( − )−1.
63. No caso de Jacobi, em termos das equações do sistema = , a igualdade +1 = −1 − + + se escreve, para cada linha = 1, … , : +1 = 1 ෍ =1, ≠ ( − ) . Esta fórmula pode ser obtida seguindo o algoritmo descrito na próxima transparência.
64. 1. Considere cada uma das equações do sistema: ෍ =1 = , = 1, … , 2. Em cada equação (para cada valor de ): a) separe os termos da diagonal: + ෍ =1,≠ = , b) Passe os termos do somatório para o segundo membro e divida tudo por : = / − ෍ =1,≠ ( / ) , c) Imponha que os do lado esquerdo dessa igualdade sejam os elementos do vetor +1: +1 = / − ෍ =1,≠ ( / )
65. No caso de Gauss-Seidel, a igualdade +1 = ( + )−1 − + se rescreve + +1 = − + que, em termos das equações do sistema = , é ෍ ≤ +1 = − ෍ > , = 1, … , . Esta fórmula pode ser obtida seguindo o algoritmo descrito na próxima transparência.
66. 1. Considere cada uma das equações do sistema: ෍ =1 = , = 1, … , 2. Em cada equação (para cada valor de ): a) Separe os termos até a diagonal ( ≤ ) dos depois ( > ): ෍ ≤ + ෍ > = , b) Passe os termos do 2º somatório para o segundo membro: ෍ ≤ = − ෍ > , c) Imponha que os do lado esquerdo dessa igualdade sejam os elementos do vetor +1: ෍ ≥ +1 = − ෍ < . Continua
67. 1. Da 1ª, calcule 1 +1 por 1 +1 = (1 − σ =2 )/11 , 2. Da 2ª, calcule 2 +1 por 2 +1 = [ 2 −( 21 1 +1 + σ =3 ) ]/22 , 3. Da 3ª, calcule 3 +1 por 3 +1 = [ 3 −( σ=1 2 3 +1 + σ =3 ) ]/33 , ⋮ ⋮ ⋮ n. Da nª, calcule +1 por +1 = [ − σ =1 −1 +1 ]/ . Agora, considerando equação por equação:
68. Desculpe-me, Mestra, mas me atrapalho muito com somatórios. Dá para fazer um exemplo numérico simples? Ok, vamos considerar o sistema linear = onde: = 4 3 −1 1 5 2 2 3 −7 e = 0 13 −10 .
69. Repetindo as equações do sistema para facilitar: ቐ 41 + 32 − 3 = 0 1 + 52 + 23 = 13 21 + 32 − 73 = −10 Seguindo o algoritmo, no caso de Jacobi, obtemos: 1 +1 = [0 −(32 − 3 )]/4 2 +1 = [13 −(1 + 23 )]/5 3 +1 = [−10 −(21 + 32 )]/(−7)
70. Tornando a repetir as equações do sistema para facilitar: ቐ 41 + 32 − 3 = 0 1 + 52 + 23 = 13 21 + 32 − 73 = −10 E, no caso de Gauss-Seidel, obtemos: 1 +1 = [0 −(32 − 3 )]/4 2 +1 = [13 −(1 +1 + 23 )]/5 3 +1 = [−10 −(21 +1 + 32 +1)]/(−7)
71. Repeti as equações do método de Gauss-Seidel, para enfatizar o aspecto que o diferencia do Jacobi: 1. Na 2ª equação o valor 1 +1 foi o recém obtido na 1ª equação; 2. Na 3ª equação os valores 1 +1 e 2 +1 foram os recém obtidos nas 1ª e 2ª equações. 1 +1 = [0 −(32 − 3 )]/4 2 +1 = [13 −(1 +1 + 23 )]/5 3 +1 = [−10 −(21 +1 + 32 +1)]/(−7)
72. Implemente computacionalmente estes outros algoritmos, Surfista. Mestra, não vejo por quê! Basta usar os programas do Mestre.
73. É importantíssimo. O motivo fundamental é que assim, apenas os coeficientes não-nulos da matriz A são utilizados. Isto nos remete ao estudo métodos numéricos para matrizes esparsas, nosso próximo tema.
74. O que são matizes esparsas, Mestre? São matrizes com muitos zeros. Tipicamente apenas 15% a 20% de seus elementos são não-nulos.
75. Na solução de problemas envolvendo equações diferenciais, tanto pelo método das diferenças finitas como dos elementos finitos. Como e onde elas aparecem?
76. Na SciPy temos a scipy.sparse (armazenamento) e a scipy.sparse.linalg (álgebra linear). A seguir a lista completa dos módulos.
77. Aqui estão todas as classes de armazena- mento do pacote. Nosso foco inicial será sobre as 4 que marquei.
78. Algum motivo especial por essa opção Mestre? O armazenamento COO é apropriado para o método dos elementos finitos e o DIA para o método das diferenças finitas.
79. O formato CSC (armazenamento por colunas) é muito parecido com o CSR:
80. As vantagens e desvantagens do formato de armazenamento por colunas também são semelhantes àquelas do formato por linhas:
81. Para salvar em disco e carregar uma matriz esparsa.
82. Para efetuar o produto . quando é esparsa.
83. No formato COO fornecemos a posição (, ) e o valor da matriz = [ ].
84. Mais detalhes sobre o formato COO:
85. Veja as possibilidades do formato CSR (armazena por linhas), Surfista:
86. Mais detalhes sobre o formato CSR:
88. Mestre, apresente um programa envolvendo a criação de uma matriz esparsa no formato COO. Sim! E como procedemos para passar do COO para outros formatos.
89. É o que farei. Nele detalharei a conversão de tanto para o formato tradicional como para o CSR, mais apropriado para operações algébricas. Sim Mestre, faça um programa envolvendo a criação de uma matriz já no formato COO.
90. Nesta parte deste programa exibo a entrada de dados para o formato COO.
91. 2 3 4 1. Obtendo o tamanho dim _ da matriz (usando conjuntos). 2. Construindo a matriz no formato COO → _coo. 3. Convertendo para o formato usual de matriz → 4. Convertendo para o formato CSR → _csr 5. Calculando o produto ∗ , os valores singulares de e seu raio espectral. 1 5 Acompanhe a descrição passo-a-passo do resto:
92. Uma execução do programa, aceitando a sugestão.
93. A continuação da execução. Observem no código que o produto ∗ é executado no formato CSR e depois transformado no formato comum para exibição, aos meus pés.
94. Detalhes importante para o formato COO, Surfista: 1. Você pode repetir entradas para a mesma posição 2. E não precisa ordenar os valores por posição. Já vi Mestra; os valores repetidos serão somados. Vi também que o formato CSR ordena os dados por linha e em cada linha.
95. Uma execução do programa, aceitando a sugestão.
96. Detalhes importante para o formato COO, Surfista: 1. Você pode repetir entradas para a mesma posição 2. E não precisa ordenar os valores por posição. Já vi Mestra; os valores repetidos serão somados. Vi também que o formato CSR ordena os dados por linha e em cada linha.
97. Nesse programa efetuamos duas conversões de formato. Eis a lista das conversões possíveis; ela está ao final da lista de métodos.
98. Tchau, até a próxima aula!