6. Método de Newton: Sair pela Tangente sem Perder o Rumo#

Objetivos de aprendizagem

  • Compreender o funcionamento do método de Newton, seu algoritmo e as condições necessárias para sua aplicação;
  • Aplicar o método de Newton na resolução de equações não lineares;
  • Analisar o comportamento da convergência do método de Newton, sua precisão e limitações;
  • Reconhecer variantes do Método de Newton para problemas de otimização, bem como o Método de Halley.

"Nenhuma grande descoberta jamais foi feita sem um palpite ousado." (Isaac Newton)

O método de Newton, também conhecido como método de Newton-Raphson, é uma técnica iterativa poderosa para encontrar aproximações das raízes de equações não lineares. Sua ideia central é usar a reta tangente à curva da função em um ponto de aproximação inicial para estimar onde essa reta cruza o eixo \(x\), fornecendo assim uma nova aproximação. Esse processo é repetido até que a diferença entre iterações sucessivas seja suficientemente pequena, indicando convergência para uma raiz.

Apesar de ser um método rápido e eficiente com convergência quadrática quando a aproximação inicial está suficientemente próxima da raiz, ele também apresenta desafios. Ele exige que a derivada \(f’(x)\) esteja disponível e bem comportada. Além disso, se o ponto inicial estiver longe da raiz, ou se a função tiver pontos de inflexão ou derivadas próximas de zero, o método pode divergir ou oscilar. Por isso, é comum combiná-lo com outras técnicas, como a inspeção gráfica ou métodos mais estáveis (como bisseção) para estimar bons pontos iniciais.

../_images/newton-ai.png

6.1. Explicação do algoritmo#

O processo iterativo do método é dado pela seguinte função de iteração:

