Matemática Numérica III

1 Sistemas Lineares 1.4 Método do gradiente conjugado 1.6 Método GMRES

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

1.5 Métodos de projeção

A maior parte dos métodos iterativos para a resolução de sistemas lineares podem ser vistos como métodos de projeção. A ideia básica é resolver o sistema linear

A\boldsymbol{x}=\boldsymbol{b}

(1.150)

buscando uma solução aproximada $\tilde{\boldsymbol{x}}\in\mathcal{K}\subseteq\mathbb{R}^{n}$ de dimensão $m\leq n$ . $\mathcal{K}$ é chamado de subespaço de busca. Uma forma de garantir que a aproximação $\tilde{\boldsymbol{x}}$ seja boa é impor as chamadas condições de Petrov²³²³23Georgi Iwanowitsch Petrov, 1912 - 1987, engenheiro soviético. Fonte: Wikipedia.-Galerkin²⁴²⁴24Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia., em que o resíduo

\tilde{\boldsymbol{r}}=\boldsymbol{b}-A\tilde{\boldsymbol{x}}

(1.151)

é ortogonal a um dado subespaço $\mathcal{L}\subseteq\mathbb{R}^{n}$ de dimensão $m$ . Os métodos de projeção são classificados em métodos ortogonais, quando $\mathcal{L}=\mathcal{K}$ , e métodos oblíquos, quando $\mathcal{L}\neq\mathcal{K}$ .

Sejam $A$ uma matriz $n\times n$ não singular, $\mathcal{K}$ e $\mathcal{L}$ dois subespaços de $\mathbb{R}^{n}$ de dimensão $m\leq n$ e $\boldsymbol{x}^{(0)}$ uma aproximação inicial da solução do sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ . O método de projeção consiste em encontrar uma aproximação $\tilde{\boldsymbol{x}}\in\mathcal{K}$ da solução do sistema, tal que

\tilde{\boldsymbol{x}}=\boldsymbol{x}^{(0)}+\boldsymbol{\delta},

(1.152)

com $\boldsymbol{\delta}\in\mathcal{K}$ , satisfazendo a condição de Petrov-Galerkin

\tilde{\boldsymbol{r}}=\boldsymbol{b}-A\tilde{\boldsymbol{x}}\perp\mathcal{L}.

(1.153)

ou, equivalentemente,

\left(\boldsymbol{r}^{(0)}-A\boldsymbol{\delta},\boldsymbol{w}\right)=0,% \forall\boldsymbol{w}\in\mathcal{L},

(1.154)

Em resumo, o método de projeção determina a solução aproximada $\boldsymbol{x}^{(m)}$ pelas iterações

	$\displaystyle\tilde{\boldsymbol{x}}=\boldsymbol{x}^{(0)}+\boldsymbol{\delta},% \mbox{com }\boldsymbol{\delta}\in\mathcal{K},$		(1.155)
	$\displaystyle\text{}\quad\left(\boldsymbol{r}^{(0)}-A\boldsymbol{\delta},% \boldsymbol{w}\right)=0,\forall\boldsymbol{w}\in\mathcal{L}.$		(1.156)

A seguir, vamos apresentar três métodos de projeção unidimensionais. Embora estes métodos não sejam geralmente usados na prática (por não serem eficientes), eles foram a base metodológica para o desenvolvimento de métodos mais eficientes, como o método GMRES.

1.5.1 Métodos de projeção unidimensionais

Métodos de projeção unidimensionais são aqueles em que $\dim(\mathcal{K})=\dim(\mathcal{L})=1$ . Assim, temos $\mathcal{K}=\operatorname{span}\{\boldsymbol{v}\}$ e $\mathcal{L}=\operatorname{span}\{\boldsymbol{w}\}$ , para dados vetores $v,w\in\mathbb{R}^{n}$ . Com isso, a aproximação $\boldsymbol{x}^{(m)}$ é dada por

\tilde{\boldsymbol{x}}=\boldsymbol{x}^{(0)}+\alpha\boldsymbol{v},

(1.157)

onde $\alpha\in\mathbb{R}$ é tal que

\left(\boldsymbol{r}^{(0)}-\alpha A\boldsymbol{v},\boldsymbol{w}\right)=0,

(1.158)

ou seja,

\alpha=\frac{(\boldsymbol{r}^{(0)},\boldsymbol{w})}{(A\boldsymbol{v},% \boldsymbol{w})}.

(1.159)

1.5.2 Iteração da descida mais íngrime

A iteração da descida mais íngrime (ou método do gradiente) é um método de projeção ortogonal em que $\mathcal{K}=\mathcal{L}=\operatorname{span}\{\boldsymbol{r}\}$ . Assim, a aproximação $\boldsymbol{x}^{(m)}$ é dada por

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{x% }^{(m)}=\boldsymbol{x}^{(m-1)}+\alpha_{m}\boldsymbol{r}^{(m-1)},

(1.160)

