Matemática numérica III

2 Sistemas não-lineares 2 Sistemas não-lineares 3 Otimização

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

2.1 Método de Newton

Consideramos o seguinte problema: dada a função $F:D\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ encontrar $\bm{x}^{*}\in\mathbb{R}^{n}$ tal que

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}F(\bm{x}^{*})% =0.

(2.1)

Salvo explicitado ao contrário, assumiremos que $F\in C^{1}(D)$ , i.e. $F$ é uma função continuamente diferenciável no domínio computacional $D\subset\mathbb{R}^{n}$ .

Vamos, denotar por $J_{F}(\bm{x})=[\jmath_{i,j}(\bm{x})]_{i,j=1}^{n,n}$ a matriz jacobiana¹¹1Carl Gustav Jakob Jacobi, 1804 - 1851, matemático alemão. Fonte: Wikipédia: Carl Gustav Jakob Jacobi. da F com

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\jmath_{i,j}(% \bm{x})=\frac{\partial f_{i}(\bm{x})}{\partial x_{j}},

(2.2)

onde $F(\bm{x})=(f_{1}(\bm{x}),f_{2}(\bm{x}),\dotsc,f_{n}(\bm{x}))$ e $\bm{x}=(x_{1},x_{2},\dotsc,x_{n})$ .

A iteração básica do método de Newton²²2Isaac Newton, 1642 - 1727, matemático, físico, astrônomo, teólogo e autor inglês. Fonte: Wikipédia: Isaac Newton. para sistemas de equações consiste em: dada uma aproximação inicial $\bm{x}^{(0)}\in\mathbb{R}^{n}$ ,

	resolver:		(2.3)
	$\displaystyle\text{}\quad\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}% {rgb}{0,0,1}J_{F}\left(\bm{x}^{(k)}\right)\bm{\delta}^{(k)}=-F\left(\bm{x}^{(k% )}\right)$		(2.4)
	computar:		(2.5)
	$\displaystyle\text{}\quad\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}% {rgb}{0,0,1}\bm{x}^{(k+1)}=\bm{x}^{(k)}+\bm{\delta}^{(k)}$		(2.6)

para $k=0,1,2,\ldots$ até que um critério de parada seja satisfeito.

Observação 2.1.1.(Convergência)

Para $\bm{x}^{(0)}$ suficientemente próximo da solução $\bm{x}^{*}$ , o método de Newton é quadraticamente convergente. Mais precisamente, este resultado de convergência local requer que $J_{F}^{-1}(\bm{x}^{*})$ seja não singular e que $J_{F}$ seja Lipschitz³³3Rudolf Otto Sigismund Lipschitz, 1832 - 1903, matemático alemão. Fonte: Wikipédia: Rudolf Lipschitz. contínua. Consulte [10, Seção 7.1] para mais detalhes.

Código 12: newton.py

⬇

1from numpy.linalg import solve, norm

3def newton(f, jac, x0, tol=1e-10, max_iter=100):

4 x = x0.copy()

5 info = 0

6 for i in range(max_iter):

7 # jacobiana

8 J = jac(x)

9 # resíduo

10 F = f(x)

11 # método linear

12 delta = solve(J, -F)

13 # atualização

14 x = x + delta

15 # critério de parada

16 if norm(delta) < tol:

17 info = 1

18 return x, info

20 return x, info

Exemplo 2.1.1.

Consideremos a equação de Burgers⁴⁴4Jan Burgers, 1895 - 1981, físico neerlandês. Fonte: Wikipédia: Jan Burgers.

\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu\frac{\partial% ^{2}u}{\partial x^{2}}

(2.7)

com $\nu>0$ , condição inicial

u(0,x)=-\operatorname{sen}(\pi x)

(2.8)

e condições de contorno de Dirichlet⁵⁵5Johann Peter Gustav Lejeune Dirichlet, 1805 - 1859, matemático alemão. Fonte: Wikipédia: Johann Peter Gustav Lejeune Dirichlet. homogêneas

u(t,-1)=u(t,1)=0.

(2.9)

Aplicando o método de Rothe⁶⁶6Erich Hans Rothe, 1895 - 1988, matemático alemão. Fonte: Wikipédia: Erich Rothe. com aproximação de Euler⁷⁷7Leonhard Paul Euler, 1707-1783, matemático e físico suíço. Fonte: Wikipédia: Leonhard Euler. implícita, obtemos

	$\displaystyle\frac{u(t+h_{t},x)-u(t,x)}{h_{t}}+$
	$\displaystyle\text{}\quad u(t+h_{t},x)u_{x}(t+h_{t},x)\approx\nu u_{xx}(t+h_{t% },x),$		(2.10)