\[ x^{(k)} = \phi(x^{(k-1)}) = x^{(k-1)} - \dfrac{ f(x^{(k-1)} )}{ f'(x^{(k-1)}) }, \ \ x^{(0)} = x_0, \ \ k = 1,2, \ldots < \infty, \ \ f'(x^{(k-1)}) \neq 0, \ \ \forall k. \]

onde:

  • \(x^{(k-1)}\) é a aproximação atual.

  • \(f(x^{(k-1)})\) é o valor atual da função.

  • \(f'(x^{(k-1)})\) é o valor atual da primeira derivada da função.

A ideia é que a interseção da tangente à curva \(f(x)\) no ponto \(x^{(k-1)}\) com o eixo \(x\) forneça uma nova aproximação \(x^{(k)}\). O Método de Newton é uma ferramenta essencial no arsenal de métodos numéricos, oferecendo uma abordagem eficiente e poderosa para a solução de uma ampla variedade de problemas matemáticos.

O algoritmo funciona da seguinte forma:

  1. Inicialize \(x^{(0)}\) (chute inicial)

  2. Para \(k\) de 1 até \(N\):

    a. Calcule \(f(x^{(k)})\)

    b. Calcule \(f'(x^{(k)})\) (derivada da função no ponto \(x^{(k)}\))

    c. Se \(f'(x^{(k)}) = 0\), retorne erro (derivada zero, o método falha)

    d. Calcule \(x^{(k)} = \phi(x^{(k-1)})\) (como o processo iterativo acima)

    e. Se \(ER(x^{(k)}, x^{(k-1)}) < \epsilon\), então a solução foi encontrada com a precisão desejada. Retorne \(x^{(k)}\).

  3. Se o número máximo de iterações for atingido, retorne a última aproximação \(x^{(k-1)}\).

6.2. Implementação do algoritmo#

No algoritmo proposto abaixo, a função newton requererá 5 argumentos:

  • a estimativa inicial, representada pela variável x0;

  • a função \(f(x)\) propriamente dita, representada por f;

  • a derivada \(f'(x)\), representada por df;

  • o erro relativo assumido, representado por tol;

  • o número máximo de iterações \(N\) para tentativa de solução, representado por nmax.

import inspect, re, numpy as np
from matplotlib.pyplot import plot
from prettytable import PrettyTable as pt

def newton(x0,f,df,tol,N):
    """Algoritmo para determinação de raízes pelo método de Newton.

    Parâmetros: 
        x0: estimativa inicial
        f: string dependendo de uma variável, i.e., a função-alvo
            (e.g., 'x**2 - 1', 'x**2*cos(x)', etc.) 
        df: string dependendo de uma variável, i.e., a derivada da função-alvo
        tol: erro desejado (tolerância)
        N: número máximo de iterações a repetir

    Retorno: 
        x: aproximação para a raiz da função    
    """

    # construtor de tabela
    table = pt()
    
    # substitui expressões da string pelas chamadas das funções do numpy
    f = re.sub('(sin|sinh|cos|cosh|tan|tanh|exp|log|sqrt|log10|arcsin|arccos|arctan|arcsinh|arccosh|arctanh)', r'np.\1', f)
    df = re.sub('(sin|sinh|cos|cosh|tan|tanh|exp|log|sqrt|log10|arcsin|arccos|arctan|arcsinh|arccosh|arctanh)', r'np.\1', df)
    
    # identifica a variável independente em f
    var = re.search(r'([a-zA-Z][\w]*) ?([\+\-\/*]|$|\))', f).group(1)

    
    # cria função
    f = eval('lambda ' + var + ' :' + f)
    
    # checa se a função é de uma variável, senão lança erro        
    try: 
        len(inspect.getfullargspec(f).args) - 1 > 0
    except:
        raise ValueError('O código é válido apenas para uma variável.')
    finally:
        # cria função derivada
        df = eval('lambda ' + var + ' :' + df)
    
    it = 0 # contador de iteracoes
    
    # cabeçalho de tabela
    table.field_names = ['i','x','f(x)','f\'(x)','ER']

    # imprime estimativa inicial
    print(f'Estimativa inicial: x0 = {x0:.6f}\n')  

    # Loop 
    for i in range(0,N):
        
        x = x0 - f(x0)/df(x0) # funcao de iteracao 
        
        e = abs(x-x0)/abs(x) # erro
        
        # tabela
        # impressão de tabela
        table.add_row([i,np.round(x,8),np.round(f(x),8),np.round(df(x),4),f'{e:e}'])
        table.align = 'c';      
        
        if e < tol:
            break
        x0 = x                
        
    print(table)
       
    if i == N:
        print(f'Solução não obtida em {N:d} iterações')
    else:
        print(f'Solução obtida: x = {x:.6f}')


    return x

6.3. Playground interativo#

Utilize o playground interativo abaixo para testar o método da bisseção para funções não lineares quaisquer. Como dado de entrada para \(f(x)\), utilize funções escritas nos moldes de Python científico como um tipo str. Para os parâmetros do intervalo inicial e de erro, utilize float.

6.4. Aplicações#

Exemplo: Resolva o problema \(f(x) = 0\), para \(f(x) = -\text{arccos}(x) + 4\text{sen}(x) + 1.7\), no intervalo \(-0.2 \le x \le 1.0\) e \(\epsilon = 10^{-3}\).

# Chamada da função
xm = newton(-0.1,
            '-arccos(x) + 4*sin(x) + 1.7',
            '4*cos(x) - 1/(1 - x**2)**1/2',
            1e-3,
            30)
Estimativa inicial: x0 = -0.100000

+----+-------------+-------------+--------+--------------+
| i  |      x      |     f(x)    | f'(x)  |      ER      |
+----+-------------+-------------+--------+--------------+
| 0  |  0.00656144 |  0.16201076 | 3.4999 | 1.624055e+01 |
| 1  | -0.03972877 | -0.06940881 | 3.4961 | 1.165156e+00 |
| 2  | -0.01987529 |  0.02983116 | 3.499  | 9.989026e-01 |
| 3  | -0.02840088 | -0.01278928 | 3.498  | 3.001876e-01 |
| 4  | -0.02474469 |  0.00548778 | 3.4985 | 1.477564e-01 |
| 5  | -0.02631332 |  -0.0023538 | 3.4983 | 5.961326e-02 |
| 6  | -0.02564047 |  0.00100976 | 3.4984 | 2.624162e-02 |
| 7  | -0.02592911 | -0.00043314 | 3.4983 | 1.113178e-02 |
| 8  | -0.02580529 |  0.00018581 | 3.4983 | 4.798023e-03 |
| 9  |  -0.0258584 |  -7.97e-05  | 3.4983 | 2.053976e-03 |
| 10 | -0.02583562 |  3.419e-05  | 3.4983 | 8.818630e-04 |
+----+-------------+-------------+--------+--------------+
Solução obtida: x = -0.025836

Exemplo: Resolva o problema \(h(z) = 0\), para \(h(z) = -z(1-2z)^{-1} - \text{tan}(z+1)\), no intervalo \([1,8]\) e \(\epsilon = 10^{-5}\).

Como no exemplo anterior, para utilizarmos o método de Newton é preciso saber a derivada da função \(h(z)\). Vamos encontrá-la utilizando o módulo de computação simbólica sympy.

# Importa variável z como símbolo
from sympy.abc import z 
import sympy as sym

# Derivada
dh = sym.diff(-z/(1 - 2*z) - sym.tan(z+1))

# Impressão
print(dh)
-2*z/(1 - 2*z)**2 - tan(z + 1)**2 - 1 - 1/(1 - 2*z)

A partir daí, utilizamos a expressão normalmente na função.

zm = newton(5,
            'z/(1 - 2*z) - tan(z+1)',
            '-2*z/(1 - 2*z)**2 - tan(z + 1)**2 - 1 - 1/(1 - 2*z)',
            1e-5,30)
Estimativa inicial: x0 = 5.000000

+---+------------+------------+---------+--------------+
| i |     x      |    f(x)    |  f'(x)  |      ER      |
+---+------------+------------+---------+--------------+
| 0 | 4.75884953 | 0.01963205 | -1.3483 | 5.067411e-02 |
| 1 | 4.77341064 | 0.00056164 | -1.3262 | 3.050462e-03 |
| 2 | 4.77383412 | 1.173e-05  | -1.3256 | 8.870850e-05 |
| 3 | 4.77384296 |  2.4e-07   | -1.3256 | 1.852869e-06 |
+---+------------+------------+---------+--------------+
Solução obtida: x = 4.773843

Compare a quantidade de iterações obtidas com o mesmo exemplo resolvido com o algoritmo da bisseção.

6.4.1. Derivação para otimização#

O método de Newton, originalmente desenvolvido para encontrar raízes de equações, pode ser adaptado com grande eficácia para problemas de otimização. Nesse contexto, seu objetivo é encontrar os pontos críticos de uma função — ou seja, valores onde a derivada é nula e que podem corresponder a mínimos, máximos ou pontos de sela. A ideia central é usar uma aproximação quadrática da função por meio da série de Taylor de segunda ordem. O processo pode ser resumido nos passos a seguir:

  1. Função Objetivo (FO): seja \(f(x)\) uma função escalar diferenciável cujo mínimo queremos encontrar.

  2. Expansão de Taylor de Segunda Ordem: expandimos \(f(x)\) em torno de um ponto \(x^{(k)}\) usando a série de Taylor até a segunda ordem:

    \[ f(x) \approx f(x^{(k)}) + f'(x^{(k)}) (x - x^{(k)}) + \frac{1}{2} f''(x^{(k)}) (x - x^{(k)})^2, \]

    onde \(f'(x^{(k)})\) é a primeira derivada (gradiente) de \(f\) em \(x^{(k)}\) e \(f''(x^{(k)})\) é a segunda derivada (Hessiana) de \(f\) em \(x^{(k)}\).

  3. Condição de otimização: para encontrar o ponto \(x\) que minimiza \(f(x)\), devemos encontrar um ponto onde a derivada da função seja zero. Definimos a direção de descida \(t\) como \(x - x^{(k)}\):

    \[ f(x^{(k+1)}) \approx f(x^{(k)}) + f'(x^{(k)}) t + \frac{1}{2} f''(x^{(k)}) t^2 \]
  4. Derivada da expansão: para minimizar \(f(x)\), tomamos a derivada da expressão acima em relação a \(t\) e igualamos a zero:

    \[ \frac{d}{dt}\left[ f(x^{(k)}) + f'(x^{(k)}) t + \frac{1}{2} f''(x^{(k)}) t^2 \right] = 0 \]

    Uma vez que \(f(x^{(k)})\) não depende de \(t\), é uma constante. Assim obtemos:

    \[ f'(x^{(k)}) + f''(x^{(k)}) t = 0 \]
  5. _solução para a direção \(t\): resolvendo para \(t\), temos:

    \[ t = -\frac{f'(x^{(k)})}{f''(x^{(k)})} \]
  6. Atualização da Iteração: a próxima aproximação \(x^{(k+1)}\) é obtida somando \(t\) ao ponto atual \(x^{(k)}\):

    \[ x^{(k+1)} = x^{(k)} + t = x^{(k)} - \frac{f'(x^{(k)})}{f''(x^{(k)})} \]

A derivada \(f'(x^{(k)})\) indica a inclinação da função em \(x^{(k)}\). A segunda derivada \(f''(x^{(k)})\) fornece informações sobre a curvatura da função. Se \(f''(x^{(k)}) > 0\), estamos em um ponto onde a função está curvando para cima (mínimo local), enquanto \(f''(x^{(k)}) < 0\) indica um ponto onde a função está curvando para baixo (máximo local). O passo \(t\) determina a magnitude e a direção do ajuste necessário para aproximar-se do ponto de mínimo (ou máximo).

6.5. Método de Halley#

O método de Halley é uma extensão do método de Newton e utiliza até a terceira derivada da função. Enquanto o método de Newton tem convergência quadrática, o método de Halley converge mais rapidamente e alcança convergência cúbica, pois usa mais informações sobre a curvatura da função. De maneira similar ao que já fizemos, vamos derivar a fórmula partindo diretamente da expansão de Taylor:

  1. Expansão de Taylor: expandimos a função \(f(x)\) em torno de um ponto \(x^{(k)}\):

    \[ f(x) \approx f(x^{(k)}) + f'(x^{(k)}) (x - x^{(k)}) + \frac{1}{2} f''(x^{(k)}) (x - x^{(k)})^2 + \frac{1}{6} f'''(x^{(k)})(x - x^{(k)})^3, \]
  2. Estimativa da raiz: supondo que \(x^{(k+1)} = x^{(k)} + \Delta x\) é a nova estimativa da raiz, então queremos \(f(x^{(k+1)}) = 0 \). Substituindo \(\Delta x = x^{(k+1)} - x^{(k)}\) na expansão de Taylor e igualando a zero, temos:

    \[ f(x^{(k)}) + f'(x^{(k)})\Delta x + \frac{1}{2} f''(x^{(k)})(\Delta x)^2 + \frac{1}{6} f'''(x^{(k)})(\Delta x)^3 = 0 \]
  3. Isolamento do passo: a solução para \(\Delta x\) (que nos dá a atualização para \(x\)) é dada pela fórmula:

    \[ \Delta x = -\frac{f(x^{(k)})}{f'(x^{(k)})} \left( \frac{2}{2 - \frac{f(x^{(k)}) f''(x^{(k)})}{f'(x^{(k)})^2}} \right) \]
  4. Iteração do Método de Halley: a fórmula iterativa do método de Halley é então:

    \[ x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} \left( \frac{2}{2 - \frac{f(x^{(k)}) f''(x^{(k)})}{f'(x^{(k)})^2}} \right) \]

As vantagens do método de Halley sobre o método de Newton são a taxa de convergência e o alcance de uma solução mais precisa com menos iterações em comparação. Por outro lado, o cálculo das derivadas de ordem superior pode ser computacionalmente custoso e complexo, especialmente para funções complicadas. A precisão do método também depende da precisão das derivadas de segunda e terceira ordem. Portanto, erros nessas derivadas podem afetar a convergência e a precisão da solução.

6.6. Problemas propostos#

  • Crie um código genérico que implemente o algoritmo do método de Newton de modo que a derivada seja calculada diretamenta por computação simbólica, sem intervenção manual, quando for possível. Dica: use sympy.lambdify.

  • Faça uma implementação do método de Newton para otimização e uma para o método de Halley. Em seguida, incorpore a capacidade de cálculo das derivadas de maneira automática.