콘텐츠로 이동

근을 찾는 뉴턴 방법

이것은 Isaac Newton이 1664년경에 고안한 반복법입니다. 그러나 이 방법은 Raphson 방법이라고 불리기도 합니다. Raphson은 Newton보다 몇 년 늦게 같은 알고리즘을 고안했지만, 그의 논문이 훨씬 먼저 출판되었기 때문입니다.

문제는 다음과 같습니다. 다음 방정식이 주어졌다고 합시다.

$$f(x) = 0$$

우리는 이 방정식을 풀고자 합니다. 더 정확히 말하면, 이 방정식의 근 중 하나를 찾고자 합니다(근이 존재한다고 가정합니다). $f(x)$는 구간 $[a, b]$에서 연속이고 미분 가능하다고 가정합니다.

알고리즘

알고리즘의 입력 매개변수는 함수 $f(x)$뿐만 아니라 초기 근삿값, 즉 알고리즘이 시작하는 어떤 $x_0$로도 이루어집니다.

plot_f(x)

이미 $x_i$를 계산했다고 가정하고, $x_{i+1}$을 다음과 같이 계산합니다. 점 $x = x_i$에서 함수 $f(x)$의 그래프에 접선을 긋고, 이 접선과 $x$축의 교점을 찾습니다. $x_{i+1}$을 찾은 점의 $x$좌표로 두고, 전체 과정을 처음부터 반복합니다.

다음 공식을 얻는 것은 어렵지 않습니다.

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

먼저 $f(x)$의 도함수인 기울기 $f'(x)$를 계산한 다음, 접선의 방정식을 구하면 다음과 같습니다.

$$ y - f(x_i) = f'(x_i)(x - x_i) $$

접선은 좌표 $y = 0$, $x = x_{i+1}$에서 x축과 만납니다.

$$ - f(x_i) = f'(x_i)(x_{i+1} - x_i) $$

이제 이 방정식을 풀면 $x_{i+1}$의 값을 얻습니다.

함수 $f(x)$가 "좋고"(매끄럽고) $x_i$가 근에 충분히 가깝다면, $x_{i+1}$은 구하려는 근에 더욱 가까워진다는 것은 직관적으로 명확합니다.

수렴 속도는 이차적입니다. 이는 대략적으로 말하면 근삿값 $x_i$의 정확한 자릿수가 반복할 때마다 두 배가 된다는 뜻입니다.

제곱근 계산에의 응용

뉴턴 방법의 예로 제곱근 계산을 사용해 봅시다.

$f(x) = x^2 - n$을 대입하고 식을 정리하면 다음을 얻습니다.

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

문제의 첫 번째 전형적인 변형은 유리수 $n$이 주어지고, 그 근을 어떤 정확도 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;
}

문제의 또 다른 흔한 변형은 정수 제곱근을 계산해야 하는 경우입니다(주어진 음이 아닌 $n$에 대해 $x^2 \le n$을 만족하는 가장 큰 $x$를 찾습니다). 이 경우 $x$가 답 근처에서 "튀기" 시작할 수 있으므로, 알고리즘의 종료 조건을 약간 바꿀 필요가 있습니다. 따라서 이전 단계에서 값 $x$가 감소했는데 현재 단계에서 증가하려고 하면 알고리즘을 멈추어야 한다는 조건을 추가합니다.

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

마지막으로 세 번째 변형은 큰 수(bignum) 연산의 경우입니다. 수 $n$이 충분히 클 수 있으므로, 초기 근삿값에 주의를 기울이는 것이 좋습니다. 당연히 초기 근삿값이 근에 가까울수록 결과를 더 빨리 얻을 수 있습니다. 초기 근삿값을 수 $2^{\textrm{bits}/2}$로 잡는 것이 충분히 간단하면서도 효과적입니다. 여기서 $\textrm{bits}$는 수 $n$의 비트 수입니다. 다음은 이 변형을 보여 주는 Java 코드입니다.

public static BigInteger isqrtNewton(BigInteger n) {
    if (n.signum() == 0)
        return BigInteger.ZERO;
    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;
}

예를 들어 이 코드는 $n = 10^{1000}$에 대해 $60$밀리초 만에 실행되며, 개선된 초기 근삿값 선택을 제거하면(단순히 $1$에서 시작하면) 약 $120$밀리초가 걸립니다.

연습 문제