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\bm{x}=\bm{b},

(1.124)

com $A$ matriz $n\times n$ simétrica positiva definida e $\bm{b}\in\mathbb{R}^{n}$ . Foi desenvolvido por Magnus Hestenes²⁰²⁰20Magnus R. Hestenes, 1906 - 1991, matemático americano estadunidense. Fonte: Wikipédia: Magnus Hestenes. e Eduard Stiefel²¹²¹21Eduard Stiefel, 1909 - 1978, matemático suíço. Fonte: Wikipédia: Eduard Stiefel. 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(\bm{x})=% \frac{1}{2}(\bm{x},A\bm{x})-(\bm{x},\bm{b}),

(1.125)

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

Note que $J(\bm{x})$ é estritamente convexa, pois $A$ é simétrica positiva definida. Assim, $J(\bm{x})$ possui um único ponto de mínimo global, que é a solução do sistema linear $A\bm{x}=b$ . De fato, seja $\bm{x}$ e $\bm{v}\neq 0$ dois vetores quaisquer e $t\in\mathbb{R}$ , temos

	$\displaystyle J(\bm{x}+t\bm{v})=\frac{1}{2}(\bm{x}+t\bm{v},A(\bm{x}+t\bm{v}))-% (\bm{x}+t\bm{v},\bm{b})$		(1.126)
	$\displaystyle\text{}\quad=\frac{1}{2}(\bm{x},A\bm{x})+\frac{t}{2}(\bm{x},A\bm{% v})+\frac{t}{2}(\bm{v},A\bm{x})+\frac{t^{2}}{2}(\bm{v},A\bm{v})-(\bm{x},\bm{b}% )-t(\bm{v},\bm{b})$		(1.127)
	$\displaystyle\text{}\quad=\frac{1}{2}(\bm{x},A\bm{x})-(\bm{x},\bm{b})+t(\bm{v}% ,A\bm{x}-\bm{b})+\frac{t^{2}}{2}(\bm{v},A\bm{v})$		(1.128)

donde, obtemos que

\displaystyle J(\bm{x}+t\bm{v})=J(\bm{x})-t(\bm{v},\bm{b}-A\bm{x})+\frac{t^{2}% }{2}(\bm{v},A\bm{v})

(1.129)

Para $\bm{x}$ e $\bm{v}$ fixos, $J(\bm{x}+t\bm{v})$ é uma função quadrática de $t$ , que atinge seu mínimo no ponto crítico

	$\displaystyle\frac{d}{dt}J(\bm{x}+t\bm{v})=-(\bm{v},\bm{b}-A\bm{x})+t(\bm{v},A% \bm{v})=0$		(1.130)
	$\displaystyle\hat{t}=\frac{(\bm{v},\bm{b}-A\bm{x})}{(\bm{v},A\bm{v})}.$		(1.131)

Substituindo este valor de $t$ em (1.129), obtemos

\displaystyle J(\bm{x}+\hat{t}\bm{v})=J(\bm{x})-\frac{(\bm{v},\bm{b}-A\bm{x})^% {2}}{2(\bm{v},A\bm{v})}\leq J(\bm{x}),

(1.132)

Assim, se $\bm{x}^{*}$ é a solução do sistema linear, temos que $\bm{b}-A\bm{x}^{*}=\bm{0}$ e, portanto,

\displaystyle J(\bm{x}^{*}+\hat{t}\bm{v})=J(\bm{x}^{*}),

(1.133)

i.e. $\bm{x}^{*}$ é o único ponto de mínimo global de $J(\bm{x})$ .

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

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

(1.134)

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

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

(1.135)

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

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

(1.136)

onde $\bm{r}^{(k)}=\bm{b}-A\bm{x}^{(k)}$ é o resíduo na iteração $k$ e $\alpha_{k}>0$ é um tamanho do passo escolhido, com $\bm{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(\bm{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(\bm{x}^{(k+1)})$ ao longo da direção do resíduo $\bm{r}^{(k)}$ . Isto é, escolhemos $\alpha_{k}$ como

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

(1.137)

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

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

(1.138)

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

	$\displaystyle\bm{x}^{(k+1)}=\bm{x}^{(k)}+\alpha_{k}\bm{r}^{(k)},$		(1.139)
	$\displaystyle\mbox{onde}\quad\alpha_{k}=\frac{(\bm{r}^{(k)},\bm{r}^{(k)})}{(% \bm{r}^{(k)},A\bm{r}^{(k)})}.$		(1.140)

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

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

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

(1.141)

evitando o cálculo de $A\bm{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 $\bm{x}^{(k)}$ na forma

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

(1.142)

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

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

(1.143)

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

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

(1.144)

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

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

(1.145)

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{(\bm{r}^{(k)},\bm{r}^{(k)})}{(\bm{p}^{(k)},A\bm{p}^{(k)})}.

(1.146)

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

	$\displaystyle\bm{x}^{(k+1)}=\bm{x}^{(k)}+\alpha_{k}\bm{p}^{(k)},$		(1.147)
	$\displaystyle\bm{r}^{(k+1)}=\bm{r}^{(k)}-\alpha_{k}A\bm{p}^{(k)},$		(1.148)
	$\displaystyle\bm{p}^{(k+1)}=\bm{r}^{(k+1)}+\beta_{k}\bm{p}^{(k)},$		(1.149)
	$\displaystyle\mbox{onde}\quad\alpha_{k}=\frac{(\bm{r}^{(k)},\bm{r}^{(k)})}{(% \bm{p}^{(k)},A\bm{p}^{(k)})},$		(1.150)
	$\displaystyle\text{}\quad\quad\beta_{k}=\frac{(\bm{r}^{(k+1)},\bm{r}^{(k+1)})}% {(\bm{r}^{(k)},\bm{r}^{(k)})}.$		(1.151)

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²²²²22Simé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.152)
	$\displaystyle u=0,\quad\text{em }\partial\Omega,$		(1.153)

onde a fonte é dada por

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

(1.154)

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

1.4.2 Exercícios

E. 1.4.1.

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

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

(1.155)

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

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

(1.156)

e $\bm{b}=(3,4)$ . Usando uma aproximação inicial $\bm{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(\bm{x})$ e mostre as iterações do método do gradiente conjugado sobre o gráfico.

E. 1.4.3.

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

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

(1.157)

então,

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

(1.158)

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+\bm{a}\cdot\nabla u=8(x-x^{2})(y-y^{2}),\leavevmode\nobreak\ (x,y)% \in(0,1)^{2}.

(1.159)

Aplique o método de diferenças finitas com fórmulas de diferenças centrais para obter um problema discreto $A\bm{u}=\bm{b}$ que aproxime a solução. Então, aplique o método do gradiente conjugado para resolver o sistema linear. Considere $\bm{a}=(1,1)$ e $\bm{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\bm{u}=\bm{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\bm{u}=A^{T}\bm{b}.

(1.160)

Analise a convergência dos métodos para diferentes tamanhos 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