onde $h_{t}>0$ é o passo no tempo. Agora, aplicamos diferenças finitas para obter

	$\displaystyle\frac{u(t+h_{t},x_{i})-u(t,x_{i})}{h_{t}}$
	$\displaystyle\text{}\quad+u(t+h_{t},x_{i})\frac{u(t+h_{t},x_{i+1})-u(t+h_{t},x% _{i-1})}{2h_{x}}$
	$\displaystyle\text{}\quad\approx\nu\frac{u(t+h_{t},x_{i-1})-2u(t+h_{t},x_{i})+% u(t+h_{t},x_{i+1})}{h_{x}^{2}},$		(2.11)

onde, $x_{i}=(i-1)h_{x}$ , $i=1,\dotsc,n_{x}$ e $h_{x}=1/(n_{x}-1)$ é o tamanho da malha.

Rearranjando os termos e denotando $u^{(k)}_{i}\approx u(t_{k},x_{i})$ , $t_{k}=(k-1)h$ , obtemos o seguinte sistema de equações não-lineares, que consiste no problema discreto associado: para cada $k=0,1,2,\ldots$ , encontrar $u^{(k+1)}_{i}$ , $i=1,\dotsc,n_{x}$ , tais que

	$\displaystyle u^{(k+1)}_{1}=0$		(2.12)
	$\displaystyle\frac{1}{h_{t}}u^{(k+1)}_{i}-\frac{1}{h_{t}}u^{(k)}_{i}+\frac{1}{% 2h_{x}}u^{(k+1)}_{i}u^{(k+1)}_{i+1}$
	$\displaystyle\text{}\quad-\frac{1}{2h_{x}}u^{(k+1)}_{i}u^{(k+1)}_{i-1}-\frac{% \nu}{h_{x}^{2}}u^{(k+1)}_{i-1}$
	$\displaystyle\text{}\quad+2\frac{\nu}{h_{x}^{2}}u^{(k+1)}_{i}-\frac{\nu}{h_{x}% ^{2}}u^{(k+1)}_{i+1}=0,$		(2.13)
	$\displaystyle u^{(k+1)}_{n_{x}}=0,$		(2.14)

sendo $u^{(0)}_{i}=-\operatorname{sen}(\pi x_{i})$ , $i=1,\dotsc,n_{x}$ .

Este problema pode ser reescrito como segue: para cada $k=0,1,\ldots$ , encontrar $\bm{w}\in\mathbb{R}^{n_{x}}$ , tal que

F\left(\bm{w};\bm{w}^{(0)}\right)=0,

(2.15)

onde $w_{i}=u^{(k+1)}_{i}$ , $w_{i}^{(0)}=u^{(k)}_{i}$ e

	$\displaystyle f_{1}\left(\bm{w};\bm{w}^{(0)}\right)=w_{1},$		(2.16)
	$\displaystyle f_{i}\left(\bm{w};\bm{w}^{(0)}\right)=\frac{1}{h_{t}}w_{i}-\frac% {1}{h_{t}}w^{(0)}_{i}+\frac{1}{2h_{x}}w_{i}w_{i+1}$
	$\displaystyle\quad-\frac{1}{2h_{x}}w_{i}w_{i-1}-\frac{\nu}{h_{x}^{2}}w_{i-1}+2% \frac{\nu}{h_{x}^{2}}w_{i}-\frac{\nu}{h_{x}^{2}}w_{i+1},$		(2.17)
	$\displaystyle f_{n_{x}}\left(\bm{w};\bm{w}^{(0)}\right)=w_{n_{x}}.$		(2.18)

A matriz jacobiana associada $J=[\jmath_{i,j}]_{i,j}^{n_{x},n_{x}}$ contém

	$\displaystyle\jmath_{i,j}=0,\quad j\neq i-1,i,i+1,$		(2.19)
	$\displaystyle\jmath_{1,1}=1,$		(2.20)
	$\displaystyle\jmath_{1,2}=0,$		(2.21)
	$\displaystyle\jmath_{i,i-1}=\frac{1}{2h_{x}}w_{i}-\frac{\nu}{h_{x}^{2}},$		(2.22)
	$\displaystyle\jmath_{i,i}=\frac{1}{h_{t}}+\frac{1}{2h_{x}}w_{i+1}-\frac{2}{h_{% x}}w_{i-1}+\frac{2\nu}{h_{x}^{2}},$		(2.23)
	$\displaystyle\jmath_{i,i+1}=\frac{1}{2h_{x}}w_{i}-\frac{\nu}{h_{x}^{2}},$		(2.24)
	$\displaystyle\jmath_{n_{x},n_{x}-1}=0$		(2.25)
	$\displaystyle\jmath_{n_{x},n_{x}}=1.$		(2.26)

