Matemática Numérica III

1 Sistemas Lineares 1.3 Métodos iterativos básicos 1.5 Métodos de projeção

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

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A\boldsymbol{% x}=\boldsymbol{b},

(1.119)

com $A$ matriz simétrica positiva definida $n\times n$ e $b\in\mathbb{R}^{n}$ . Foi desenvolvido por Magnus Hestenes²⁰²⁰endnote: ²⁰Magnus R. Hestenes, 1906 - 1991, matemático americano estadunidense. Fonte: Wikipedia. e Eduard Stiefel²¹²¹endnote: ²¹Eduard 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

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}J(\boldsymbol% {x})=\frac{1}{2}(\boldsymbol{x},A\boldsymbol{x})-(\boldsymbol{x},\boldsymbol{b% }),

(1.120)

onde $(\cdot,\cdot)$ denota o produto interno usual em $\mathbb{R}^{n}$ .

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

	$\displaystyle J(\boldsymbol{x})=\frac{1}{2}(\boldsymbol{x},A\boldsymbol{x})-% \boldsymbol{x}^{T}b+\frac{1}{2}(\boldsymbol{x}^{},A\boldsymbol{x}^{})-\left(% \boldsymbol{x}^{*},b\right)$		(1.121)
	$\displaystyle\text{}\quad=\frac{1}{2}(\boldsymbol{x}-\boldsymbol{x}^{},A(% \boldsymbol{x}-\boldsymbol{x}^{}))-\frac{1}{2}\left(\boldsymbol{x}^{},A% \boldsymbol{x}^{}\right)$		(1.122)

Logo, como $A$ é simétrica positiva definida, temos que $J(\boldsymbol{x})$ é mínimo se, e somente se, $\boldsymbol{x}-\boldsymbol{x}^{*}=\boldsymbol{0}$ , i.e. $\boldsymbol{x}=\boldsymbol{x}^{*}$ .

Ainda, observamos que o gradiente de $J(\boldsymbol{x})$ é dado por

\nabla J(\boldsymbol{x})=A\boldsymbol{x}-\boldsymbol{b},

(1.123)

i.e. $\nabla J(\boldsymbol{x})$ é igual ao vetor oposto do resíduo do sistema linear

\boldsymbol{r}(\boldsymbol{x})=\boldsymbol{b}-A\boldsymbol{x}.

(1.124)

A solução $\boldsymbol{x}^{*}$ é o único vetor tal que $\nabla J(\boldsymbol{x}^{*})=\boldsymbol{0}$ . 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(\boldsymbol{x})$ . Isto nos motiva a considerar o método do gradiente (ou método da descida mais íngrime), que consiste nas iterações

\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\alpha_{k}\boldsymbol{r}^{(k)},

(1.125)

onde $\boldsymbol{r}^{(k)}=\boldsymbol{b}-A\boldsymbol{x}^{(k)}$ é o resíduo na iteração $k$ e $\alpha_{k}>0$ é um tamanho do passo escolhido, com $\boldsymbol{x}^{(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 $\left(\boldsymbol{x}^{(k)}\right)_{k\in\mathbb{N}}$ da solução do sistema, movendo-se na direção do resíduo do sistema linear.

Podemos escolher o tamanho do passo $\alpha_{k}$ de modo a minimizar $J(\boldsymbol{x}^{(k+1)})$ ao longo da direção do resíduo $\boldsymbol{r}^{(k)}$ . Isto é, escolhemos $\alpha_{k}$ como

\alpha_{k}=\min_{\alpha>0}J(\boldsymbol{x}^{(k)}+\alpha\boldsymbol{r}^{(k)}),

(1.126)

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

\alpha_{k}=\frac{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{r}^% {(k)},A\boldsymbol{r}^{(k)})}.

(1.127)

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

	$\displaystyle\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\alpha_{k}\boldsymbol% {r}^{(k)},$		(1.128)
	$\displaystyle\mbox{onde}\quad\alpha_{k}=\frac{(\boldsymbol{r}^{(k)},% \boldsymbol{r}^{(k)})}{(\boldsymbol{r}^{(k)},A\boldsymbol{r}^{(k)})}.$		(1.129)

Observação 1.4.1.(Detalhe de implementação)

Observamos que o cálculo de $\alpha_{k}$ requer o produto $A\boldsymbol{r}^{(k)}$ e que o resíduo $\boldsymbol{r}^{(k)}$ requer o produto $A\boldsymbol{x}^{(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

\boldsymbol{r}^{(k+1)}=\boldsymbol{r}^{(k)}-\alpha_{k}A\boldsymbol{r}^{(k)},

(1.130)

evitando o cálculo de $A\boldsymbol{x}^{(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)}")

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)

17 print(f"iter {k+1}, ||r|| = {np.sqrt(r2)}")

19 if np.sqrt(r2) <= max(rtol*norm_b, atol):

20 info = 1

21 break

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 $\boldsymbol{x}^{(k)}$ na forma

\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\alpha_{k}\boldsymbol{p}^{(k)},

(1.131)

onde $\boldsymbol{p}^{(k)}$ é uma direção de busca conjugada as direções anteriores $\boldsymbol{p}^{(0)},\boldsymbol{p}^{(1)},\cdots,\boldsymbol{p}^{(k-1)}$ . Isto é, as direções de busca satisfazem

(\boldsymbol{p}^{(i)},A\boldsymbol{p}^{(j)})=0,

(1.132)

para todo $i\neq j$ . A primeira direção de busca é escolhida como o resíduo inicial, i.e. $\boldsymbol{p}^{(0)}=\boldsymbol{r}^{(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.

\boldsymbol{p}^{(k)}=\boldsymbol{r}^{(k)}+\beta_{k}\boldsymbol{p}^{(k-1)},

(1.133)

onde $\beta_{k}$ é escolhido de modo a garantir a conjugação das direções. Pode-se mostrar (consulte [9]) que a escolha

\beta_{k}=\frac{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{r}^{% (k-1)},\boldsymbol{r}^{(k-1)})}.

(1.134)

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

\alpha_{k}=\frac{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{p}^% {(k)},A\boldsymbol{p}^{(k)})}.

(1.135)

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

	$\displaystyle\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\alpha_{k}\boldsymbol% {p}^{(k)},$		(1.136)
	$\displaystyle\boldsymbol{r}^{(k+1)}=\boldsymbol{r}^{(k)}-\alpha_{k}A% \boldsymbol{p}^{(k)},$		(1.137)
	$\displaystyle\boldsymbol{p}^{(k+1)}=\boldsymbol{r}^{(k+1)}+\beta_{k}% \boldsymbol{p}^{(k)},$		(1.138)
	$\displaystyle\mbox{onde}\quad\alpha_{k}=\frac{(\boldsymbol{r}^{(k)},% \boldsymbol{r}^{(k)})}{(\boldsymbol{p}^{(k)},A\boldsymbol{p}^{(k)})},$		(1.139)
	$\displaystyle\text{}\quad\quad\beta_{k}=\frac{(\boldsymbol{r}^{(k+1)},% \boldsymbol{r}^{(k+1)})}{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}.$		(1.140)

Código 7: mgc.py

⬇

1import numpy as np

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)

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)

