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

e-maxx algorithm оригинал: C/C++ #algorithm #emaxx #math #numerical-methods #search

Static

Источник: e-maxx.ru/algo, страница PDF 420.

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

существует). Предполагается, что

непрерывна и дифференцируема на отрезке

Алгоритм

Входным параметром алгоритма, кроме функции

, является также начальное приближение —

некоторое

, от которого алгоритм начинает идти.

Пусть уже вычислено

, вычислим

следующим образом. Проведём касательную к графику функции

в

точке

, и найдём точку пересечения этой касательной с осью абсцисс.

положим равным найденной точке,

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

Интуитивно ясно, что если функция

достаточно "хорошая" (гладкая), а

находится достаточно близко от корня,

то

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

в приближенном значении

удваивается с каждой итерацией.

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

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

Если подставить

, то после упрощения выражения получаем:

Первый типичный вариант задачи — когда дано дробное число

, и нужно подсчитать его корень с некоторой

точностью

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

Другой распространённый вариант задачи — когда требуется посчитать целочисленный корень (для данного

найти наибольшее

такое, что

). Здесь приходится немного изменять условие останова алгоритма,

поскольку может случиться, что

начнёт "скакать" возле ответа. Поэтому мы добавляем условие, что если значение

на предыдущем шаге уменьшилось, а на текущем шаге пытается увеличиться, то алгоритм надо остановить.

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;

Наконец, приведём ещё третий вариант — для случая длинной арифметики. Поскольку число

может быть

достаточно большим, то имеет смысл обратить внимание на начальное приближение. Очевидно, что чем оно ближе к корню, тем быстрее будет достигнут результат. Достаточно простым и эффективным будет брать в качестве

начального приближения число

, где

— количество битов в числе

. Вот код на языке

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;

}

Например, этот вариант кода выполняется для числа

за

миллисекунд, а если убрать улучшенный

выбор начального приближения (просто начинать с

), то будет выполняться примерно

миллисекунд.

C# решение

auto-draft, проверить перед отправкой

using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    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;
    }
    Console.Write ("%.15lf", 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;
    }
    Console.WriteLine( x;
    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;
    }
}

C++ решение

сопоставлено/оригинал

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

Java решение

auto-draft, проверить перед отправкой

import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    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;
    }
    System.out.print ("%.15lf", x);
    int n;
    cin >> n;
    int x = 1;
    boolean decreased = false;
    for (;;) {
            int nx = (x + n / x) >> 1;
            if (x == nx || nx > x && decreased)  break;
            decreased = nx < x;
            x = nx;
    }
    System.out.println( x;
    BigInteger n; // входные данные
    BigInteger a = BigInteger.ONE.shiftLeft (n.bitLength() / 2);
    booleanean 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;
    }
}

Материал разбит как алгоритмическая задача: изучить постановку, понять асимптотику и реализовать алгоритм на выбранном языке.

Вакансии для этой задачи

Показаны активные вакансии с пересечением по тегам задачи.

Все вакансии

Активных вакансий пока нет.