Código 13: burgers_newton.py

⬇

1from numpy import linspace, sin, pi, zeros

3# parâmetros

4nu = 0.01

6# malha espacial

7N = 100

8L = 1.0

9dx = L/(N+1)

10x = linspace(0, L, N+1)

12# malha temporal

13dt = 0.001

14T = 1.0

15nt = int(T/dt)

17# c.i.

18t0 = 0.0

19u0 = sin(pi*x)

21# track

22U = [u0.copy()]

23nt_track = 100

25def f_burgers(u):

26 F = zeros(N+1)

27 F[0] = u[0] # u(0,t) = 0

28 for i in range(1, N):

29 F[i] = (u[i] - u0[i])/dt + u[i]* (u[i+1] - u[i-1])/(2*dx) - nu*(u[i+1] - 2*u[i] + u[i-1])/(dx*dx)

30 F[N] = u[N] # u(1,t) = 0

31 return F

33def jac_burgers(u):

34 J = zeros((N+1, N+1))

35 J[0, 0] = 1.0 # u(0,t) = 0

36 for i in range(1, N):

37 J[i, i] = 1/dt + (u[i+1] - u[i-1])/(2*dx) + 2*nu/(dx*dx)

38 J[i, i+1] = u[i]/(2*dx) - nu/(dx*dx)

39 J[i, i-1] = -u[i]/(2*dx) - nu/(dx*dx)

40 J[N, N] = 1.0 # u(1,t) = 0

41 return J

43# solução temporal

44u = u0.copy()

45for n in range(nt):

46 u_new, info = newton(f_burgers, jac_burgers, u)

47 u0 = u.copy() # atualizar c.i.

48 u = u_new.copy() # atualizar solução

50 if (n+1) % nt_track == 0:

51 U.append(u.copy()) # armazenar solução

Exercícios

E. 2.1.1.

Faça uma implementação pelo método de diferenças finitas com método de Euler implícito para resolver o seguinte problema de Burgers:

	$\displaystyle u_{t}+uu_{x}=\nu u_{xx},\leavevmode\nobreak\ (t,x)\in(0,t_{f}]% \times(-1,1),$		(2.27)
	$\displaystyle u(0,x)=2\nu\pi\frac{\operatorname{sen}(\pi x)}{2+\cos(\pi x)},% \leavevmode\nobreak\ x\in[-1,1],$		(2.28)
	$\displaystyle u(t,-1)=u(t,1)=0,\leavevmode\nobreak\ t\in[0,t_{f}],$		(2.29)

Teste seu código para $\nu=1.,0.1,0.01,0.001,0.0001$ e $t_{f}=1.0$ . Dica: a solução analítica é dada por [Wood2006]

u(t,x)=2\nu\pi\frac{e^{-\nu\pi^{2}t}\operatorname{sen}(\pi x)}{2+e^{-\nu\pi^{2% }t}\cos(\pi)x}.

(2.30)

Dica: $\displaystyle u(t,x)=2\nu\pi\frac{e^{-\nu\pi^{2}t}\operatorname{sen}(\pi x)}{2% +e^{-\nu\pi^{2}t}\cos(\pi)x}$ .

E. 2.1.2.

Faça uma implementação pelo método de diferenças finitas com método de Euler implícito para resolver o seguinte problema de Burgers:

	$\displaystyle u_{t}+uu_{x}=\nu u_{xx},\leavevmode\nobreak\ (t,x)\in(0,t_{f}]% \times(-1,1),$		(2.31)
	$\displaystyle u(0,x)=-\operatorname{sen}(\pi x),\leavevmode\nobreak\ x\in[-1,1],$		(2.32)
	$\displaystyle u(t,-1)=u(t,1)=0,\leavevmode\nobreak\ t\in[0,t_{f}],$		(2.33)

Teste seu código para $\nu=1.,0.1,0.01,0.001,0.0001$ e $t_{f}=1.0$ . Dica> a solução analítica é dada por [Basdevant1986]

u(t,x)=\frac{-\int_{-\infty}^{\infty}\operatorname{sen}(\pi(x-\eta))f(x-\eta)e% ^{-\frac{\eta^{2}}{4\nu t}}d\eta}{\int_{-\infty}^{\infty}f(x-\eta)e^{-\frac{% \eta^{2}}{4\nu t}}d\eta},

