Método de Newton para encontrar raízes

Este é um método iterativo inventado por Isaac Newton por volta de 1664. No entanto, esse método também é chamado de método Raphson, já que Raphson inventou o mesmo algoritmo alguns anos depois de Newton, mas seu artigo foi publicado muito antes.

A tarefa é a seguinte. Dada a equação:

$$f(x) = 0$$

Queremos resolver a equação. Mais precisamente, queremos encontrar uma de suas raízes (presume-se que a raiz exista). Assumindo que $f(x)$ seja contínua e diferenciável em um intervalo $[a, b]$.

Algoritmo

Os parâmetros de entrada do algoritmo consistem não apenas na função $f(x)$ mas também na aproximação inicial - algum $x_0$, com o qual o algoritmo inicia.

Suponha que já calculamos $x_i$, agora calcularemos $x_{i+1}$ da seguinte maneira. Desenhe a tangente no gráfico da função $f(x)$ no ponto $x = x_i$, e encontre o ponto de interseção dessa tangente com o eixo $x$. $x_{i+1}$ é definido igual à coordenada $x$ do ponto encontrado, e repetimos todo o processo desde o início.

Obtemos a seguinte fórmula:

$$ x_{i+1} = x_i - \frac{f(x_i)}{f^\prime(x_i)} $$

É intuitivo que se a função $f(x)$ for "boa"(simples), e $x_i$ estiver próximo o suficiente da raiz, então $x_{i+1}$ estará ainda mais próximo da raiz desejada .

A taxa de convergência é quadrática, isso significa que o número de dígitos exatos no valor aproximado de $x_i$ dobra a cada iteração.

Aplicação para calcular a raiz quadrada

Vamos usar o cálculo da raiz quadrada como um exemplo do método de Newton.

Se substituirmos $f(x) = x^2 - n$, depois de simplificar a expressão, obtemos:

$$ x_{i+1} = \frac{x_i + \frac{n}{x_i}}{2} $$

A primeira variante típica do problema é quando um número racional $n$ é fornecido e sua raiz deve ser calculada com alguma precisão eps:

double sqrt_newton(double n) {
    const double eps = 1E-15;
    double x = 1;
    for (;;) {
        double nx = (x + n / x) / 2;
        if (abs(x - nx) < eps)
            break;
        x = nx;
    }
    return x;
}

Outra variante comum do problema é quando precisamos calcular a raiz inteira (para o $n$ fornecido, encontre o maior $x$ de modo que $x^2 \le n$). Aqui é necessário alterar levemente a condição de terminação do algoritmo, pois pode acontecer que $x$ comece a "pular" perto da resposta. Portanto, adicionamos uma condição de que, se o valor $x$ diminuiu na etapa anterior e ela tenta aumentar na etapa atual, o algoritmo deve ser parado.

int isqrt_newton(int n) {
    int x = 1;
    bool decreased = false;
    for (;;) {
        int nx = (x + n / x) >> 1;
        if (x == nx || nx > x && decreased)
            break;
        decreased = nx < x;
        x = nx;
    }
    return x;
}

Finalmente, recebemos a terceira variante - para o caso da aritmética bignum ou modular rápida. Como o número $n$ pode ser grande o suficiente, faz sentido prestar atenção à aproximação inicial. Obviamente, quanto mais próximo da raiz, mais rápido o resultado será alcançado. É bastante simples e eficaz usar a aproximação inicial como o número $2^{\textrm{bits}/2}$, onde $\textrm{bits}$ é o número de bits no número $n$. Aqui está o código em Java que demonstra essa variante:

public static BigInteger isqrtNewton(BigInteger n) {
    BigInteger a = BigInteger.ONE.shiftLeft(n.bitLength() / 2);
    boolean p_dec = false;
    for (;;) {
        BigInteger b = n.divide(a).add(a).shiftRight(1);
        if (a.compareTo(b) == 0 || a.compareTo(b) < 0 && p_dec)
            break;
        p_dec = a.compareTo(b) > 0;
        a = b;
    }
    return a;
}

Por exemplo, esse código é executado em $60$ milissegundos quando $n = 10^{1000}$, e, se removermos a seleção aprimorada da aproximação inicial (apenas começando com $1$), ele será executado em cerca de $120$ milisegundos.

Problemas