Matemática Numérica III

2 Sistemas Não Lineares 2 Sistemas Não Lineares 3 Sistemas Não Lineares e 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 $\boldsymbol{x}^{*}\in\mathbb{R}^{n}$ tal que

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}F(\boldsymbol% {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}(\boldsymbol{x})=[\jmath_{i,j}(\boldsymbol{x})]_{i,j=1}^{n,n}$ a matriz jacobiana³⁴³⁴endnote: ³⁴Carl 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}(% \boldsymbol{x})=\frac{\partial f_{i}(\boldsymbol{x})}{\partial x_{j}},

(2.2)

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

A iteração básica do método de Newton³⁵³⁵endnote: ³⁵Isaac 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 $\boldsymbol{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(\boldsymbol{x}^{(k)}\right)\boldsymbol{\delta}^{(k)}=-F% \left(\boldsymbol{x}^{(k)}\right)$		(2.4)
	computar:		(2.5)
	$\displaystyle\text{}\quad\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}% {rgb}{0,0,1}\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\boldsymbol{\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 $\boldsymbol{x}^{(0)}$ suficientemente próximo da solução $\boldsymbol{x}^{*}$ , o método de Newton é quadraticamente convergente. Mais precisamente, este resultado de convergência local requer que $J_{F}^{-1}(\boldsymbol{x}^{*})$ seja não singular e que $J_{F}$ seja Lipschitz³⁶³⁶endnote: ³⁶Rudolf Otto Sigismund Lipschitz, 1832 - 1903, matemático alemão. Fonte: Wikipédia. contínua. Consulte [8, 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³⁷³⁷endnote: ³⁷Jan 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³⁸³⁸endnote: ³⁸Johann 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³⁹³⁹endnote: ³⁹Erich Hans Rothe, 1895 - 1988, matemático alemão. Fonte: Wikipédia. com aproximação de Euler⁴⁰⁴⁰endnote: ⁴⁰Leonhard Paul Euler, 1707-1783, matemático e físico suíço. Fonte: Wikipédia: Ronald Fisher. 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 $\boldsymbol{w}\in\mathbb{R}^{n_{x}}$ , tal que

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

(2.15)

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

	$\displaystyle f_{1}\left(\boldsymbol{w};\boldsymbol{w}^{(0)}\right)=w_{1},$		(2.16)
	$\displaystyle f_{i}\left(\boldsymbol{w};\boldsymbol{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(\boldsymbol{w};\boldsymbol{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 [12]

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 [1]

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$ [2]:

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⁴¹⁴¹endnote: ⁴¹Isaac 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⁴²⁴²endnote: ⁴²Carl 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 Matriz 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 3.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 3.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 Matriz Jacobiana

E. 2.1.4.

Envie seu comentário