(2.34)

onde $f(y)=\exp\left(\frac{\cos(\pi y)}{2\pi\nu}\right)$ e $\nu=0.01/\pi$ .

E. 2.1.3.

Faça uma implementação pelo método de diferenças finitas com método de Euler implícito para resolver o seguinte problema de Burgers:

	$\displaystyle u_{t}+uu_{x}=\nu u_{xx},\leavevmode\nobreak\ (t,x)\in(0,t_{f}]% \times(-1,1),$		(2.35)
	$\displaystyle u(0,x)=2\operatorname{sign}(x),\leavevmode\nobreak\ x\in[-1,1],$		(2.36)
	$\displaystyle u(t,-1)=u_{a},u(t,1)=u_{b},\leavevmode\nobreak\ t\in[0,t_{f}],$		(2.37)

Teste seu código para $\nu=1.,0.1,0.01,0.001,0.0001$ e $t_{f}=1.0$ . Dica: é conhecida a seguinte solução analítica para $\nu=1$ [Benton1972]:

u(t,x)=2\frac{G(t,x)-G(t,-x)}{G(t,x)+G(t,-x)},

(2.38)

onde

G(t,x)=\frac{1}{2}e^{t-x}\operatorname{erfc}\frac{2t-x}{2\sqrt{t}}.

(2.39)

2.1.1 Método tipo Newton

Em revisão

Existem várias modificações do Método de Newton⁸⁸8Isaac Newton, 1642 - 1727, matemático, físico, astrônomo, teólogo e autor inglês. Fonte: Wikipédia: Isaac Newton. que buscam reduzir o custo computacional. Há estratégias para simplificar as computações da matriz jacobiana⁹⁹9Carl Gustav Jakob Jacobi, 1804 - 1851, matemático alemão. Fonte: Wikipédia: Carl Gustav Jakob Jacobi. e para reduzir o custo nas computações das atualizações de Newton.

Atualização cíclica da jacobiana

Em revisão

Geralmente, ao simplificarmos a matriz jacobina $J_{F}$ ou aproximarmos a atualização de Newton $\delta^{(k)}$ , perdemos a convergência quadrática do método (consulte a Observação 2.1.1). A ideia é, então, buscar uma convergência pelo menos super-linear, i.e.

\|e^{(k+1)}\|\approx\rho_{k}\|e^{(k)}\|,

(2.40)

com $\rho_{k}\to 0$ quando $k\to\infty$ . Aqui, $e^{(k)}:=x^{*}-x^{(k)}$ . Se a convergência é superlinear, então podemos usar a seguinte aproximação

\rho_{k}\approx\frac{\|\delta^{(k)}\|}{\|\delta^{(k-1)}\|}

(2.41)

ou, equivalentemente,

\rho_{k}\approx\left(\frac{\|\delta^{(k)}\|}{\|\delta^{(0)}\|}\right)^{\frac{1% }{k}}

(2.42)

Isso mostra que podemos acompanhar a convergência das iterações pelo fator

\beta_{k}=\frac{\|\delta^{(k)}\|}{\|\delta^{(0)}\|}.

(2.43)

Ao garantirmos $0<\beta_{k}<1$ , deveremos ter uma convergência superlinear.

Vamos, então, propor o seguinte método tipo Newton de atualização cíclica da matriz jacobiana.

1.

Escolha $0<\beta<1$
2.

$J:=J_{F}(x^{(0)})$
3.

$J\delta^{(0)}=-F(x^{(0)})$
4.

$x^{(1)}=x^{(0)}+\delta^{(0)}$
5.
Para $k=1,\ldots$ até critério de convergência:
1. (a)
  
  $J\delta^{(k)}=-F(x^{(k)})$
2. (b)
  
  $x^{(k+1)}=x^{(k)}+\delta^{(k)}$
3. (c)
  Se $\|\delta^{(k)}\|/\|\delta^{(0)}\|>\beta$ :
  1. i.
    
    $J:=J_{F}(x^{(k)})$

E. 2.1.4.

Implemente uma versão do método tipo Newton apresentado acima e aplique-o para simular o problema discutido no Exemplo 2.1.1 para $\nu=1.,0.1,0.01,0.001,0.0001$ . Faça uma implementação com suporte para matrizes esparsas.

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

2.1 Método de Newton

Observação 2.1.1.(Convergência)

Exemplo 2.1.1.

Exercícios

E. 2.1.1.

E. 2.1.2.

E. 2.1.3.

2.1.1 Método tipo Newton

Atualização cíclica da jacobiana

E. 2.1.4.

Envie seu comentário