Matemática numérica I

3 Métodos para sistemas lineares 3.3 Métodos de Jacobi e de Gauss-Seidel 3.5 Método do gradiente conjugado

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

3.4 Método do gradiente

Começamos observando que se $A$ é uma matriz $n\times n$ positiva definida¹¹¹¹11 $A$ é simétrica e $\bm{x}^{T}A\bm{x}>0$ para todo $\bm{x}\neq 0$ ., temos que $\bm{x}\in\mathbb{R}^{n}$ é solução de

A\bm{x}=\bm{b}

(3.293)

se, e somente se, $\bm{x}$ é solução do seguinte problema de minimização

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\min_{\bm{y}% \in\mathbb{R}^{n}}\left\{f(\bm{y}):=\frac{1}{2}\bm{y}^{T}A\bm{y}-\bm{y}^{T}\bm% {b}\right\}.

(3.294)

De fato, sejam $\bm{x}$ a solução de (3.293) e $\bm{y}$ a solução de (3.294), então

	$\displaystyle f(\bm{y}):=\frac{1}{2}\bm{y}^{T}A\bm{y}-\bm{y}^{T}\bm{b}$		(3.295)
	$\displaystyle\text{}\quad=\frac{1}{2}\bm{y}^{T}A\bm{y}-\bm{y}^{T}\bm{b}+{% \color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\frac{1}{2}% \bm{x}^{T}A\bm{x}-+\frac{1}{2}\bm{x}^{T}A\bm{x}}$		(3.296)
	$\displaystyle\text{}\quad=\frac{1}{2}(\bm{x}-\bm{y})^{T}A(\bm{x}-\bm{y})-\frac% {1}{2}\bm{x}^{T}A\bm{x}$		(3.297)

O segundo termo é independente de $\bm{y}$ e, como $A$ é positiva definida, temos

\frac{1}{2}(\bm{x}-\bm{y})^{T}A(\bm{x}-\bm{y})\geq 0.

(3.298)

Logo, o mínimo de $f$ ocorre quando $\bm{x}-\bm{y}=\bm{0}$ , i.e. $\bm{y}=\bm{x}$ .

A iteração do Método do Gradiente tem a forma

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \bm{x}^{(0)}=\text{aprox. inicial},$		(3.299)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \bm{x}^{(k+1)}=\bm{x}^{(k)}+\alpha_{k}\bm{d}^{(k)},$		(3.300)

para $k=0,1,2,\ldots$ , onde $\alpha_{k}>0$ é o tamanho do passo e $\bm{d}^{(k)}$ é o vetor direção de busca.

Para escolhermos a direção $\bm{d}^{(k)}$ , tomamos a fórmula de Taylor de $f$ em torno da aproximação $\bm{x}^{(k)}$

f\left(\bm{x}^{(k+1)}\right)=f\left(\bm{x}^{(k)}\right)+\alpha^{(k)}\nabla f% \left(\bm{x}^{(k)}\right)\cdot\bm{d}^{(k)}+\epsilon,

(3.301)

onde $\nabla f$ denota o gradiente de $f$ , i.e.

	$\displaystyle\nabla f(\bm{x})$	$\displaystyle:=\left(\frac{\partial f}{\partial x_{1}}(\bm{x}),\frac{\partial f% }{\partial x_{2}}(\bm{x}),\dotsc,\frac{\partial f}{\partial x_{n}}(\bm{x})\right)$		(3.302)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \nabla f(\bm{x})$	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% =A\bm{x}^{(k)}-\bm{b}.$		(3.303)

De , segue que se

\nabla f\left(\bm{x}^{(k)}\right)\cdot\bm{d}^{(k)}<0,

(3.304)

então $f\left(\bm{x}^{(k+1)}\right)<f\left(\bm{x}^{(k)}\right)$ , para $\alpha^{(k)}$ suficientemente pequeno. Em particular, podemos escolher

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{d}^{(k)}=% -\nabla f\left(\bm{x}^{(k)}\right),

(3.305)

se $\nabla f(\bm{x}^{(k)})\neq 0$ .

Refer to caption — Figura 3.1: Iterações do Método do Gradiente para sistemas lineares. Linha: $\|b-Ax\|$ . Pontos: $\bm{x}^{(k)}$ .

Do exposto acima, temos a iteração do Método do Gradiente

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \bm{x}^{(0)}=\text{aprox. inicial}$		(3.306)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \bm{x}^{(k+1)}=\bm{x}^{(k)}+\alpha^{(k)}\bm{r}^{(k)}$		(3.307)

