Fractais e o método de Newton

Apesar do método de Newton datar do século XVII, o teorema de convergência do método de Newton, em sua forma mais geral, foi demonstrado apenas em 1948, por Leonid Kantorovich . No vídeo, apresentei o teorema mais simples, cuja demonstração está bem feita em Dennis e Schnabel (1996).

Teorema: Sejam $f$ uma função com segunda derivada contínua $(f\in C^2)$ e $x_*$ um zero de $f$ tal que $f'(x_*) \neq 0.$ Então existe uma vizinhança $V$ de $x_*$ tal que para todo $x_0 \in V$ a sequência gerada pelo método de Newton converge quadraticamente para $x_*,$ isto é, \begin{equation} \lim_{k\to\infty} \frac{|x_{k+1}-x_*|}{|x_k-x_*|^2} = C_*, \label{eqn:convrate} \end{equation} onde $C_* = {1\over 2}|f''(x_*) / f'(x_*)|.$

Duas observações importantes devemos fazer a respeito deste teorema. A primeira é que \eqref{eqn:convrate} significa que para qualquer $C > C_*,$ existe $K$ tal que para todo $k\ge K,$ $$|x_{k+1}-x_*| \le C |x_k-x_*|^2, $$ ou seja, a partir de certo ponto, o erro reduzirá a um fator do erro da iteração passada ao quadrado. Essa taxa de convergência quadrática é muito boa, porém a sutileza do resultado é que esta taxa é garantida no limite, ou seja, pode ser necessário que muitas iterações sejam computadas até que a convergência quadrática seja observada.

Como exemplo, considere $f(x) = -{1\over 2}x^2+3x-4.$ Os zeros de $f$ são $2$ e $4.$ Veja que $f'(x) = -x+3$ e $f''(x) = -1.$ Sendo assim, para $x_* =2,$ temos $C_* = {1\over 2}|f''(x_*) / f'(x_*)| = {1\over 2}.$ Tome $x_0 = 0$ (será que $0$ está na tal vizinhança de $x_*$ que o teorema exige?). A sequência gerada pelo método de Newton é $$ \begin{array}{ll} x_0 = 0, & |x_0-x_*| = 2 \\ x_1 = 1.3333, & |x_1-x_*| = 6.6667\cdot 10^{-1}\\ x_2 = 1.8667, & |x_2-x_*| = 1.3333\cdot 10^{-1}\\ x_3 = 1.9922, & |x_3-x_*| = 7.8431\cdot 10^{-3}\\ x_4 = 1.99996948, & |x_4-x_*| = 3.0518\cdot 10^{-5} \approx {1\over 2}|x_3-x_*|^2 \\ \end{array} $$

Outra observação sobre o teorema é que está assegurada a existência dessa vizinhança $V$ de $x_*,$ dentro da qual há garantia de convergência da sequência gerada pelo método de Newton. Porém, o teorema não dá qualquer indicação de como determinar essa vizinhança. Isto é, sabemos que ela existe, só não sabemos qual é.

Após a metade do século XX, com o advento dos computadores, pesquisadores começaram a identificar numericamente (e visualmente) essas regiões de convergência. O experimento era escolher diversos pontos iniciais para o método de Newton e observar para qual zero da função o método convergia. O ponto inicial era então marcado com a cor que identificava o zero atingido. O teorema nos garante o que acontece quando $x_0$ está próximo o suficiente de um zero. Mas ninguém imaginava o que aconteceria quando $x_0$ estava fora dessas vizinhanças preditas pelo teorema. Os experimentos produziam figuras mais interessantes quando realizados no plano, por exemplo considerando funções de variáveis complexas. Foi assim que os primeiros fractais relacionados ao método de Newton foram construídos.

Para ir além, recomendo o livro Greenbaum e Chartier (2012).

Fractal no Octave

Podemos fazer um experimento no Octave para tentar produzir um fractal como exibido no vídeo. Naquele exemplo, estávamos interessados em estimar os zeros de $f(z)=z^3-1.$ Sabemos que esse polinômio, definido nos complexos, tem 3 zeros, sendo eles $1,$ $(\sqrt{3} + i)/2,$ e $(\sqrt{3} -i)/2.$ Vamos tentar descobrir para qual deles o método de Newton convergirá quando iniciado em diferentes de pontos do plano.

