Métodos Iterativos para sistemas lineares e matrizes esparsas

Transcript

1. LNCC UFRJ Métodos iterativos para sistemas lineares e matrizes esparsas Prof. Paulo Bordoni
2. Para , , ∈ : 1. , ≥ 0 – não-negatividade 2. , = 0 ↔ = 3. , = , – simetria 4. , ≤ , + , − desigualdade triangular Um espaço métrico é uma entidade matemática constituída por um conjunto M e uma função : × → ℝ, chamada métrica do espaço, (para medir distâncias) que satisfaz as propriedades óbvias:
3. Todo espaço vetorial normado é um espaço métrico. Basta definir a métrica por , = − . O exemplo que mais utilizaremos é ℝ com , = − . Em seguida vem os ℝ com a distância definida através das normas que já vimos.
4. Seja : → uma função definida num espaço métrico M . Um ponto ∈ é um ponto-fixo de quando, e só quando, () = .
5. A função : ℝ → ℝ, é definida por = 3 possui três pontos-fixo. É verdade Mestra. Os pontos = 1, = 0 e = −1 são pontos-fixo. Todos satisfazem = .
6. Ah, Loirinha, o gráfico de () só corta a reta = nesses três pontos. Meu programa mostra isso. Mas, atenção para a escala! Fazendo as contas está claro, mas são só esses três?
7. O programinha do Surfista, que permite achar pontos-fixo por inspeção visual.
8. Uma função : → , um espaço métrico é uma contração quando, e apenas quando, existe uma constante ∈ [0,1) tal que , ≤ , , ∀, ∈ . f () y f () − − () ≤ − com ∈ [0,1) : ℝ → ℝ
9. Surfista, eis um exemplo de contração. O ponto-fixo é a origem 0, 0 : ↦ 3/4 1 0 1 0 3 2 1 0 4 1 2 3 3 2 1 0 1 2 3 4
10. Um dos resultados mais importantes de análise é o Teorema do ponto-fixo de Banach: Num espaço métrico completo toda contração : → admite um único ponto-fixo ∈ , (isto é = ).
11. Se : → é uma contração, podemos determinar seu ponto-fixo p iterativamente, construindo uma sequência 0 , 1 , … , , … tal que lim →∞ = . Sim, basta definir +1 = ( ), com 0 ∈ .
12. Uma das propriedade importantes é que (, ) ≤ 1 − (1 , 0 ) Que permite estimar a velocidade de convergência.
13. Para localizar as informações sobre ponto- fixo na scipy.optimize.
14. As informações sobre ponto-fixo na scipy.optimize.
15. O gráfico, mostrado lá atrás pelo Sherlock, foi obtido com este programa, que usa a função fixed_point( ) da scipy.
16. Este programa mostra, e desenha, algumas iterações ao ponto-fixo de uma função : [, ] → ℝ
17. Observem a convergência das iterações ao ponto-fixo = 0.5 de = ()/2. Claramente, 0.5 = 0.5.
18. Neste caso, também temos a convergência das iterações ao ponto-fixo de = 3−.
19. Observando o gráfico abaixo, podemos afirmar que − () ≤ 0.85 − , ∀ ∈ [0,1], i. é, que f é uma contração nessa grande vizinhança do ponto-fixo = 0.5478 … − () + ()
20. Por quê, Mestre? Simplesmente porque o gráfico de f fica inteiramente contido no cone de convergência. Acompanhe meu raciocínio:
21. Considere um ponto ∈ 0,1 , > . Observe na figura que − < < + (), i. é, − − + < < − + . Portanto − − < − < ( − ) ou, − () < ( − ), já que = (). Claro que o mesmo vale para ∈ 0,1 , < . E em p temos a igualdade, assim − () ≤ − , ∀ ∈ [0,1].
22. A convergência ao ponto fixo de = 0.6[1 + 2 cos()].
23. Por inspeção visual, podemos afirmar que = 0, = 1, = 2 e = 3 são pontos fixos de = + 1 2 (), cujo gráfico desenhamos abaixo.
24. Entretanto, mesmo escolhendo 0 = 1.95, muito mais próximo do ponto- fixo = 2, a convergência “vai para” o ponto-fixo = 1 Confira!
25. Mestre, tentei com 0 = 2.05, e também não convergiu para o ponto-fixo = 2. Mas convergiu para o ponto-fixo = 3.
26. O túnel é iluminado do início ao fim pelo teorema do ponto-fixo de Banach. A condição para convergência é que f seja uma contração local entorno de p, i. é: para algum ∈ [0,1), − () ≤ − , para todo nas proximidades do ponto-fixo p. Mestres, há alguma luz no fim do túnel?
27. Se: 1. ∈ [, ], 2. , ⊆ [, ], 3. f é derivável em (a, b) 4. ∃ ∈ 0, 1 tal que ′() ≤ , ∀ ∈ (, ) então, definindo uma sequência 0 , 1 , … , , … por +1 = ( ), com 0 ∈ [, ] poderemos afirmar que: a. lim →∞ = e p é o único ponto fixo de f em [, ] b. − ≤ 0 − , − 0 c. − ≤ Τ [ (1 − )] 1 − 0 , ∀ > 1. Às páginas 58,59 do Análise Numérica de Burden & Faires, 8ª ed. encontramos uma demonstração do Teorema:
28. 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:
29. Eles são fundamentais! Vejam os métodos iterativos no livro do Golub:
30. E no livro do Trefethen/Bau:
31. Outros dois livros fundamentais sobre o assunto, além do Golub e do Trefethen:
32. O “Templates” é um resumo mais focado na parte computacional do tema .
33. Os pré-condicionadores e sua importância.
34. Agora, passaremos a examinar os métodos iterativos mais simples para resolver sistemas lineares.
35. 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 ().
36. Para tudo funcionar, precisaremos garantir que a () seja uma contração. E como podemos garantir que essa : ℝ → ℝ é uma contração, Mestra?
37. Como antes, precisamos garantir que − () < − com 0 < < 1. Ora, se = + é linear, temos − = ( − ) e assim ( − ) ≤ ∙ − . Portanto se 0 < < 1 teremos que é uma contração.
38. 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 = + .
39. 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.
40. 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.
41. Por exemplo, se = 4 3 −1 1 5 2 2 3 −7 , então: = 4 5 −7 = 1 2 3 = 3 −1 2
42. 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.
43. 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
44. 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.
45. Assuma que é ponto-fixo da . Então: = = + = = −−1 + + −1 Portanto = − + + e + + = isto é = .
46. 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.
47. 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 = .
48. 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.
49. + (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
50. = 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
51. Clamo por exemplos concretos! Atenda ao Surfista, Mestre.
52. Para construir os métodos de Jacobi, Gauss-Seidel e Gauss-Seidel com relaxação, precisaremos de:
53. Detalhes da rotina diag( ):
54. Detalhes da rotina diagflat( ):
55. Detalhes das rotinas tril( ) e triu( ):
56. Apresentarei agora programas para os métodos de Jacobi, Gauss-Seidel e relaxação e exemplos de execução.
57. Eis aqui um programa para o método de Jacobi. Construção da matriz do método de Jacobi.
58. A continuação:
59. Uma execução com convergência.
60. E uma outra execução sem convergência. Veja o raio espectral da matriz .
61. O código para Gauss-Seidel: Marquei onde ele difere do Jacobi.
62. O mesmo problema, e uma convergência bem mais rápida, ~65% mais.
63. 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 .
64. Agora o programa para o método de Gauss-Seidel com relaxação. Aos meus pés a parte de relaxação:
65. A continuação dele:
66. Uma execução do programa, com sobre- relaxação.
67. 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%.
68. Os resultados da sub-relaxação. O ganho com relação a Gauss- Seidel foi ~ 27%. Já com relação ao Jacobi foi ~75%.
69. Mestre, e como você conseguiu esses valores de ? Ganhou na loteria? Não Loirinha! Existem alguns resultados teóricos para ajudar na escolha de , mas na prática trabalhamos tentativa e erro.
70. Agora vou apresentar a forma prática de implementar os três métodos, trabalhando diretamente com as equações, sem usar as funções diag( ), tril( ), triu( ) da NumPy.
71. Não entendo o motivo! Basta usar os programas mostrados pelo Mestre. A importância dos métodos iterativos reside em utilizar apenas os coeficientes não-nulos da matriz A, o que não acontece com as diag( ), tril( ), triu( ).
72. 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.
73. 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.
74. 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,≠ ( / )
75. 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.
76. 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
77. 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:
78. 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 .
79. 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)
80. 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)
81. 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)
82. Para implementar computacionalmente estes outros algoritmos precisaremos estudar métodos numéricos para matrizes esparsas, nosso próximo tema.
83. 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.
84. 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?
85. Na SciPy temos a scipy.sparse (armazenamento) e a scipy.sparse.linalg (álgebra linear). A seguir a lista completa dos módulos.
86. Aqui estão todas as classes de armazena- mento do pacote. Nosso foco inicial será sobre as 4 que marquei.
87. 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.
88. O formato CSC (armazenamento por colunas) é muito parecido com o CSR:
89. As vantagens e desvantagens do formato de armazenamento por colunas também são semelhantes àquelas do formato por linhas:
90. Para salvar em disco e carregar uma matriz esparsa.
91. Para efetuar o produto . quando é esparsa.
92. No formato COO fornecemos a posição (, ) e o valor da matriz = [ ].
93. Mais detalhes sobre o formato COO:
94. Veja as possibilidades do formato CSR (armazena por linhas), Surfista:
95. Mais detalhes sobre o formato CSR:
97. 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.
98. É 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.
99. Nesta parte deste programa exibo a entrada de dados para o formato COO.
100. 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:
101. Uma execução do programa, aceitando a sugestão.
102. 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.
103. 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.
104. Uma execução do programa, aceitando a sugestão.
105. 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.
106. 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.
107. Tchau, até a próxima aula!