Algoritmo para método de Newton

O núcleo do algoritmo para o método de Newton é bem simples. Ele consiste apenas de computar sucessivamente uma aproximação linear para a função e resolver o problema de encontrar o zero dessa aproximação, tomando-o como nova estimativa para o zero da função.

Entretanto, ao sair do campo teórico para o computacional, devemos nos atentar a algumas questões importantes. Vimos que a convergência do método de Newton é assegurada apenas quando partimos de um ponto inicial próximo do zero da função. E se este não for o caso?

Quando dissemos que as aproximações são computadas sucessivamente, precisamos especificar até quando ou por quanto tempo, visto que na prática não é possível iterar ao infinito.

Por fim, como a fórmula de iteração do Método de Newton realiza uma divisão, pode acontecer do denominador, $f'(x_k),$ ser zero. O que fazer então nessa situação?

Vamos começar lidando com as considerações sobre o ponto inicial. Como o método de Newton precisa de uma boa estimativa inicial, uma estratégia pode ser utilizar um método menos exigente, que funcione mesmo que essa aproximação inicial não seja tão boa. Depois desse método melhorá-la um pouco, então emprega-se o método de Newton. O método da bissecção poderia ser utilizado com essa finalidade. Claro que aí, recairíamos nas dificuldades do método da bissecção.

Outra estratégia é globalizar o método de Newton, ou seja, modificá-lo para que o método se torne robusto em relação ao ponto de partida. As duas principais formas de globalizar o método de Newton são através da inserção de uma busca linear na direção apontada pelo método de Newton ou através da definição de uma região de confiança, dentro da qual a aproximação linear é válida. Ambas as estratégias são melhor discutidas no contexto de métodos de otimização e fogem ao escopo desta aula. Para quem tem interesse, recomendo o livro de Dennis e Schnabel (1996).

Passando à segunda questão, um algoritmo computacional não pode deixar de explicitar o critério de parada de execução. O desejável seria pedir que o erro ficasse abaixo de um certo nı́vel aceitável, ou seja, dado um $\epsilon \gt 0,$ gostaríamos de iterar até que $$ |e_k| \le \epsilon, $$ onde $e_k = (x_k-x_*)$ é o erro no passo $k.$ Porém, o problema é como medir o erro absoluto, sem o conhecimento de $x_*.$

Um dos critérios que sim podemos aplicar é sobre o valor da $f(x_k).$ A lógica é: se $f(x_*)=0$ então $|f(x_k)|$ pequeno deve significar que $x_k$ está próximo de $x_*.$ Será mesmo?