O script a seguir inicia definindo a função $f$, sua derivada, e uma malha sobre a região $[-1.5,1.5]\times[-1.5,1.5].$ Depois, utiliza cada ponto desta malha como aproximação inicial para o método de Newton aplicado à função $f$. Para rodar esse script, você precisa ter a função newton.m, que foi trabalhada na aula passada.

f = @(z) z^3-1;   # definição da função f: ℂ → ℂ
g = @(z) 3*z^2;   # derivada da função g(z) = f'(z)

# Vetor com os zero de f
sol = [1; (-1 + sqrt(3)*i)/2; (-1 - sqrt(3)*i)/2];

# Quantidade de pontos em cada direção da malha
n = 200;

# Região do plano onde pontos iniciais serão tomados
# (essa região deve conter os zeros de f)
x = linspace(-1.5, 1.5, n);
y = linspace(-1.5, 1.5, n);

# Matriz para guardar os resultados
C = zeros(n,n);

# Tolerância para o critério de parada do Newton
tol = 1.0e-6;

# Loop sobre todos os pontos da malha
for ix = 1:n
   for iy = 1:n
      z0 = x(ix) + y(iy) * i;           # ponto inicial
      [z,fz] = newton(f,g,z0,tol,20);

      if abs(fz) <= tol                 # testa se houve convergência
	 [err, I] = min(abs(z - sol));  # qual zero está mais próximo de z?
	 C(iy,ix) = I;                  # guarda o índice do zero mais próximo
      else
	 C(iy,ix) = NaN;                # NaN, no caso de não convergência
      endif
   endfor
endfor

imagesc(x,y,C)
axis("xy");
colorbar;

Experimente o código acima e veja a figura produzida. Depois brinque com ele. Por exemplo, altere o valor n para gerar uma figura com mais detalhes. Experimente outras funções $f$ (só não esqueça de trocar também a definição da derivada e o vetor que contém os zeros de $f$).

Referência

J.E. Dennis, Jr. e R.B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM, 1996.

Anne Greenbaum e Timothy P. Chartier. Numerical Methods: Design, Analysis, and Computer Implementation of Algorithm. Princeton University Press, 2012.

1. Construa a iteração de Newton para estimar o zero da função $f(x) = x^3.$ Partindo de $x_0 = 1,$ observe como o erro reduz a cada iteração. A sequência aparenta estar convergindo quadraticamente? Por quê?

2. A convergência do método de Newton só pode ser garantida numa vizinhança da solução. Grafique as funções abaixo e observe o comportamento das iterações produzidas pelo método de Newton a partir do ponto inicial $x_0$ dado. Para cada função, exiba outro ponto inicial, para o qual a iteração de Newton de fato converge para o zero da função.

1. $f(x) = xe^{-x},$ $x_0 = 2$
2. $g(x) = x^3 - x - 3,$ $x_0 = 0$
3. $h(x) = \arctan(x),$ $x_0 = 1.45$

3. Seja $f(x) = (1+e^x)^{-1} - {1\over 2}.$ O único zero de $f$ está em $x=0.$

1. Faça o gráfico de $f$ e tente "adivinhar" o que acontece quando o método de Newton é usado para estimar esse zero, para diferentes escolhas de $x_0.$
2. Compute a sequência produzida pelo método de Newton, partindo de $x_0 = 1.$
3. Compute a sequência produzida pelo método de Newton, partindo de $x_0 = 3.$
4. Tente empiricamente descobrir o valor de $x_0$ que está no limite entre os dois casos observados acima.
5. Mostre que para $x_0,$ solução de $2x_0 = \sinh(x_0)$, a sequência gerada pelo método de Newton é $x_0,-x_0,x_0,\ldots$

4. Seja $f: \mathbb{R} \to \mathbb{R}.$ Se $x_∗$ for um zero de $f,$ defina $I(x^∗)$ como o maior intervalo contendo $x^∗$ tal que para todo $x_0 \in I(x^∗),$ a sequência gerada pelo método de Newton, partindo de $x_0$ , converge para $x^∗.$ Considere $f(x) = x^2 \sin x + x \cos(2x) − 3.$

1. Aplique o método de Newton para encontrar os primeiros dois zeros positivos de $f,$ digamos $x^∗_1$ e $x^∗_2$ , e os primeiros dois zeros negativos de $f,$ digamos $x^∗_{−1}$ e $x^∗_{−2}.$
2. Estime numericamente $I(x^∗_{−1})$ e $I(x^∗_1).$