Métodos iterativos estacionários para sistemas lineares

...

Transcript

1. LNCC UFRJ Sistemas lineares, métodos iterativos estacionários Prof. Paulo Bordoni
2. Os conceitos de ponto-fixo e de normas matriciais serão fundamentais para o estudo da convergência dos métodos iterativos para sistemas lineares. Vamos rever nosso estudo.
3. Seja : → uma função definida num espaço métrico M . Um ponto ∈ é um ponto-fixo de quando, e só quando, () = . Sim, vamos recordar o que vimos na aula passada: Começando pelo conceito de ponto-fixo:
4. As informações sobre ponto-fixo na scipy.optimize.
5. Também vimos que: 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) : ℝ → ℝ
6. Você se lembra Surfista, da sua contração: O ponto-fixo é a origem 0, 0 : 3 2 1 0 4 1 2 3 3 2 1 0 1 2 3 4 ↦ 3/4 1 0 0 1
7. Vimos também que: 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 ∈ .
8. E ainda a associação entre ponto-fixo e raiz de equação: “p é um ponto-fixo de uma função () se, e somente se, p é raiz de () = − ()” E a recíporoca: “r é raiz de uma função se, e somente se, r é ponto-fixo de () = − ()”.
9. Sim, lembro bem, Mestra: = ⟹ = − = 0 e também = 0 ⟹ 0 = − ⟹ = E que eu falei: “Deixe a fala do Sherlock por minha conta, Surfista!”
10. Vejam no programa a seguir. Obteremos o ponto-fixo de uma função = () usando o método oferecido pela SciPy. Automaticamente teremos a raiz da equação − = 0. Portanto, poderemos usar o método do ponto-fixo não apenas para achar raízes de equações.
11. Este programa mostra, numérica e visualmente, que: se () = − então = ⟺ = 0.
12. É mesmo, não vou esquecer jamais!
13. Na realidade, para resolver uma equação () = 0, usando o método do ponto-fixo, tudo que precisamos é escrever = () onde é alguma função construída manipulando algebricamente (). Complicou tudo, Filósofo.
14. Ficou com medo da liberdade, Surfista? Por exemplo, para a equação 3 + 42 − 10 = 0, poderemos buscar pontos fixos de: a. = 3 + 42 + − 10 b. = Τ 1 2 10 − 3 c. = − Τ ′ - o método de Newton-Raphson
15. Agora, passaremos a examinar métodos iterativos para resolver sistemas lineares.
16. 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, onde = − , procuraremos por pontos-fixo de alguma função linear = + obtida através de manipulação algébrica de ().
17. Para tudo funcionar, precisaremos garantir que a () seja uma contração. E como podemos garantir que essa : ℝ → ℝ é uma contração, Mestra?
18. Para responder tua pergunta, Loirinha apresentaremos antes resultados importantes sobre normas vetoriais e matriciais e limites de sequências em ℝ geradas por funções : ℝ → ℝ do tipo = + .
19. Prova-se que: • 2 = 1 = , onde 1 é o maior valor singular de • () ≤ , para qualquer norma natural .
20. Gerando matrizes com raio espectral = = 1.
21. Confiram que = 1 :
22. Confiram que lim →∞ = 0. Listei as 10 primeiras potências. Façam mais 40.
23. Uma matriz A é dita convergente quando lim →∞ = 0 As afirmações abaixo são equivalentes: 1. A é uma matriz convergente, 2. < 1, 3. lim →∞ = 0
24. No exemplo houve convergência apesar de = 1, mas é facílimo confirmar que se < 1 então A é uma matriz convergente. Sim Mestra! Prova-se que ∙ ≤ ∙ . Segue daí que ≤ e portanto que lim →∞ ≤ lim →∞ = 0 .
25. Para qualquer 0 ∈ ℝ seja a sequência de vetores de ℝ definida por +1 = + , ≥ 0. Então: 1. () é convergente se, e só se, < 1, 2. lim →∞ = , com satisfazendo = + .
26. 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.
27. 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.
28. Por exemplo, se = 4 3 −1 1 5 2 2 3 −7 , então: = 4 5 −7 = 1 2 3 = 3 −1 2
29. 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.
30. 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
31. 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.
32. Assuma que é ponto-fixo da . Então: = = + = = −−1 + + −1 Portanto = − + + e + + = isto é = .
33. 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.
34. 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 = .
35. 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 −
36. + (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
37. = 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
38. Clamo por exemplos concretos! Atenda ao Surfista, Mestre.
39. Para construir os métodos de Jacobi, Gauss-Seidel e Gauss-Seidel com relaxação, precisaremos de:
40. Detalhes da rotina diag( ):
41. Detalhes da rotina diagflat( ):
42. Detalhes das rotinas tril( ) e triu( ):
43. Apresentarei agora programas para os métodos de Jacobi, Gauss-Seidel e relaxação e exemplos de execução.
44. Eis aqui meu programa para o método de Jacobi. Construção da matriz do método de Jacobi.
45. A continuação:
46. Uma execução com convergência.
47. E uma outra execução sem convergência. Veja o raio espectral da matriz .
48. O código para Gauss-Seidel: Marquei onde ele difere do Jacobi.
49. O mesmo problema, e uma convergência bem mais rápida, ~65% mais.
50. 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 .
51. Agora o programa para o método de Gauss-Seidel com relaxação. Aos meus pés a parte de relaxação:
52. A continuação dele:
53. Uma execução do programa, com sobre- relaxação.
54. 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%.
55. Os resultados da sub-relaxação. O ganho com relação a Gauss- Seidel foi ~ 27%. Já com relação ao Jacobi foi ~75%.
56. 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
57. 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.
58. 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.
59. 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.
60. 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.
61. 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.
62. 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).
63. Agora o programa. Analise-o com atenção Surfista.
64. A sub-relaxação.
65. A sobre-relaxação.
66. A parte gráfica.
67. Faremos mais um exemplo para ver o que se repete. Mestres, os resultados do exemplo anterior se repetem sempre?
68. 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
69. O método de Jacobi convergiu com 36 iterações.
70. E o método de Gauss-Seidel convergiu com 11 iterações. Bem mais rápido.
71. Eis os resultados gráficos obtidos:
72. Vejam que o número mínimo de iterações da sobre-relaxação e Gauss-Seidel são iguais!
73. Com a sub-relaxação, o ganho relativamente ao Gauss- Seidel foi de apenas 1 iteração!
74. Agora vou apresentar uma outra forma de implementar os três métodos, trabalhando diretamente com as equações.
75. 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.
76. 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.
77. 1. Considere cada uma das equações do sistema ( = 1, … , ): ෍ =1 = 2. Em cada equação, separe os termos da diagonal: + ෍ =1,≠ = , 3. Passe os termos do somatório para o segundo membro e divida tudo por , obtendo: = ൘ − ෍ =1,≠ , 4. Imponha que os do lado esquerdo dessa igualdade sejam os elementos +1 vetor +1 e os do lado direito os elementos do vetor e pronto: = ൘ − ෍ =1,≠ .
78. 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.
79. 1. Considere cada uma das equações do sistema ( = 1, … , ): ෍ =1 = 2. Em cada equação, separe os termos da diagonal: ෍ ≥ + ෍ < = , 3. Passe os termos do 2º somatório para o segundo membro: ෍ ≥ = − ෍ < , 4. Imponha que os do lado esquerdo dessa igualdade sejam os elementos +1 vetor +1 e os do lado direito os elementos do vetor e pronto: ෍ ≥ +1 = − ෍ < . Continua
80. 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:
81. 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 .
82. 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)
83. 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)
84. Implemente computacionalmente estes outros algoritmos, Surfista. Não vejo porquê, Mestra! Basta usar os programas do Mestre.
85. É 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. Um assunto para o futuro.
86. Tchau, até a próxima aula!