MAXimal

добавлено: 11 Jun 2008 11:14
редактировано: 21 Sep 2010 0:54

Метод Ньютона (касательных) для поиска корней

Это итерационный метод, изобретённый Исааком Ньютоном (Isaak Newton) около 1664 г. Впрочем, иногда этот метод называют методом Ньютона-Рафсона (Raphson), поскольку Рафсон изобрёл тот же самый алгоритм на несколько лет позже Ньютона, однако его статья была опубликована намного раньше.

Задача заключается в следующем. Дано уравнение:

 f(x) = 0.

Требуется решить это уравнение, точнее, найти один из его корней (предполагается, что корень существует). Предполагается, что f(x) непрерывна и дифференцируема на отрезке [a;b].

Алгоритм

Входным параметром алгоритма, кроме функции f(x), является также начальное приближение — некоторое x_0, от которого алгоритм начинает идти.

Пусть уже вычислено x_i, вычислим x_{i+1} следующим образом. Проведём касательную к графику функции f(x) в точке x = x_i, и найдём точку пересечения этой касательной с осью абсцисс. x_{i+1} положим равным найденной точке, и повторим весь процесс с начала.

Нетрудно получить следующую формулу:

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

Интуитивно ясно, что если функция f(x) достаточно "хорошая" (гладкая), а x_i находится достаточно близко от корня, то x_{i+1} будет находиться ещё ближе к искомому корню.

Скорость сходимости является квадратичной, что, условно говоря, означает, что число точных разрядов в приближенном значении x_i удваивается с каждой итерацией.

Применение для вычисления квадратного корня

Рассмотрим метод Ньютона на примере вычисления квадратного корня.

Если подставить f(x) = \sqrt{x}, то после упрощения выражения получаем:

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

Первый типичный вариант задачи — когда дано дробное число n, и нужно подсчитать его корень с некоторой точностью \rm EPS:

double n;
cin >> n;
const double EPS = 1E-15;
double x = 1;
for (;;) {
	double nx = (x + n / x) / 2;
	if (abs (x - nx) < EPS)  break;
	x = nx;
}
printf ("%.15lf", x);

Другой распространённый вариант задачи — когда требуется посчитать целочисленный корень (для данного n найти наибольшее x такое, что x^2 \le n). Здесь приходится немного изменять условие останова алгоритма, поскольку может случиться, что x начнёт "скакать" возле ответа. Поэтому мы добавляем условие, что если значение x на предыдущем шаге уменьшилось, а на текущем шаге пытается увеличиться, то алгоритм надо остановить.

int n;
cin >> 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;
}
cout << x;

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

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

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