Quando a sequência ${x_k}$ converge para $x_*,$ esta relação de fato pode ser estabelecida. Porém, a constante de proporcionalidade está relacionada ao valor de $f'(x_∗).$ Para ver isso, considere o polinômio de Taylor de segunda ordem para $f$ em torno de $x_*,$ supondo que $f$ tem até segunda derivada contínua. Observe que $$ f(x_k) = f(x_*) + f'(x_*)(x_k-x_*) + {f''(\xi)\over 2}(x_k-x_*)^2. $$ Se $x_k$ estiver convergindo para $x_*,$ e como $f(x_*) = 0$, temos que $$ f(x_k) \approx f'(x_*)(x_k-x_*). $$ Logo $$ |e_k| = |x_k-x_*| \approx \left|{f(x_k)\over f'(x_*)}\right| \lessapprox {\epsilon \over |f'(x_*)|}, $$ se $|f(x_k)| \le \epsilon.$

Desta forma, usar como critério de parada $|f(x_k)|\le \epsilon$ pode ser muito restritivo, se $|f'(x_*)|$ for grande, ou muito permissivo, se $|f'(x_*)|$ for pequeno demais. Quando o valor de $f'(x_*)\approx 1,$ teríamos a melhor situação, no sentido de que o valor de $|f(x_k)|$ seria um bom parâmetro para estimar o valor de $|e_k|.$ Isto é a essência do que visualmente está representado na figura a seguir.

Na prática, os critérios mais usados, além de um limite máximo para o número de iterações, são:

• $|f(x_k)|\le \epsilon |f(x_0)|,$ ou seja, uma redução relativa no valor de função. Assim o parâmetro $\epsilon$ pode ser um número adimensional. O problema deste critério é que se $|f(x_0)|$ for muito grande, o método pode parar ainda com $|f(x_k)|$ muito grande. Por outro lado, se $|f(x_0)|$ já for muito pequeno, então o método vai se esforçar exageradamente para reduzir demais o valor de $|f(x_k)|.$
• $|f(x_k)| \leq \epsilon_1 |f(x_0)| + \epsilon_2,$ ou seja, uma combinação entre uma tolerância relativa e uma absoluta para o valor de função. A escolha apropriada de $\epsilon_1$ e $\epsilon_2$ evita as mazelas do critério anterior. Por exemplo, tente entender o que acontece quando $\epsilon_1 = 10^{-6}$ e $\epsilon_2 = 10^{-6},$ com relação às diferentes possibilidades para $|f(x_0)|.$
• $|s_k| \leq \epsilon |s_0|,$ para $s_k = x_{k+1} - x_k$, ou seja, uma redução relativa no comprimento dos passos dados pelo método. Isso evita insistir nas iterações quando o método parece estagnar.

Ainda sobre o critério de parada, um último comentário é que, seja qual for sua escolha, sempre é importante limitar também o número de máximo de iterações. Esta é a única garantia de que o algoritmo de fato vai parar em algum momento.

Por fim, vamos discutir o último aspecto levantado no início desta aula, o caso em que $f'(x_k) = 0.$ A iteração de Newton para ser computada precisa de uma divisão por $f'(x_k).$ Se essa derivada for nula, a única possibilidade é evitar isso, quer seja pela escolha de outro ponto, quer seja pela troca momentânea do método de iteração.

Claro que a chance de você acertar um ponto, no curso das iterações de Newton, que tenha a derivada exatamente nula, é muitíssima pequena, nula até. Será que podemos nos tranquilizar então?

Na prática, a derivada ser quase nula já é um problema. Como a reta tangente num ponto onde a derivada é muito pequena é praticamente paralela ao eixo horizontal, sua intersecção com o eixo (próximo iterando do método de Newton) estará muito distante. Com isto, o iterando pode se afastar da região onde procurávamos por um zero da função, levando à convergência para outro zero ou mesmo à divergência. Veja o exemplo da figura abaixo, onde a função tem o eixo $x$ como uma assíntota. Apesar do zero da função estar próximo de $x_0,$ a sequência gerada pelo método de Newton se afastou desse zero, inclusive divergindo.

Para uma discussão mais aprofundada sobre a implementação do método de Newton, recomendo o livro Kelley (2003).

Referências

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

C. T. Kelley. Solving Nonlinear Equations with Newton's Method. SIAM, 2003.

Considere o seguinte protótipo para uma implementação do método de Newton em Octave.

  function [x, fx, n] = newton(f, g, x0, tol, N)
   x = x0;
   fx = f(x);
   n = 0;
   s = tol + 1;

   while (s > tol & n <= N)
      gx = g(x);
      s = - fx / gx;
      x = x + s;
      s = abs(s);
      fx = f(x);
      n = n + 1;
   endwhile
endfunction

1. No modelo de implementação acima, identifique qual o papel de cada variável da função.

2. Utilize essa função no Octave, para aproximar os zeros das funções abaixo. Inicialmente, utilize como ponto inicial o $x_0$ fornecido. Faço o gráfico das funções e observe o comportamento da sequência gerada pelo método de Newton. Experimente outros valores para o ponto inicial.

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

3. Sobre a implementação acima, qual o critério de parada que está sendo utilizado? Como a implementação deveria ser modificada para utilizar algum dos outros critérios apresentados? Quais pontos da discussão dessa aula não foram contemplados nessa implementação. Tente melhorar a função apresentada, tornando-a mais robustas, ou seja, menos suscetível a erros ou fracassos.

4. Baixe a rotina newtsol.m, distribuída juntamente com o livro Kelley (2003). Utilize-a nos mesmos testes do exercício 2. Estude como a essa rotina foi implementada. Tente identificar os pontos importantes da implementação.