Метод Ньютона для пошуку коренів

Це ітераційний метод, винайдений Ісааком Ньютоном приблизно у 1664 році. Однак цей метод інколи також називають методом Рафсона, оскільки Рафсон винайшов той самий алгоритм через кілька років після Ньютона, але його статтю опублікували значно раніше.

Задача формулюється так. Маємо таке рівняння:

$$f(x) = 0$$

Ми хочемо розв'язати це рівняння. Точніше, ми хочемо знайти один із його коренів (вважаємо, що корінь існує). Припускаємо, що $f(x)$ є неперервною і диференційовною на відрізку $[a, b]$.

Коли підходить цей алгоритм?

• Функція $f(x)$ неперервна й диференційовна, а похідну $f'(x)$ можна обчислити (аналітично чи чисельно)?
• Потрібна швидка квадратична збіжність і є достатньо близьке початкове наближення $x_0$? (якщо похідна недоступна або поведінка непевна → Бінарний пошук кореня)
• Шукаєте саме корінь рівняння $f(x)=0$, а не екстремум функції? (якщо потрібен екстремум унімодальної функції → Тернарний пошук)

Алгоритм

Вхідними параметрами алгоритму є не лише функція $f(x)$, а й початкове наближення — деяке $x_0$, з якого алгоритм починає роботу.

Припустимо, що ми вже обчислили $x_i$; обчислимо $x_{i+1}$ так. Проведемо дотичну до графіка функції $f(x)$ у точці $x = x_i$ і знайдемо точку перетину цієї дотичної з віссю $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)$$

Дотична перетинає вісь $x$ у точці з координатами $y = 0$ і $x = x_{i+1}$:

$$- 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:

C++

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;
}

Python

def sqrt_newton(n: float) -> float:
    eps = 1e-15
    x = 1.0
    while True:
        nx = (x + n / x) / 2
        if abs(x - nx) < eps:
            break
        x = nx
    return x

TypeScript

function sqrtNewton(n: number): number {
  const eps = 1e-15;
  let x = 1;
  for (;;) {
    const nx = (x + n / x) / 2;
    if (Math.abs(x - nx) < eps) break;
    x = nx;
  }
  return x;
}

Go

func sqrtNewton(n float64) float64 {
	const eps = 1e-15
	x := 1.0
	for {
		nx := (x + n/x) / 2
		if math.Abs(x-nx) < eps {
			break
		}
		x = nx
	}
	return x
}

Інший поширений варіант задачі — коли нам потрібно обчислити цілочисловий корінь (для заданого $n$ знайти найбільше $x$ таке, що $x^2 \le n$). Тут необхідно трохи змінити умову завершення алгоритму, оскільки може статися, що $x$ почне «стрибати» біля відповіді. Тому ми додаємо умову: якщо значення $x$ зменшилося на попередньому кроці, а на поточному кроці намагається зрости, то алгоритм потрібно зупинити.

C++

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;
}

Python

def isqrt_newton(n: int) -> int:
    x = 1
    decreased = False
    while True:
        nx = (x + n // x) >> 1
        if x == nx or (nx > x and decreased):
            break
        decreased = nx < x
        x = nx
    return x

TypeScript

function isqrtNewton(n: number): number {
  let x = 1;
  let decreased = false;
  for (;;) {
    // ділення з округленням донизу для додатних чисел
    const nx = Math.floor((x + Math.floor(n / x)) / 2);
    if (x === nx || (nx > x && decreased)) break;
    decreased = nx < x;
    x = nx;
  }
  return x;
}

Go

func isqrtNewton(n int) int {
	x := 1
	decreased := false
	for {
		nx := (x + n/x) >> 1
		if x == nx || (nx > x && decreased) {
			break
		}
		decreased = nx < x
		x = nx
	}
	return x
}

Нарешті, розгляньмо третій варіант — для випадку довгої арифметики. Оскільки число $n$ може бути доволі великим, має сенс звернути увагу на початкове наближення. Очевидно, що чим воно ближче до кореня, тим швидше буде досягнуто результату. Достатньо просто й ефективно взяти початкове наближення у вигляді числа $2^{\textrm{bits}/2}$, де $\textrm{bits}$ — кількість бітів у числі $n$. Ось код мовою Java, який демонструє цей варіант:

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;
}

Наприклад, цей код виконується за $60$ мілісекунд для $n = 10^{1000}$, а якщо прибрати покращений вибір початкового наближення (просто почавши з $1$), то він виконуватиметься приблизно за $120$ мілісекунд.

Задачі для практики

• UVa 10428 - The Roots

Відеоматеріали