| | | |

Matemática Numérica III

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

1.4 Método do gradiente conjugado

O método do gradiente conjugado é um dos métodos mais eficientes para a resolução de sistemas lineares

A𝒙=𝒃, (1.119)

com A matriz simétrica positiva definida n×n e bn. Foi desenvolvido por Magnus Hestenes2020endnote: 20Magnus R. Hestenes, 1906 - 1991, matemático americano estadunidense. Fonte: Wikipedia. e Eduard Stiefel2121endnote: 21Eduard Stiefel, 1909 - 1978, matemático suíço. Fonte: Wikipédia. em 1952.

A ideia é encontrar a solução do sistema linear como o ponto de mínimo da função quadrática

J(𝒙)=12(𝒙,A𝒙)(𝒙,𝒃), (1.120)

onde (,) denota o produto interno usual em n.

Note que J(𝒙) é estritamente convexa, pois A é simétrica positiva definida. Assim, J(𝒙) possui um único ponto de mínimo global, que é a solução do sistema linear A𝒙=b. De fato, supondo 𝒙 solução do sistema, temos

J(𝒙)=12(𝒙,A𝒙)𝒙Tb+12(𝒙,A𝒙)(𝒙,b) (1.121)
=12(𝒙𝒙,A(𝒙𝒙))12(𝒙,A𝒙) (1.122)

Logo, como A é simétrica positiva definida, temos que J(𝒙) é mínimo se, e somente se, 𝒙𝒙=𝟎, i.e. 𝒙=𝒙.

Ainda, observamos que o gradiente de J(𝒙) é dado por

J(𝒙)=A𝒙𝒃, (1.123)

i.e. J(𝒙) é igual ao vetor oposto do resíduo do sistema linear

𝒓(𝒙)=𝒃A𝒙. (1.124)

A solução 𝒙 é o único vetor tal que J(𝒙)=𝟎. Mais ainda, lembrando que o gradiente aponta na direção de maior crescimento da função, temos que o resíduo aponta na direção de maior decrescimento de J(𝒙). Isto nos motiva a considerar o método do gradiente (ou método da descida mais íngrime), que consiste nas iterações

𝒙(k+1)=𝒙(k)+αk𝒓(k), (1.125)

onde 𝒓(k)=𝒃A𝒙(k) é o resíduo na iteração k e αk>0 é um tamanho do passo escolhido, com 𝒙(0) uma aproximação inicial da solução do sistema. Em outras palavras, o método do gradiente constrói uma sequência de aproximações (𝒙(k))k da solução do sistema, movendo-se na direção do resíduo do sistema linear.

Podemos escolher o tamanho do passo αk de modo a minimizar J(𝒙(k+1)) ao longo da direção do resíduo 𝒓(k). Isto é, escolhemos αk como

αk=minα>0J(𝒙(k)+α𝒓(k)), (1.126)

o que é conhecido como busca linear exata. Fazendo as contas, obtemos

αk=(𝒓(k),𝒓(k))(𝒓(k),A𝒓(k)). (1.127)

Com tudo isso, temos a iteração do método do gradiente

𝒙(k+1)=𝒙(k)+αk𝒓(k), (1.128)
ondeαk=(𝒓(k),𝒓(k))(𝒓(k),A𝒓(k)). (1.129)
Observação 1.4.1.(Detalhe de implementação)

Observamos que o cálculo de αk requer o produto A𝒓(k) e que o resíduo 𝒓(k) requer o produto A𝒙(k). Essas multiplicações matriz-vetor são os passos computacionalmente mais custosos da iteração. Assim, para economizar uma multiplicação por iteração, podemos atualizar o resíduo como

𝒓(k+1)=𝒓(k)αkA𝒓(k), (1.130)

evitando o cálculo de A𝒙(k+1).

Código 6: mg.py
1import numpy as np
2def mg(A, b, x0, rtol=1e-5, atol=0.0, maxiter=100000):
3 x = x0.copy()
4 r = b - A @ x
5 r2 = np.dot(r, r)
6 norm_b = np.linalg.norm(b)
7 print(f"iter {0}, ||r|| = {np.sqrt(r2)}")
8
9 info = 0
10 for k in range(maxiter):
11 Ar = A @ r
12 alpha = r2 / np.dot(r, Ar)
13 x = x + alpha * r
14 r = r - alpha * Ar
15 r2 = np.dot(r, r)
16
17 print(f"iter {k+1}, ||r|| = {np.sqrt(r2)}")
18
19 if np.sqrt(r2) <= max(rtol*norm_b, atol):
20 info = 1
21 break
22
23 return x, info, k+1