com $k=0,1,2,\ldots$ , onde $\bm{r}^{(k)}$ é o resíduo

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{r}^{(k)}=% \bm{b}-A\bm{x}^{(k)}.

(3.308)

Exemplo 3.4.1.

Consideramos o sistema $Ax=b$ com

	$\displaystyle A=\begin{bmatrix}2&-1&0&0\\ -1&2&-1&0\\ 0&-1&2&-1\\ 0&0&-1&2\end{bmatrix},$		(3.309)
	$\displaystyle b=\begin{bmatrix}-3\\ 2\\ 2\\ -3\end{bmatrix}.$		(3.310)

Na Tabela 3.3 temos os resultados do emprego do método do gradiente com $\bm{x}^{(1)}=(0,0,0,0)$ e com passo constante $\alpha^{(k)}\equiv 0.5$ .

Tabela 3.3: Resultados referentes ao Exemplo 3.4.1.

k	$\bm{x}^{(k)}$	$\\|A\bm{x}^{(k)}-\bm{b}\\|$
0	$(0.0,0.0,0.0,0.0)$	$5.1\mathrm{e}\!+\!0$
1	$(-1.5,1.0,1.0,-1.5)$	$1.6\mathrm{e}\!+\!0$
2	$(-1.0,0.8,0.8,-1.0)$	$5.0\mathrm{e}\!-\!1$
3	$(-1.1,0.9,0.9,-1.1)$	$1.8\mathrm{e}\!-\!1$
4	$(-1.1,0.9,0.9,-1.1)$	$8.8\mathrm{e}\!-\!2$
5	$(-1.1,0.9,0.9,-1.1)$	$6.2\mathrm{e}\!-\!2$
6	$(-1.0,0.9,0.9,-1.0)$	$4.9\mathrm{e}\!-\!2$
7	$(-1.0,0.9,0.9,-1.0)$	$4.0\mathrm{e}\!-\!2$
8	$(-1.0,0.9,0.9,-1.0)$	$3.2\mathrm{e}\!-\!2$
9	$(-1.0,1.0,1.0,-1.0)$	$2.6\mathrm{e}\!-\!2$
10	$(-1.0,1.0,1.0,-1.0)$	$2.1\mathrm{e}\!-\!2$

Código 10: mg.py

⬇

1import numpy as np

2import numpy.linalg as npla

4def mg(A, b, x0, alpha=1e-2,

5 maxiter=100, atol=1.49e-8, rtol=1.49e-8):

7 n = b.size

8 x = x0.copy()

9 res = b - A@x

10 nres = npla.norm(res)

11 print(f'0: {x}, nres = {nres:.1e}')

12 info = -1

13 for k in range(maxiter):

14 x = x + alpha*res

16 res = b - A@x

17 nres = npla.norm(res)

19 print(f'{k+1}: {x}, nres = {nres:.1e}')

21 if (nres <= max(atol, rtol*npla.norm(b))):

22 info = 0

23 break

25 return x, info

27# matriz coefs

28A = np.array([[2., -1., 0., 0.],

29 [-1., 2., -1., 0.],

30 [0., -1., 2., -1.],

31 [0., 0., -1., 2.]])

32# vetor constante

33b = np.array([-3., 2., 2., -3.])

35# aprox. inicial

36x0 = np.zeros_like(b)

38x, info = mg(A, b, x0, alpha=0.5)

3.4.1 Escolha do passo

Para a escolha do passo, podemos usar o Método da Pesquisa Linear. A ideia é escolher o passo $\alpha^{(k)}$ tal que

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}f\left(\bm{x}% ^{(k)}+\alpha^{(k)}\bm{r}^{(k)}\right)=\min_{\alpha>0}f\left(\bm{x}^{(k)}+% \alpha\bm{r}^{(k)}\right).

(3.311)