onde

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\alpha_{m}=% \frac{(\boldsymbol{r}^{(m-1)},\boldsymbol{r}^{(m-1)})}{(A\boldsymbol{r}^{(m-1)% },\boldsymbol{r}^{(m-1)})}.

(1.161)

Notemos que, para uma matriz $A$ simétrica positiva definida, o método da descida mais íngrime é equivalente ao método do gradiente discutido na Seção 1.4, observando que $\boldsymbol{r}=-\nabla f(\boldsymbol{x})$ , onde

	$\displaystyle f(\boldsymbol{x})$	$\displaystyle=\\|\boldsymbol{x}-\boldsymbol{x}^{*}\\|_{A}^{2}$		(1.162)
		$\displaystyle\text{}\quad=(\boldsymbol{x}-\boldsymbol{x}^{},A(\boldsymbol{x}-% \boldsymbol{x}^{})),$		(1.163)

onde $\boldsymbol{x}^{*}$ é a solução do sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ .

Para detalhes sobre a implementação, consulte o Código 6.

1.5.3 Iteração do mínimo resíduo

Para uma matriz apenas positiva definida (i.e. $A+A^{T}$ simétrica positiva definida), a iteração do mínimo resíduo é um método de projeção oblíquo em que $\mathcal{K}=\operatorname{span}\{\boldsymbol{r}\}$ e $\mathcal{L}=\operatorname{span}\{A\boldsymbol{r}\}$ . Assim, a aproximação $\boldsymbol{x}^{(m)}$ é dada por

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{x% }^{(m)}=\boldsymbol{x}^{(m-1)}+\alpha_{m}\boldsymbol{r}^{(m-1)},

(1.164)

onde

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\alpha_{m}=% \frac{(A\boldsymbol{r}^{(m-1)},\boldsymbol{r}^{(m-1)})}{(A\boldsymbol{r}^{(m-1% )},A\boldsymbol{r}^{(m-1)})}.

(1.165)

Aqui, cada iterada minimiza

f(\boldsymbol{x})=\|b-A\boldsymbol{x}\|^{2},

(1.166)

na direção do resíduo $\boldsymbol{r}$ .

Código 8: imr.py

⬇

1from numpy import dot

2from numpy.linalg import norm

4def imr(A, b, x0, rtol=1e-5, atol=0.0, maxiter=100000):

5 x = x0.copy()

6 r = b - A @ x

7 norm_b = norm(b)

8 r2 = dot(r, r)

10 info = 0

11 for k in range(maxiter):

12 Ar = A @ r

13 alpha = dot(Ar, r) / dot(Ar, Ar)

14 x = x + alpha * r

15 r = r - alpha * Ar

16 r2 = dot(r, r)

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

19 info = 1

20 break

22 return x, info, k+1

Exemplo 1.5.1.(Problema de difusão-advecção 2D)

Consideremos o seguinte problema de difusão-advecção 2D

	$\displaystyle-\epsilon\Delta u+\boldsymbol{a}\cdot\nabla u=f,\quad\text{em }% \Omega=(0,1)\times(0,1),$		(1.167)
	$\displaystyle u=0,\quad\text{em }\partial\Omega,$		(1.168)

onde $\epsilon=1$ é o coeficiente de difusão, $\boldsymbol{a}=(20,1)$ é o campo de advecção e a fonte é dada por

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

(1.169)

Assumimos uma malha espacial uniforme de $n\times n$ , com tamanho de malha $h=1/(n-1)$ . Denotamos $u_{i,j}\approx u(x_{i},y_{i})$ , onde $x_{i}=ih$ e $y_{j}=jh$ , para $i,j=1,2,\dotsc,n$ . Aplicando um espuema upwind para a advecção e diferenças finitas centrais para a difusão, obtemos o seguinte esquema de diferenças finitas

	$\displaystyle-\epsilon\left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}+\frac{u% _{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}\right)$		(1.170)
	$\displaystyle\text{}\quad+a_{1}\frac{u_{i,j}-u_{i-1,j}}{h}+a_{2}\frac{u_{i,j+1% }-u_{i,j}}{h}=f_{i,j},$		(1.171)

para $i,j=2,3,\dotsc,n-1$ , onde $f_{i,j}=f(x_{i},y_{j})$ . As condições de contorno são dadas por $u_{1,j}=u_{n,1}=u_{i,1}=u_{i,n}=0$ , para $i,j=1,2,\dotsc,n-1$ . Por fim, consideramos a enumeração dos nodos da malha $k=i-1+(j-2)(n-2)$ para obtermos o sistema linear

A\boldsymbol{u}=\boldsymbol{f},

(1.172)

onde $A$ é uma matriz esparsa de ordem $N=(n-2)^{2}$ , $\boldsymbol{u}\in\mathbb{R}^{N}$ é o vetor das incógnitas e $\boldsymbol{f}\in\mathbb{R}^{N}$ é o vetor da fonte.