1.4.1 Método do gradiente conjugado

O método do gradiente conjugado é uma melhoria do método do gradiente, que constrói as aproximações 𝒙(k) na forma

𝒙(k+1)=𝒙(k)+αk𝒑(k), (1.131)

onde 𝒑(k) é uma direção de busca conjugada as direções anteriores 𝒑(0),𝒑(1),,𝒑(k1). Isto é, as direções de busca satisfazem

(𝒑(i),A𝒑(j))=0, (1.132)

para todo ij. A primeira direção de busca é escolhida como o resíduo inicial, i.e. 𝒑(0)=𝒓(0). As demais direções são construídas como combinações lineares do resíduo atual e da direção de busca anterior, i.e.

𝒑(k)=𝒓(k)+βk𝒑(k1), (1.133)

onde βk é escolhido de modo a garantir a conjugação das direções. Pode-se mostrar (consulte [6]) que a escolha

βk=(𝒓(k),𝒓(k))(𝒓(k1),𝒓(k1)). (1.134)

garante a conjugação das direções de busca. O tamanho do passo αk é escolhido como no método do gradiente, i.e.

αk=(𝒓(k),𝒓(k))(𝒑(k),A𝒑(k)). (1.135)

Em resumo, a iteração do método do gradiente conjugado é dada por

𝒙(k+1)=𝒙(k)+αk𝒑(k), (1.136)
𝒓(k+1)=𝒓(k)αkA𝒑(k), (1.137)
𝒑(k+1)=𝒓(k+1)+βk𝒑(k), (1.138)
ondeαk=(𝒓(k),𝒓(k))(𝒑(k),A𝒑(k)), (1.139)
βk=(𝒓(k+1),𝒓(k+1))(𝒓(k),𝒓(k)). (1.140)
Código 7: mgc.py
1import numpy as np
2
3def mgc(A, b, x0, rtol=1e-5, atol=0.0, maxiter=100000):
4 x = x0.copy()
5 r = b - A @ x
6 p = r.copy()
7 r2 = np.dot(r, r)
8 norm_b = np.linalg.norm(b)
9
10 info = 0
11 for k in range(maxiter):
12 Ap = A @ p
13 alpha = r2 / np.dot(p, Ap)
14 x = x + alpha * p
15 r = r - alpha * Ap
16 r2_new = np.dot(r, r)
17
18 if np.sqrt(r2_new) <= max(rtol*norm_b, atol):
19 info = 1
20 break
21
22 beta = r2_new / r2
23 p = r + beta * p
24 r2 = r2_new
25
26 return x, info, k+1
Exemplo 1.4.1.(Equação de Poisson 2D)

Considere o problema de Poisson2222endnote: 22Siméon Denis Poisson, 1781 - 1840, matemático francês. Fonte: Wikipédia:Siméon Denis Poisson. 2D

Δu=f,em Ω=(0,1)×(0,1), (1.141)
u=0,em Ω, (1.142)

onde a fonte é dada por