Observando que $f(\bm{x}^{(k)}+\alpha\bm{r}^{(k)})$ é função apenas de $\alpha$ , temos que seu mínimo ocorre em seu ponto crítico, i.e.

	$\displaystyle\frac{d}{d\alpha}f\left(\bm{x}^{(k)}+\alpha\bm{r}^{(k)}\right)=0$		(3.312)
	$\displaystyle\nabla f\left(\bm{x}^{(k+1)}\right)\cdot\frac{d}{d\alpha}\left(% \bm{x}^{(k)}+\alpha\bm{r}^{(k)}\right)=0$		(3.313)
	$\displaystyle\nabla f(\bm{x}^{(k+1)})\cdot\bm{r}^{(k)}=0$		(3.314)
	$\displaystyle\left[A\left(\bm{x}^{(k)}+\alpha^{(k)}\bm{r}^{(k)}\right)-b\right% ]\cdot\bm{r}^{(k)}=0$		(3.315)
	$\displaystyle\left(A\bm{x}^{(k)}-\bm{b}\right)\cdot\bm{r}^{(k)}+\alpha^{(k)}% \bm{r}^{(k)}\cdot A\bm{r}^{(k)}=0,$		(3.316)

donde

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\alpha^{(k)}=% \frac{\bm{r}^{(k)}\cdot\bm{r}^{(k)}}{\bm{r}^{(k)}\cdot A\bm{r}^{(k)}}.

(3.317)

Exemplo 3.4.2.

Consideramos o sistema $Ax=b$ com

	$\displaystyle A=\begin{bmatrix}2&-1&0&0\\ -1&2&-1&0\\ 0&-1&2&-1\\ 0&0&-1&2\end{bmatrix},$		(3.318)
	$\displaystyle b=\begin{bmatrix}-3\\ 2\\ 2\\ -3\end{bmatrix}.$		(3.319)

Na Tabela 3.4 temos os resultados do emprego do método do gradiente com $\bm{x}^{(1)}=(0,0,0,0)$ e com passo escolhido conforme (3.317).

Tabela 3.4: Resultados referentes ao Exemplo LABEL:cap_sislin_sec_metg:ex:metg_alpha.

k	$\bm{x}^{(k)}$	$\alpha^{(k)}$	$\\|A\bm{x}^{(k)}-\bm{b}\\|$
$0$	$(0.0,0.0,0.0,0.0)$	$3.8\mathrm{e}\!-\!1$	$5.1\mathrm{e}\!+\!0$
$1$	$(-1.1,0.8,0.8,-1.1)$	$2.6\mathrm{e}\!+\!0$	$1.5\mathrm{e}\!-\!1$
$2$	$(-1.0,1.0,1.0,-1.0)$	$3.8\mathrm{e}\!-\!1$	$3.0\mathrm{e}\!-\!2$
$3$	$(-1.0,1.0,1.0,-1.0)$	$2.6\mathrm{e}\!+\!0$	$8.8\mathrm{e}\!-\!4$
$4$	$(-1.0,1.0,1.0,-1.0)$	$3.8\mathrm{e}\!+\!0$	$1.8\mathrm{e}\!-\!4$
$5$	$(-1.0,1.0,1.0,-1.0)$	$2.6\mathrm{e}\!+\!0$	$5.2\mathrm{e}\!-\!6$

⬇

1import numpy as np

2import numpy.linalg as npla

4def mg_pl(A, b, x0,

5 maxiter=100, atol=1.49e-8, rtol=1.49e-8):

7 n = b.size

8 x = x0.copy()

9 res = b - A@x

10 alpha = np.dot(res, res)/np.dot(res, A@res)

11 nres = npla.norm(res)

12 print(f'0: {x}, alpha = {alpha}, nres = {nres:.1e}')

13 info = -1

14 for k in range(maxiter):

15 x = x + alpha*res

17 res = b - A@x

18 alpha = np.dot(res, res)/np.dot(res, A@res)

19 nres = npla.norm(res)

21 print(f'{k+1}: {x}, alpha = {alpha}, nres = {nres:.1e}')

23 if (nres <= max(atol, rtol*npla.norm(b))):

24 info = 0

25 break

27 return x, info

29# matriz coefs

30A = np.array([[2., -1., 0., 0.],

31 [-1., 2., -1., 0.],

32 [0., -1., 2., -1.],

33 [0., 0., -1., 2.]])

34# vetor constante

35b = np.array([-3., 2., 2., -3.])

37# aprox. inicial

38x0 = np.zeros_like(b)

40x, info = mg_pl(A, b, x0)

3.4.2 Exercícios

E. 3.4.1.

Considere o sistema linear $A\bm{x}=\bm{b}$ com

	$\displaystyle A=\begin{bmatrix}2&-1&0&0\\ -1&2&-1&0\\ 0&-1&2&-1\\ 0&0&-1&2\end{bmatrix},$		(3.320)
	$\displaystyle b=\begin{bmatrix}5\\ -7\\ 6\\ -1\end{bmatrix}.$		(3.321)