18 if np.sqrt(r2_new) <= max(rtol*norm_b, atol):

19 info = 1

20 break

22 beta = r2_new / r2

23 p = r + beta * p

24 r2 = r2_new

26 return x, info, k+1

Exemplo 1.4.1.(Equação de Poisson 2D)

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

	$\displaystyle-\Delta u=f,\quad\text{em }\Omega=(0,1)\times(0,1),$		(1.141)
	$\displaystyle u=0,\quad\text{em }\partial\Omega,$		(1.142)

onde a fonte é dada por

f(x)=\begin{cases}100,&\text{se }0.6\leq x,y\leq 0.7,\\ 0,&\text{caso contrário}.\end{cases}

(1.143)

Discretizando o problema com diferenças finitas centrais em uma malha uniforme de $n\times n$ pontos, obtemos um sistema linear $A\boldsymbol{u}=\boldsymbol{b}$ , onde $A$ é uma matriz esparsa simétrica positiva definida de ordem $N=n^{2}$ e $\boldsymbol{u}$ é o vetor que contém os valores aproximados de $u$ nos nodos da malha, $\boldsymbol{b}\in\mathbb{R}^{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 $\textrm{rtol}=10^{-5}$ e $\textrm{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\boldsymbol{x}=\boldsymbol{b}$ em no máximo $n$ iterações (desconsiderando-se erros de arredondamentos), onde $n$ é a ordem da matriz $A$ [3, Teorema 7.32].

1.4.2 Exercícios

E. 1.4.1.

Aplique o método do gradiente para resolver o sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ , com

A=\begin{bmatrix}2&1\\ 1&3\end{bmatrix}

(1.144)

e $\boldsymbol{b}=(3,4)$ . Usando uma aproximação inicial $\boldsymbol{x}^{(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(\boldsymbol{x})$ 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\boldsymbol{x}=\boldsymbol{b}$ , com

A=\begin{bmatrix}2&1\\ 1&3\end{bmatrix}

(1.145)

e $\boldsymbol{b}=(3,4)$ . Usando uma aproximação inicial $\boldsymbol{x}^{(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(\boldsymbol{x})$ 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

\alpha_{k}=\min_{\alpha>0}J(\boldsymbol{x}^{(k)}+\alpha\boldsymbol{p}^{(k)}),

(1.146)

então,

\alpha_{k}=\frac{(\boldsymbol{r}^{(k)},\boldsymbol{r}^{(k)})}{(\boldsymbol{p}^% {(k)},A\boldsymbol{p}^{(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

-\Delta u+\boldsymbol{a}\cdot\nabla u=8(x-x^{2})(y-y^{2}),\leavevmode\nobreak% \ (x,y)\in(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\boldsymbol{u}=\boldsymbol{b}$ que aproxime a solução. Então, aplique o método do gradiente conjugado para resolver o sistema linear. Considere $\boldsymbol{a}=(1,1)$ e $\boldsymbol{a}=(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\boldsymbol{u}=\boldsymbol{b}$ não é simétrico. Entretanto, como $A$ é não-singular, $A^{T}A$ é simétrica positiva definida. Assim, aplique os métodos do gradiente e do gradiente conjugado para resolver o sistema normal

A^{T}A\boldsymbol{u}=A^{T}\boldsymbol{b}.

(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!

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.

Política de uso de dados