f(x)={100,se 0.6x,y0.7,0,caso contrário. (1.143)

Discretizando o problema com diferenças finitas centrais em uma malha uniforme de n×n pontos, obtemos um sistema linear A𝒖=𝒃, onde A é uma matriz esparsa simétrica positiva definida de ordem N=n2 e 𝒖 é o vetor que contém os valores aproximados de u nos nodos da malha, 𝒃N.

A seguinte tabela mostra o número de iterações necessárias para a convergência do método do gradiente (MG) e do método do gradiente conjugado (MGC), para diferentes tamanhos de malha. O critério de parada é dado por rtol=105 e atol=0.

n MG MGC
11 203 11
21 808 48
41 3190 97
81 12702 196
Observação 1.4.2.(Aspectos de convergência)

Para A matriz simétrica positiva definida, o método do gradiente conjugado converge para a solução do sistema linear A𝒙=𝒃 em no máximo n iterações (desconsiderando-se erros de arredondamentos), onde n é a ordem da matriz A [1, Teorema 7.32].

1.4.2 Exercícios

E. 1.4.1.

Aplique o método do gradiente para resolver o sistema linear A𝒙=𝒃, com

A=[2113] (1.144)

e 𝒃=(3,4). Usando uma aproximação inicial 𝒙(0)=(0,0), faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno do funcional J(𝒙) e mostre as iterações do método do gradiente sobre o gráfico.

E. 1.4.2.

Aplique o método do gradiente conjugado para resolver o sistema linear A𝒙=𝒃, com

A=[2113] (1.145)

e 𝒃=(3,4). Usando uma aproximação inicial 𝒙(0)=(0,0), faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno do funcional J(𝒙) e mostre as iterações do método do gradiente sobre o gráfico.

E. 1.4.3.

Para A matriz simétrica positiva definida, mostre que, se

αk=minα>0J(𝒙(k)+α𝒑(k)), (1.146)

então,

αk=(𝒓(k),𝒓(k))(𝒑(k),A𝒑(k)). (1.147)
E. 1.4.4.

Considere a seguinte equação diferencial parcial de difusão-advecção com condições de contorno de Dirichlet homogêneas

Δu+𝒂u=8(xx2)(yy2),(x,y)(0,1)2. (1.148)

Aplique o método de diferenças finitas com fórmulas de diferenças centrais para obter um problema discreto A𝒖=𝒃 que aproxime a solução. Então, aplique o método do gradiente conjugado para resolver o sistema linear. Considere 𝒂=(1,1) e 𝒂=(10,10). Compare os resultados.

E. 1.4.5.

Use um esquema upwind para a discretização do problema de difusão-advecção do Exemplo 1.4.4. I.e., use a fórmula de diferenças central para a discretização do termo de difusão e a fórmula de diferenças progressiva para o termo de advecção. Verifique que o sistema linear resultante A𝒖=𝒃 não é simétrico. Entretanto, como A é não-singular, ATA é simétrica positiva definida. Assim, aplique os métodos do gradiente e do gradiente conjugado para resolver o sistema normal

ATA𝒖=AT𝒃. (1.149)

Analise a convergência dos métodos para diferentes tamanho de malha. Por fim, compare os resultados com os obtidos no Exemplo 1.4.4.


Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a política de uso de dados para mais informações.

Licença Creative Commons
Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Matemática Numérica III

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

1.4 Método do gradiente conjugado

O método do gradiente conjugado é um dos métodos mais eficientes para a resolução de sistemas lineares

A𝒙=𝒃, (1.119)

com A matriz simétrica positiva definida n×n e bn. Foi desenvolvido por Magnus Hestenes2020endnote: 20Magnus R. Hestenes, 1906 - 1991, matemático americano estadunidense. Fonte: Wikipedia. e Eduard Stiefel2121endnote: 21Eduard Stiefel, 1909 - 1978, matemático suíço. Fonte: Wikipédia. em 1952.

A ideia é encontrar a solução do sistema linear como o ponto de mínimo da função quadrática

J(𝒙)=12(𝒙,A𝒙)(𝒙,𝒃), (1.120)

onde (,) denota o produto interno usual em n.

Note que J(𝒙) é estritamente convexa, pois A é simétrica positiva definida. Assim, J(𝒙) possui um único ponto de mínimo global, que é a solução do sistema linear A𝒙=b. De fato, supondo 𝒙 solução do sistema, temos

J(𝒙)=12(𝒙,A𝒙)𝒙Tb+12(𝒙,A𝒙)(𝒙,b) (1.121)
=12(𝒙𝒙,A(𝒙𝒙))12(𝒙,A𝒙) (1.122)

Logo, como A é simétrica positiva definida, temos que J(𝒙) é mínimo se, e somente se, 𝒙𝒙=𝟎, i.e. 𝒙=𝒙.

Ainda, observamos que o gradiente de J(𝒙) é dado por

J(𝒙)=A𝒙𝒃, (1.123)

i.e. J(𝒙) é igual ao vetor oposto do resíduo do sistema linear

𝒓(𝒙)=𝒃A𝒙. (1.124)

A solução 𝒙 é o único vetor tal que J(𝒙)=𝟎. Mais ainda, lembrando que o gradiente aponta na direção de maior crescimento da função, temos que o resíduo aponta na direção de maior decrescimento de J(𝒙). Isto nos motiva a considerar o método do gradiente (ou método da descida mais íngrime), que consiste nas iterações

𝒙(k+1)=𝒙(k)+αk𝒓(k), (1.125)

onde 𝒓(k)=𝒃A𝒙(k) é o resíduo na iteração k e αk>0 é um tamanho do passo escolhido, com 𝒙(0) uma aproximação inicial da solução do sistema. Em outras palavras, o método do gradiente constrói uma sequência de aproximações (𝒙(k))k da solução do sistema, movendo-se na direção do resíduo do sistema linear.

Podemos escolher o tamanho do passo αk de modo a minimizar J(𝒙(k+1)) ao longo da direção do resíduo 𝒓(k). Isto é, escolhemos αk como

αk=minα>0J(𝒙(k)+α𝒓(k)), (1.126)

o que é conhecido como busca linear exata. Fazendo as contas, obtemos

αk=(𝒓(k),𝒓(k))(𝒓(k),A𝒓(k)). (1.127)

Com tudo isso, temos a iteração do método do gradiente

𝒙(k+1)=𝒙(k)+αk𝒓(k), (1.128)
ondeαk=(𝒓(k),𝒓(k))(𝒓(k),A𝒓(k)). (1.129)
Observação 1.4.1.(Detalhe de implementação)

Observamos que o cálculo de αk requer o produto A𝒓(k) e que o resíduo 𝒓(k) requer o produto A𝒙(k). Essas multiplicações matriz-vetor são os passos computacionalmente mais custosos da iteração. Assim, para economizar uma multiplicação por iteração, podemos atualizar o resíduo como

𝒓(k+1)=𝒓(k)αkA𝒓(k), (1.130)

evitando o cálculo de A𝒙(k+1).

Código 6: mg.py
1import numpy as np
2def mg(A, b, x0, rtol=1e-5, atol=0.0, maxiter=100000):
3 x = x0.copy()
4 r = b - A @ x
5 r2 = np.dot(r, r)
6 norm_b = np.linalg.norm(b)
7 print(f"iter {0}, ||r|| = {np.sqrt(r2)}")
8
9 info = 0
10 for k in range(maxiter):
11 Ar = A @ r
12 alpha = r2 / np.dot(r, Ar)
13 x = x + alpha * r
14 r = r - alpha * Ar
15 r2 = np.dot(r, r)
16
17 print(f"iter {k+1}, ||r|| = {np.sqrt(r2)}")
18
19 if np.sqrt(r2) <= max(rtol*norm_b, atol):
20 info = 1
21 break
22
23 return x, info, k+1

1.4.1 Método do gradiente conjugado

O método do gradiente conjugado é uma melhoria do método do gradiente, que constrói as aproximações 𝒙(k) na forma

𝒙(k+1)=𝒙(k)+αk𝒑(k), (1.131)

onde 𝒑(k) é uma direção de busca conjugada as direções anteriores 𝒑(0),𝒑(1),,𝒑(k1). Isto é, as direções de busca satisfazem

(𝒑(i),A𝒑(j))=0, (1.132)

para todo ij. A primeira direção de busca é escolhida como o resíduo inicial, i.e. 𝒑(0)=𝒓(0). As demais direções são construídas como combinações lineares do resíduo atual e da direção de busca anterior, i.e.

𝒑(k)=𝒓(k)+βk𝒑(k1), (1.133)

onde βk é escolhido de modo a garantir a conjugação das direções. Pode-se mostrar (consulte [6]) que a escolha

βk=(𝒓(k),𝒓(k))(𝒓(k1),𝒓(k1)). (1.134)

garante a conjugação das direções de busca. O tamanho do passo αk é escolhido como no método do gradiente, i.e.

αk=(𝒓(k),𝒓(k))(𝒑(k),A𝒑(k)). (1.135)

Em resumo, a iteração do método do gradiente conjugado é dada por

𝒙(k+1)=𝒙(k)+αk𝒑(k), (1.136)
𝒓(k+1)=𝒓(k)αkA𝒑(k), (1.137)
𝒑(k+1)=𝒓(k+1)+βk𝒑(k), (1.138)
ondeαk=(𝒓(k),𝒓(k))(𝒑(k),A𝒑(k)), (1.139)
βk=(𝒓(k+1),𝒓(k+1))(𝒓(k),𝒓(k)). (1.140)
Código 7: mgc.py
1import numpy as np
2
3def mgc(A, b, x0, rtol=1e-5, atol=0.0, maxiter=100000):
4 x = x0.copy()
5 r = b - A @ x
6 p = r.copy()
7 r2 = np.dot(r, r)
8 norm_b = np.linalg.norm(b)
9
10 info = 0
11 for k in range(maxiter):
12 Ap = A @ p
13 alpha = r2 / np.dot(p, Ap)
14 x = x + alpha * p
15 r = r - alpha * Ap
16 r2_new = np.dot(r, r)
17
18 if np.sqrt(r2_new) <= max(rtol*norm_b, atol):
19 info = 1
20 break
21
22 beta = r2_new / r2
23 p = r + beta * p
24 r2 = r2_new
25
26 return x, info, k+1
Exemplo 1.4.1.(Equação de Poisson 2D)

Considere o problema de Poisson2222endnote: 22Siméon Denis Poisson, 1781 - 1840, matemático francês. Fonte: Wikipédia:Siméon Denis Poisson. 2D

Δu=f,em Ω=(0,1)×(0,1), (1.141)
u=0,em Ω, (1.142)

onde a fonte é dada por

f(x)={100,se 0.6x,y0.7,0,caso contrário. (1.143)

Discretizando o problema com diferenças finitas centrais em uma malha uniforme de n×n pontos, obtemos um sistema linear A𝒖=𝒃, onde A é uma matriz esparsa simétrica positiva definida de ordem N=n2 e 𝒖 é o vetor que contém os valores aproximados de u nos nodos da malha, 𝒃N.

A seguinte tabela mostra o número de iterações necessárias para a convergência do método do gradiente (MG) e do método do gradiente conjugado (MGC), para diferentes tamanhos de malha. O critério de parada é dado por rtol=105 e atol=0.

n MG MGC
11 203 11
21 808 48
41 3190 97
81 12702 196
Observação 1.4.2.(Aspectos de convergência)

Para A matriz simétrica positiva definida, o método do gradiente conjugado converge para a solução do sistema linear A𝒙=𝒃 em no máximo n iterações (desconsiderando-se erros de arredondamentos), onde n é a ordem da matriz A [1, Teorema 7.32].

1.4.2 Exercícios

E. 1.4.1.

Aplique o método do gradiente para resolver o sistema linear A𝒙=𝒃, com

A=[2113] (1.144)

e 𝒃=(3,4). Usando uma aproximação inicial 𝒙(0)=(0,0), faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno do funcional J(𝒙) e mostre as iterações do método do gradiente sobre o gráfico.

E. 1.4.2.

Aplique o método do gradiente conjugado para resolver o sistema linear A𝒙=𝒃, com

A=[2113] (1.145)

e 𝒃=(3,4). Usando uma aproximação inicial 𝒙(0)=(0,0), faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno do funcional J(𝒙) e mostre as iterações do método do gradiente sobre o gráfico.

E. 1.4.3.

Para A matriz simétrica positiva definida, mostre que, se

αk=minα>0J(𝒙(k)+α𝒑(k)), (1.146)

então,

αk=(𝒓(k),𝒓(k))(𝒑(k),A𝒑(k)). (1.147)
E. 1.4.4.

Considere a seguinte equação diferencial parcial de difusão-advecção com condições de contorno de Dirichlet homogêneas

Δu+𝒂u=8(xx2)(yy2),(x,y)(0,1)2. (1.148)

Aplique o método de diferenças finitas com fórmulas de diferenças centrais para obter um problema discreto A𝒖=𝒃 que aproxime a solução. Então, aplique o método do gradiente conjugado para resolver o sistema linear. Considere 𝒂=(1,1) e 𝒂=(10,10). Compare os resultados.

E. 1.4.5.

Use um esquema upwind para a discretização do problema de difusão-advecção do Exemplo 1.4.4. I.e., use a fórmula de diferenças central para a discretização do termo de difusão e a fórmula de diferenças progressiva para o termo de advecção. Verifique que o sistema linear resultante A𝒖=𝒃 não é simétrico. Entretanto, como A é não-singular, ATA é simétrica positiva definida. Assim, aplique os métodos do gradiente e do gradiente conjugado para resolver o sistema normal

ATA𝒖=AT𝒃. (1.149)

Analise a convergência dos métodos para diferentes tamanho de malha. Por fim, compare os resultados com os obtidos no Exemplo 1.4.4.


Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a política de uso de dados para mais informações.

Licença Creative Commons
Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Pedro H A Konzen
| | | |