Por tentativa e erro, encontre um valor para $\alpha$ tal que o Método do Gradiente converge para solução do sistema em menos de $80$ iterações. Use

\bm{x}^{(0)}=\bm{0}

(3.322)

como aproximação inicial e assuma o critério de parada

\|\bm{r}\|\leq\max\{\texttt{tol},\texttt{tol}\|\bm{b}\|\},

(3.323)

onde $\bm{r}$ é o resíduo do sistema e $\texttt{tol}=1.49\mathrm{e}\!-\!8$ .

$alpha=4.9\mathrm{e}\!-\!1$ , $\bm{x}=\left(2.0,-1.0,3.0,1.0\right)$

E. 3.4.2.

Considere o sistema linear $A\bm{x}=\bm{b}$ com

	$\displaystyle A=\begin{bmatrix}2&-1&1&0\\ -1&2&-1&0\\ 1&-1&3&-1\\ 0&0&-1&2\end{bmatrix},$		(3.324)
	$\displaystyle b=\begin{bmatrix}-8\\ 9\\ -10\\ 5\end{bmatrix}.$		(3.325)

Por tentativa e erro, encontre um valor para $\alpha$ tal que o Método do Gradiente converge para solução do sistema em menos de $50$ iterações. Use

\bm{x}^{(0)}=\bm{0}

(3.326)

como aproximação inicial e assuma o critério de parada

\|\bm{r}\|\leq\max\{\texttt{tol},\texttt{tol}\|\bm{b}\|\},

(3.327)

onde $\bm{r}$ é o resíduo do sistema e $\texttt{tol}=1.49\mathrm{e}\!-\!8$ .

$alpha=3.6\mathrm{e}\!-\!1$ , $\bm{x}=\left(-2.0,3.0,-1.0,2.0\right)$

E. 3.4.3.

Considere o sistema linear dado no E.3.4.1. Utilizando a mesma aproximação inicial e tolerância, aplique o Método do Gradiente com Pesquisa Linear. Quantas iterações são necessárias até a convergência e qual o valor médio de $\alpha^{(k)}$ utilizado durante as iterações?

$75$ , $\overline{\alpha}=5.0\mathrm{e}\!-\!1$

E. 3.4.4.

Considere o sistema linear dado no E.3.4.2. Utilizando a mesma aproximação inicial e tolerância, aplique o Método do Gradiente com Pesquisa Linear. Quantas iterações são necessárias até a convergência e qual o valor médio de $\alpha^{(k)}$ utilizado durante as iterações?

$35$ , $\overline{\alpha}=3.7\mathrm{e}\!-\!1$

E. 3.4.5.

Considere o problema de Laplace

	$\displaystyle-$	$\displaystyle u_{xx}=2,\quad 0<x<1,$		(3.328)
		$\displaystyle u(0)=u(1)=0.$		(3.329)

A discretização pelo Método das Diferenças Finitas em uma malha uniforme $x_{i}=ih$ , $i=0.1,\dotsc,n$ , $h=1/(n-1)$ , leva ao seguinte sistema linear

	$\displaystyle u_{1}=0$		(3.330)
	$\displaystyle-\frac{1}{h^{2}}u_{i-1}+\frac{2}{h^{2}}u_{i}-\frac{1}{h^{2}}u_{i+% 1}=2,$		(3.331)
	$\displaystyle u_{n}=0$		(3.332)

onde $u_{i}\approx u(x_{i})$ . Com $u^{(0)}=\bm{0}$ , aplique o Método do Gradiente com Pesquisa Linear para computar a solução deste sistema quando $n=10,20,40,80$ . Quantas iterações são necessárias para obter-se a convergência do método com critério de convergência

\|\bm{r}\|\leq\texttt{tol},

(3.333)

onde, $\bm{r}:=\bm{b}-A\bm{x}$ é o resíduo e $\texttt{tol}=1e-4$ .

$n$	$k$
$10$	$214$
$20$	$918$
$40$	$3840$
$80$	$15910$

Análise numérica

E. 3.4.6.

Sendo $f$ dada em (3.294), mostre que

\nabla f\left(\bm{x}\right)=-\bm{r},

(3.334)

com, $\bm{r}:=b-A\bm{x}$ o resíduo do sistema $A\bm{x}=\bm{b}$ .

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