Observando que $A$ é apenas positiva definida (i.e. $A+A^{T}$ é simétrica positiva definida), aplicamos a iteração do mínimo resíduo para resolver o sistema linear. A seguinte tabela mostra o número de iterações necessárias para a convergência do método, para diferentes tamanhos de malha. O critério de parada é dado por $\textrm{rtol}=10^{-5}$ e $\textrm{atol}=0$ . Verifique!

$n$	IMR
11	129
21	469
41	1761
81	6770

1.5.4 Iteração da descida mais íngrime do resíduo

A iteração da descida mais íngrime do resíduo é um método de projeção oblíquo em que $\mathcal{K}=\operatorname{span}\{A^{T}\boldsymbol{r}\}$ e $\mathcal{L}=\operatorname{span}\{AA^{T}\boldsymbol{r}\}$ . Assim, a aproximação $\boldsymbol{x}^{(m)}$ é dada por

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{x% }^{(m)}=\boldsymbol{x}^{(m-1)}+\alpha_{m}A^{T}\boldsymbol{r}^{(m-1)},

(1.173)

onde

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\alpha_{m}=% \frac{(A^{T}\boldsymbol{r}^{(m-1)},A^{T}\boldsymbol{r}^{(m-1)})}{(AA^{T}% \boldsymbol{r}^{(m-1)},AA^{T}\boldsymbol{r}^{(m-1)})}.

(1.174)

Notemos que esta iteração é equivalente à do mínimo resíduo aplicada às equações normais

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

(1.175)

Logo, se $A$ é não-singular, o método converge para a solução do sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ .

Código 9: idir.py

⬇

1def idir(A, b, x0, rtol=1e-5, atol=0.0, maxiter=100000):

2 x = x0.copy()

3 r = b - A @ x

4 norm_b = norm(b)

5 r2 = dot(r, r)

7 info = 0

8 for k in range(maxiter):

9 v = A.T @ r

10 Av = A @ v

11 alpha = dot(v, v) / dot(Av, Av)

12 x = x + alpha * r

13 r = r - alpha * Av

14 r2 = dot(r, r)

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

17 info = 1

18 break

20 return x, info, k+1

1.5.5 Exercícios

E. 1.5.1.

Aplique a iteração da descida mais íngrime para resolver o sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ , com

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

(1.176)

e $\boldsymbol{b}=(3,4)$ . Usando uma aproximação inicial $\boldsymbol{x}^{(0)}=(0.5,0)$ , faça uma análise geométrica das iterações.

Dica: faça um gráfico de contorno da norma $\|\boldsymbol{r}(\boldsymbol{x})\|_{2}=\|\boldsymbol{b}-A\boldsymbol{x}\|_{2}$ e mostre as iterações do método sobre o gráfico.

E. 1.5.2.

Aplique a iteração do mínimo resíduo para resolver o sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ , com

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

(1.177)

e $\boldsymbol{b}=(3,3)$ . Usando uma aproximação inicial $\boldsymbol{x}^{(0)}=(0,0.5)$ , faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno da norma $\|\boldsymbol{r}(\boldsymbol{x})\|_{2}=\|\boldsymbol{b}-A\boldsymbol{x}\|_{2}$ e mostre as iterações do método sobre o gráfico.

E. 1.5.3.

Aplique a iteração da descida mais íngrime resíduo para resolver o sistema linear $A\boldsymbol{x}=\boldsymbol{b}$ , com

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

(1.178)

e $\boldsymbol{b}=(3,4.1)$ . Usando uma aproximação inicial $\boldsymbol{x}^{(0)}=(0.5,0)$ , faça uma análise geométrica das iterações do método.

Dica: faça um gráfico de contorno da norma $\|\boldsymbol{r}(\boldsymbol{x})\|_{2}=\|\boldsymbol{b}-A\boldsymbol{x}\|_{2}$ e mostre as iterações do método sobre o gráfico.

E. 1.5.4.

Considere o problema de difusão-advecção 2D do Exemplo 1.5.1, discretizado com o esquema upwind. Aplique a iteração da descida mais íngrime do resíduo para resolver o sistema linear resultante para os seguintes coeficientes de advecção:

a)

$\boldsymbol{a}=(-1,1)$ ;
b)

$\boldsymbol{a}=(1,-1)$ ;
c)

$\boldsymbol{a}=(-1,-1)$ .

Analise a convergência do método para diferentes tamanho de malha.

Dica: lembre-se que o esquema upwind deve ser adaptado para cada caso.

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

Política de uso de dados

Matemática Numérica III

1.5 Métodos de projeção

1.5.1 Métodos de projeção unidimensionais

1.5.2 Iteração da descida mais íngrime

1.5.3 Iteração do mínimo resíduo

Exemplo 1.5.1.(Problema de difusão-advecção 2D)

1.5.4 Iteração da descida mais íngrime do resíduo

1.5.5 Exercícios

E. 1.5.1.

E. 1.5.2.

E. 1.5.3.

E. 1.5.4.

Envie seu comentário