Метод Ньютона

Шаблон:Другие значения Метод Ньютона, алгоритм Ньютона (также известный как метод касательных) — это итерационный численный метод нахождения корня (нуля) заданной функции. Метод был впервые предложен английским физиком, математиком и астрономом Исааком Ньютоном (1643—1727). Поиск решения осуществляется путём построения последовательных приближений и основан на принципах простой итерации. Метод обладает квадратичной сходимостью. Модификацией метода является метод хорд и касательных. Также метод Ньютона может быть использован для решения задач оптимизации, в которых требуется определить ноль первой производной либо градиента в случае многомерного пространства.

Шаблон:TOChidden

Описание метода

Вывод

Чтобы численно решить уравнение $f (x) = 0$ методом простой итерации, его необходимо привести к эквивалентному уравнению: $x = φ (x)$ , где $φ$ — сжимающее отображение.

Для наилучшей сходимости метода в точке очередного приближения $x^{*}$ должно выполняться условие $φ^{'} (x^{*}) = 0$ . Решение данного уравнения ищут в виде $φ (x) = x + α (x) f (x)$ , тогда:

φ^{'} (x^{*}) = 1 + α^{'} (x^{*}) f (x^{*}) + α (x^{*}) f^{'} (x^{*}) = 0 .

В предположении, что точка приближения «достаточно близка» к корню $\tilde{x}$ и что заданная функция непрерывна $(f (x^{*}) \approx f (\tilde{x}) = 0)$ , окончательная формула для $α (x)$ такова:

α (x) = - \frac{1}{f^{'} (x)} .

С учётом этого функция $φ (x)$ определяется:

φ (x) = x - \frac{f (x)}{f^{'} (x)} .

При некоторых условиях эта функция в окрестности корня осуществляет сжимающее отображение.

Шаблон:Доказательство

В этом случае алгоритм нахождения численного решения уравнения $f (x) = 0$ сводится к итерационной процедуре вычисления:

x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})} .

По теореме Банаха последовательность приближений стремится к корню уравнения $f (x) = 0$ .

Геометрическая интерпретация

Основная идея метода заключается в следующем: задаётся начальное приближение вблизи предположительного корня, после чего строится касательная к графику исследуемой функции в точке приближения, для которой находится пересечение с осью абсцисс. Эта точка берётся в качестве следующего приближения. И так далее, пока не будет достигнута необходимая точность.

Пусть 1) вещественнозначная функция $f (x) : (a, b) \to ℝ$ непрерывно дифференцируема на интервале $(a, b)$ ;
2) существует искомая точка $x^{*} \in (a, b)$ : $f (x^{*}) = 0$ ;
3) существуют $C > 0$ и $δ > 0$ такие, что
$| f^{'} (x) | ⩾ C$ для $x \in (a, x^{*} - δ] \cup [x^{*} + δ, b)$
и $f^{'} (x) \neq 0$ для $x \in (x^{*} - δ, x^{*}) \cup (x^{*}, x^{*} + δ)$ ;
4) точка $x_{n} \in (a, b)$ такова, что $f (x_{n}) \neq 0$ .
Тогда формула итеративного приближения $x_{n}$ к $x^{*}$ может быть выведена из геометрического смысла касательной следующим образом:

f^{'} (x_{n}) = t g α_{n} = \frac{Δ y}{Δ x} = \frac{f (x_{n}) - 0}{x_{n} - x_{n + 1}} = \frac{0 - f (x_{n})}{x_{n + 1} - x_{n}},

где $α_{n}$ — угол наклона касательной прямой $y (x) = f (x_{n}) + (x - x_{n}) \cdot t g α_{n}$ к графику $f$ в точке $(x_{n}; f (x_{n}))$ .

Следовательно (в уравнении касательной прямой полагаем $y (x_{n + 1}) = 0$ ) искомое выражение для $x_{n + 1}$ имеет вид :

x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})} .

Если $x_{n + 1} \in (a, b)$ , то это значение можно использовать в качестве следующего приближения к $x^{*}$ .

Если $x_{n + 1} \notin (a, b)$ , то имеет место «перелёт» (корень $x^{*}$ лежит рядом с границей $(a, b)$ ). В этом случае надо (воспользовавшись идеей метода половинного деления) заменять $x_{n + 1}$ на $\frac{x_{n} + x_{n + 1}}{2}$ до тех пор, пока точка «не вернётся» в область поиска $(a, b)$ .

Замечания. 1) Наличие непрерывной производной даёт возможность строить непрерывно меняющуюся касательную на всей области поиска решения $(a, b)$ .
2) Случаи граничного (в точке $a$ или в точке $b$ ) расположения искомого решения $x^{*}$ рассматриваются аналогичным образом.
3) С геометрической точки зрения равенство $f^{'} (x_{n}) = 0$ означает, что касательная прямая к графику $f$ в точке $(x_{n}; f (x_{n}))$ - параллельна оси $O X$ и при $f (x_{n}) \neq 0$ не пересекается с ней в конечной части.
4) Чем больше константа $C > 0$ и чем меньше константа $δ > 0$ из пункта 3 условий, тем для $x_{n} \in (a, x^{*} - δ] \cup [x^{*} + δ, b)$ пересечение касательной к графику $f$ и оси $O X$ ближе к точке $(x^{*}; 0)$ , то есть тем ближе значение $x_{n + 1}$ к искомой $x^{*} \in (a, b)$ .

Итерационный процесс начинается с некоторого начального приближения $x_{0} \in (a, b)$ , причём между $x_{0} \in (a, b)$ и искомой точкой $x^{*} \in (a, b)$ не должно быть других нулей функции $f$ , то есть «чем ближе $x_{0}$ к искомому корню $x^{*}$ , тем лучше». Если предположения о нахождении $x^{*}$ отсутствуют, методом проб и ошибок можно сузить область возможных значений, применив теорему о промежуточных значениях.

Для предварительно заданных $ε_{x} > 0$ , $ε_{f} > 0$ итерационный процесс завершается если $| \frac{f (x_{n})}{f^{'} (x_{n})} | \approx | x_{n + 1} - x_{n} | < ε_{x}$ и $| f (x_{n + 1}) | < ε_{f}$ .
В частности, для матрицы дисплея $ε_{x}$ и $ε_{f}$ могут быть рассчитаны, исходя из масштаба отображения графика $f$ , то есть если $x_{n}$ и $x_{n + 1}$ попадают в один вертикальный, а $f (x_{n})$ и $f (x_{n + 1})$ в один горизонтальный ряд.

Алгоритм

Задается начальное приближение $x_{0}$ .
Пока не выполнено условие остановки, в качестве которого следует взять $| \frac{x_{n + 1} - x_{n}}{1 - \frac{x_{n + 1} - x_{n}}{x_{n} - x_{n - 1}}} | < ε$ , где $ε$ выполняет роль абсолютной погрешности (так как метод Ньютона является частным случаем метода простой итерации^[1]), вычисляют новое приближение: $x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}$ .

Пример

Иллюстрация применения метода Ньютона к функции $f (x) = \cos x - x^{3}$ с начальным приближением в точке $x_{0} = 0, 5$ .
Файл:Newton ex.PNG График последовательных приближений.	Файл:Newton conv.PNG График сходимости.
Согласно способу практического определения скорость сходимости может быть оценена как тангенс угла наклона графика сходимости, то есть в данном случае равна двум.

Рассмотрим задачу о нахождении положительных $x$ , для которых $\cos x = x^{3}$ . Эта задача может быть представлена как задача нахождения нуля функции $f (x) = \cos x - x^{3}$ . Имеем выражение для производной $f^{'} (x) = - \sin x - 3 x^{2}$ . Так как $\cos x ⩽ 1$ для всех $x$ и $x^{3} > 1$ для $x > 1$ , очевидно, что решение лежит между 0 и 1. Возьмём в качестве начального приближения значение $x_{0} = 0, 5$ , тогда:

\begin{matrix} x_{1} & = & x_{0} - \frac{f (x_{0})}{f^{'} (x_{0})} & = & 1, 112 141 637 097, \\ x_{2} & = & x_{1} - \frac{f (x_{1})}{f^{'} (x_{1})} & = & \underline{0,} 909 672 693 736, \\ x_{3} & = & x_{2} - \frac{f (x_{2})}{f^{'} (x_{2})} & = & \underline{0, 86} 7 263 818 209, \\ x_{4} & = & x_{3} - \frac{f (x_{3})}{f^{'} (x_{3})} & = & \underline{0, 865 47} 7 135 298, \\ x_{5} & = & x_{4} - \frac{f (x_{4})}{f^{'} (x_{4})} & = & \underline{0, 865 474 033 1} 11, \\ x_{6} & = & x_{5} - \frac{f (x_{5})}{f^{'} (x_{5})} & = & \underline{0, 865 474 033 102} . \end{matrix}

Подчёркиванием отмечены верные значащие цифры. Видно, что их количество от шага к шагу растёт (приблизительно удваиваясь с каждым шагом): от 1 к 2, от 2 к 5, от 5 к 10, иллюстрируя квадратичную скорость сходимости.

Условия применения

Файл:Newton bad.PNG

Иллюстрация расхождения метода Ньютона, применённого к функции

f (x) = x^{3} - 2 x + 2

с начальным приближением в точке

x_{0} = 0

.

Рассмотрим ряд примеров, указывающих на недостатки метода.

Контрпримеры

Если начальное приближение недостаточно близко к решению, то метод может не сойтись.

Пусть

f (x) = x^{3} - 2 x + 2 .

Тогда

x_{n + 1} = x_{n} - \frac{x_{n}^{3} - 2 x_{n} + 2}{3 x_{n}^{2} - 2} .

Возьмём ноль в качестве начального приближения. Первая итерация даст в качестве приближения единицу. В свою очередь, вторая снова даст ноль. Метод зациклится и решение не будет найдено. В общем случае построение последовательности приближений может быть очень запутанным.

Файл:Newton ex2.png

График производной функции

f (x) = x + x^{2} \sin (2 / x)

при приближении

x

к нулю справа.

Если производная не непрерывна в точке корня, то метод может расходиться в любой окрестности корня.

Рассмотрим функцию:

f (x) = {\begin{matrix} 0, & x = 0, \\ x + x^{2} \sin (\frac{2}{x}), & x \neq 0 . \end{matrix}

Тогда $f^{'} (0) = 1$ и $f^{'} (x) = 1 + 2 x \sin (2 / x) - 2 \cos (2 / x)$ всюду, кроме 0.

В окрестности корня производная меняет знак при приближении $x$ к нулю справа или слева. В то время, как $f (x) ⩾ x - x^{2} > 0$ для $0 < x < 1$ .

Таким образом $f (x) / f^{'} (x)$ не ограничено вблизи корня, и метод будет расходиться, хотя функция всюду дифференцируема, её производная не равна нулю в корне, $f$ бесконечно дифференцируема везде, кроме как в корне, а её производная ограничена в окрестности корня.

Если не существует вторая производная в точке корня, то скорость сходимости метода может быть заметно снижена.

Рассмотрим пример:

f (x) = x + x^{4 / 3} .

Тогда $f^{'} (x) = 1 + (4 / 3) x^{1 / 3}$ и $f^{″} (x) = (4 / 9) x^{- 2 / 3}$ за исключением $x = 0$ , где она не определена.

На очередном шаге имеем $x_{n}$ :

x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})} = \frac{(1 / 3) x_{n}^{4 / 3}}{(1 + (4 / 3) x_{n}^{1 / 3})} .

Скорость сходимости полученной последовательности составляет приблизительно 4/3. Это существенно меньше, нежели 2, необходимое для квадратичной сходимости, поэтому в данном случае можно говорить лишь о линейной сходимости, хотя функция всюду непрерывно дифференцируема, производная в корне не равна нулю, и $f$ бесконечно дифференцируема везде, кроме как в корне.

Если производная в точке корня равна нулю, то скорость сходимости не будет квадратичной, а сам метод может преждевременно прекратить поиск, и дать неверное для заданной точности приближение.

Пусть

f (x) = x^{2} .

Тогда $f^{'} (x) = 2 x$ и следовательно $x - f (x) / f^{'} (x) = x / 2$ . Таким образом сходимость метода не квадратичная, а линейная, хотя функция всюду бесконечно дифференцируема.

Ограничения

Пусть задано уравнение $f (x) = 0$ , где $f (x) : 𝕏 \to ℝ$ и надо найти его решение.

Ниже приведена формулировка основной теоремы, которая позволяет дать чёткие условия применимости. Она носит имя советского математика и экономиста Леонида Витальевича Канторовича (1912—1986).

Теорема Канторовича.

Если существуют такие константы $A, B, C$ , что:

$\frac{1}{| f^{'} (x) |} < A$ на $[a, b]$ , то есть $f^{'} (x)$ существует и не равна нулю;
$| \frac{f (x)}{f^{'} (x)} | < B$ на $[a, b]$ , то есть $f (x)$ ограничена;
$\exists f^{″} (x)$ на $[a, b]$ , и $| f^{″} (x) | ⩽ C ⩽ \frac{1}{2 A B}$ ;

Причём длина рассматриваемого отрезка $| a - b | < \frac{1}{A B} (1 - \sqrt{1 - 2 A B C})$ . Тогда справедливы следующие утверждения:

на $[a, b]$ существует корень $x^{*}$ уравнения $f (x) = 0 : \exists x^{*} \in [a, b] : f (x^{*}) = 0$ ;
если $x_{0} = \frac{a + b}{2}$ , то итерационная последовательность сходится к этому корню: ${x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}} \to x^{*}$ ;
погрешность может быть оценена по формуле $| x^{*} - x_{n} | ⩽ \frac{B}{2^{n - 1}} (2 A B C)^{2^{n - 1}}$ .

Из последнего из утверждений теоремы в частности следует квадратичная сходимость метода:

| x^{*} - x_{n} | ⩽ \frac{B}{2^{n - 1}} (2 A B C)^{2^{n - 1}} = \frac{1}{2} \frac{B}{2^{n - 2}} {((2 A B C)^{2^{n - 2}})}^{2} = α | x^{*} - x_{n - 1} |^{2} .

Тогда ограничения на исходную функцию $f (x)$ будут выглядеть так:

функция должна быть ограничена;
функция должна быть гладкой, дважды дифференцируемой;
её первая производная $f^{'} (x)$ равномерно отделена от нуля;
её вторая производная $f^{″} (x)$ должна быть равномерно ограничена.

Историческая справка

Метод был описан Исааком Ньютоном в рукописи «Об анализе уравнениями бесконечных рядов» (Шаблон:Lang-la), адресованной в 1669 году Барроу, и в работе «Метод флюксий и бесконечные ряды» (Шаблон:Lang-la) или «Аналитическая геометрия» (Шаблон:Lang-la) в собраниях трудов Ньютона, которая была написана в 1671 году. В своих работах Ньютон вводит такие понятия, как разложение функции в ряд, бесконечно малые и флюксии (производные в нынешнем понимании). Указанные работы были изданы значительно позднее: первая вышла в свет в 1711 году благодаря Уильяму Джонсону, вторая была издана Джоном Кользоном в 1736 году уже после смерти создателя. Однако описание метода существенно отличалось от его нынешнего изложения: Ньютон применял свой метод исключительно к полиномам. Он вычислял не последовательные приближения $x_{n}$ , а последовательность полиномов и в результате получал приближённое решение $x$ .

Этот же метод применён Ньютоном в его трактате "Математические начала" для решения уравнения Кеплера, где Ньютон предложил вполне современную аналитическую форму вычисления, записав последовательность приближений в виде переразлагаемого в каждой новой точке аналитического ряда: Шаблон:Начало цитатыряд ... сходится настолько быстро, что едва ли когда-нибудь понадобится идти в нём далее второго члена ...Шаблон:Конец цитаты

Впервые метод был опубликован в трактате «Алгебра» Джона Валлиса в 1685 году, по просьбе которого он был кратко описан самим Ньютоном. В 1690 году Джозеф Рафсон опубликовал упрощённое описание в работе «Общий анализ уравнений» (Шаблон:Lang-la). Рафсон рассматривал метод Ньютона как чисто алгебраический и ограничил его применение полиномами, однако при этом он описал метод на основе последовательных приближений $x_{n}$ вместо более трудной для понимания последовательности полиномов, использованной Ньютоном. Наконец, в 1740 году метод Ньютона был описан Томасом Симпсоном как итеративный метод первого порядка решения нелинейных уравнений с использованием производной в том виде, в котором он излагается здесь. В той же публикации Симпсон обобщил метод на случай системы из двух уравнений и отметил, что метод Ньютона также может быть применён для решения задач оптимизации путём нахождения нуля производной или градиента.

В 1879 году Артур Кэли в работе «Проблема комплексных чисел Ньютона — Фурье» (Шаблон:Lang-en) был первым, кто отметил трудности в обобщении метода Ньютона на случай мнимых корней полиномов степени выше второй и комплексных начальных приближений. Эта работа открыла путь к изучению теории фракталов.

Обобщения и модификации

Файл:Newton mod.PNG

Иллюстрация последовательных приближений метода одной касательной, применённого к функции

f (x) = e^{x} - 2

с начальным приближением в точке

x_{0} = 1, 8

.

Метод секущих

Родственный метод секущих является «приближённым» методом Ньютона и позволяет не вычислять производную. Значение производной в итерационной формуле заменяется её оценкой по двум предыдущим точкам итераций:

$f^{'} (x_{n}) \approx \frac{f (x_{n}) - f (x_{n - 1})}{x_{n} - x_{n - 1}}$ .

Таким образом, основная формула имеет вид

x_{n + 1} = x_{n} - f (x_{n}) \cdot \frac{x_{n} - x_{n - 1}}{f (x_{n}) - f (x_{n - 1})} .

Этот метод схож с методом Ньютона, но имеет немного меньшую скорость сходимости. Порядок сходимости метода равен золотому сечению — 1,618…

Замечания. 1) Для начала итерационного процесса требуются два различных значения $x_{0}$ и $x_{1}$ .
2) В отличие от «настоящего метода Ньютона» (метода касательных), требующего хранить только $x_{n}$ (и в ходе вычислений — временно $f (x_{n})$ и $f^{'} (x_{n})$ ), для метода секущих требуется сохранение $x_{n - 1}$ , $x_{n}$ , $f (x_{n - 1})$ , $f (x_{n})$ .
3) Применяется, если вычисление $f^{'} (x)$ затруднено (например, требует большого количества машинных ресурсов: времени и/или памяти).

Метод одной касательной

В целях уменьшения числа обращений к значениям производной функции применяют так называемый метод одной касательной.

Формула итераций этого метода имеет вид:

x_{n + 1} = x_{n} - \frac{1}{f^{'} (x_{0})} f (x_{n}) .

Суть метода заключается в том, чтобы вычислять производную лишь один раз, в точке начального приближения $x_{0}$ , а затем использовать это значение на каждой последующей итерации:

α (x) = α_{0} = - \frac{1}{f^{'} (x_{0})} .

При таком выборе $α_{0}$ в точке $x_{0}$ выполнено равенство:

φ^{'} (x_{0}) = 1 + α_{0} f^{'} (x_{0}) = 0,

и если отрезок, на котором предполагается наличие корня $x^{*}$ и выбрано начальное приближение $x_{0}$ , достаточно мал, а производная $φ^{'} (x)$ непрерывна, то значение $φ^{'} (x^{*})$ будет не сильно отличаться от $φ^{'} (x_{0}) = 0$ и, следовательно, график $y = φ (x)$ пройдёт почти горизонтально, пересекая прямую $y = x$ , что в свою очередь обеспечит быструю сходимость последовательности точек приближений к корню.

Этот метод является частным случаем метода простой итерации. Он имеет линейный порядок сходимости.

Метод Ньютона-Фурье

Метод Ньютона-Фурье - это расширение метода Ньютона, выведенное Жозефом Фурье для получения оценок на абсолютную ошибку аппроксимации корня, в то же время обеспечивая квадратичную сходимость с обеих сторон.

Предположим, что Шаблон:Math дважды непрерывно дифференцируема на отрезке Шаблон:Math и что Шаблон:Mvar имеет корень на этом интервале. Дополнительно положим, что Шаблон:Math на этом отрезке (например, это верно, если Шаблон:Math, Шаблон:Math, и Шаблон:Math на этом отрезке). Это гарантирует наличие единственного корня на этом отрезке, обозначим его Шаблон:Mvar. Эти рассуждения относятся к вогнутой вверх функции. Если она вогнута вниз, то заменим Шаблон:Math на Шаблон:Math, поскольку они имеют одни и те же корни.

Пусть Шаблон:Math будет правым концом отрезка, на котором мы ищем корень, а Шаблон:Math - левым концом того же отрезка. Если Шаблон:Math найдено, определим

x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})},

которое выражает обычный метод Ньютона, как описано выше. Затем определим

z_{n + 1} = z_{n} - \frac{f (z_{n})}{f^{'} (x_{n})},

где знаменатель равен Шаблон:Math, а не Шаблон:Math. Итерации Шаблон:Mvar будут строго убывающими к корню, а итерации Шаблон:Mvar - строго возрастающими к корню. Также выполняется следующее соотношение:

\lim_{n \to \infty} \frac{x_{n + 1} - z_{n + 1}}{(x_{n} - z_{n})^{2}} = \frac{f^{″} (α)}{2 f^{'} (α)}

,

таким образом, расстояние между Шаблон:Mvar и Шаблон:Mvar уменьшается квадратичным образом.

Многомерный случай

Обобщим полученный результат на многомерный случай.

Пусть необходимо найти решение системы:

{\begin{matrix} f_{1} (x_{1}, x_{2}, \dots, x_{n}) & = & 0, \\ \dots \\ f_{m} (x_{1}, x_{2}, \dots, x_{n}) & = & 0 . \end{matrix}

Выбирая некоторое начальное значение ${\vec{x}}^{[0]}$ , последовательные приближения ${\vec{x}}^{[j + 1]}$ находят путём решения систем уравнений:

f_{i} + \sum_{k = 1}^{n} \frac{\partial f_{i}}{\partial x_{k}} (x_{k}^{[j + 1]} - x_{k}^{[j]}) = 0, i = 1, 2, \dots, m,

где ${\vec{x}}^{[j]} = (x_{1}^{[j]}, x_{2}^{[j]}, \dots, x_{n}^{[j]}), j = 0, 1, 2, \dots$ .

Применительно к задачам оптимизации

Пусть необходимо найти минимум функции многих переменных $f (\vec{x}) : ℝ^{n} \to ℝ$ . Эта задача равносильна задаче нахождения нуля градиента $\nabla f (\vec{x})$ . Применим изложенный выше метод Ньютона:

\nabla f ({\vec{x}}^{[j]}) + H ({\vec{x}}^{[j]}) ({\vec{x}}^{[j + 1]} - {\vec{x}}^{[j]}) = 0, j = 1, 2, \dots, n,

где $H (\vec{x})$ — гессиан функции $f (\vec{x})$ .

В более удобном итеративном виде это выражение выглядит так:

{\vec{x}}^{[j + 1]} = {\vec{x}}^{[j]} - H^{- 1} ({\vec{x}}^{[j]}) \nabla f ({\vec{x}}^{[j]}) .

В случае квадратичной функции метод Ньютона находит экстремум за одну итерацию.

Нахождение матрицы Гессе связано с большими вычислительными затратами, и зачастую не представляется возможным. В таких случаях альтернативой могут служить квазиньютоновские методы, в которых приближение матрицы Гессе строится в процессе накопления информации о кривизне функции. Также можно использовать метод коллинеарных градиентов, метод первого порядка, который может приблизительно (но с высокой точностью) находить шаги $- H^{- 1} \nabla f$ метода Ньютона.

Метод Ньютона — Рафсона

Метод Ньютона — Рафсона является улучшением метода Ньютона нахождения экстремума, описанного выше. Основное отличие заключается в том, что на очередной итерации каким-либо из методов одномерной оптимизации выбирается оптимальный шаг:

{\vec{x}}^{[j + 1]} = {\vec{x}}^{[j]} - λ_{j} H^{- 1} ({\vec{x}}^{[j]}) \nabla f ({\vec{x}}^{[j]}),

где $λ_{j} = \arg \min_{λ} f ({\vec{x}}^{[j]} - λ H^{- 1} ({\vec{x}}^{[j]}) \nabla f ({\vec{x}}^{[j]})) .$ Для оптимизации вычислений применяют следующее улучшение: вместо того, чтобы на каждой итерации заново вычислять гессиан целевой функции, ограничиваются начальным приближением $H (f ({\vec{x}}^{[0]}))$ и обновляют его лишь раз в $m$ шагов, либо не обновляют вовсе.

Применительно к задачам о наименьших квадратах

На практике часто встречаются задачи, в которых требуется произвести настройку свободных параметров объекта или подогнать математическую модель под реальные данные. В этих случаях появляются задачи о наименьших квадратах:

F (\vec{x}) = ‖ \vec{f} (\vec{x}) ‖ = \sum_{i = 1}^{m} f_{i}^{2} (\vec{x}) = \sum_{i = 1}^{m} (φ_{i} (\vec{x}) - ℱ_{i})^{2} \to \min .

Эти задачи отличаются особым видом градиента и матрицы Гессе:

\nabla F (\vec{x}) = 2 J^{T} (\vec{x}) \vec{f} (\vec{x}),

H (\vec{x}) = 2 J^{T} (\vec{x}) J (\vec{x}) + 2 Q (\vec{x}), Q (\vec{x}) = \sum_{i = 1}^{m} f_{i} (\vec{x}) H_{i} (\vec{x}),

где $J (\vec{x})$ — матрица Якоби вектор-функции $\vec{f} (\vec{x})$ , $H_{i} (\vec{x})$ — матрица Гессе для её компоненты $f_{i} (\vec{x})$ .

Тогда очередной шаг $\vec{p}$ определяется из системы:

[J^{T} (\vec{x}) J (\vec{x}) + \sum_{i = 1}^{m} f_{i} (\vec{x}) H_{i} (\vec{x})] \vec{p} = - J^{T} (\vec{x}) \vec{f} (\vec{x}) .

Метод Гаусса — Ньютона

Шаблон:Main Метод Гаусса — Ньютона строится на предположении о том, что слагаемое $J^{T} (\vec{x}) J (\vec{x})$ доминирует над $Q (\vec{x})$ . Это требование не соблюдается, если минимальные невязки велики, то есть если норма $‖ \vec{f} (\vec{x}) ‖$ сравнима с максимальным собственным значением матрицы $J^{T} (\vec{x}) J (\vec{x})$ . В противном случае можно записать:

J^{T} (\vec{x}) J (\vec{x}) \vec{p} = - J^{T} (\vec{x}) \vec{f} (\vec{x}) .

Таким образом, когда норма $‖ Q (\vec{x}) ‖$ близка к нулю, а матрица $J (\vec{x})$ имеет полный столбцевой ранг, шаг $\vec{p}$ мало отличается от ньютоновского (с учётом $Q (\vec{x})$ ), и метод может достигать квадратичной скорости сходимости, хотя вторые производные и не учитываются. Улучшением метода является алгоритм Левенберга — Марквардта, основанный на эвристических соображениях.

Обобщение на комплексную плоскость

До сих пор в описании метода использовались функции, осуществляющие отображения в пределах множества вещественных значений. Однако метод может быть применён и для нахождения нуля функции комплексной переменной. При этом процедура остаётся неизменной:

z_{n + 1} = z_{n} - \frac{f (z_{n})}{f^{'} (z_{n})} .

Особый интерес представляет выбор начального приближения $z_{0}$ . Ввиду того, что функция может иметь несколько нулей, в различных случаях метод может сходиться к различным значениям, и вполне естественно возникает желание выяснить, какие области обеспечат сходимость к тому или иному корню. Этот вопрос заинтересовал Артура Кэли ещё в 1879 году, однако разрешить его смогли лишь в 70-х годах двадцатого столетия с появлением вычислительной техники. Оказалось, что на пересечениях этих областей (их принято называть областями притяжения) образуются так называемые фракталы — бесконечные самоподобные геометрические фигуры.

Ввиду того, что Ньютон применял свой метод исключительно к полиномам, фракталы, образованные в результате такого применения, обрели название фракталов Ньютона или бассейнов Ньютона. Шаблон:-

Реализация

Scala

object NewtonMethod {

  val accuracy = 1e-6

  @tailrec
  def method(x0: Double, f: Double => Double, dfdx: Double => Double, e: Double): Double = {
    val x1 = x0 - f(x0) / dfdx(x0)
    if (abs(x1 - x0) < e) x1
    else method(x1, f, dfdx, e)
  }

  def g(C: Double) = (x: Double) => x*x - C

  def dgdx(x: Double) = 2*x

  def sqrt(x: Double) = x match {
    case 0 => 0
    case x if (x < 0) => Double.NaN
    case x if (x > 0) => method(x/2, g(x), dgdx, accuracy) 
  }
}

Python

from math import sin, cos
from typing import Callable
import unittest


def newton(f: Callable[[float], float], f_prime: Callable[[float], float], x0: float, 
	eps: float=1e-7, kmax: int=1e3) -> float:
	"""
	solves f(x) = 0 by Newton's method with precision eps
	:param f: f
	:param f_prime: f'
	:param x0: starting point
	:param eps: precision wanted
	:return: root of f(x) = 0
	"""
	x, x_prev, i = x0, x0 + 2 * eps, 0
	
	while abs(x - x_prev) >= eps and i < kmax:
		x, x_prev, i = x - f(x) / f_prime(x), x, i + 1

	return x


class TestNewton(unittest.TestCase):
	def test_0(self):
		def f(x: float) -> float:
			return x**2 - 20 * sin(x)


		def f_prime(x: float) -> float:
			return 2 * x - 20 * cos(x)


		x0, x_star = 2, 2.7529466338187049383

		self.assertAlmostEqual(newton(f, f_prime, x0), x_star)


if __name__ == '__main__':
	unittest.main()

PHP

<?php
// PHP 5.4
function newtons_method(
	$a = -1, $b = 1, 
	$f = function($x) {
	
		return pow($x, 4) - 1;
	
	},
	$derivative_f = function($x) {

		return 4 * pow($x, 3);
	
	}, $eps = 1E-3) {

        $xa = $a;
        $xb = $b;

        $iteration = 0;

        while (abs($xb) > $eps) {

            $p1 = $f($xa);
            $q1 = $derivative_f($xa);
            $xa -= $p1 / $q1;
            $xb = $p1;
            ++$iteration;

        }

        return $xa;

}

Octave

function res = nt()
  eps = 1e-7;
  x0_1 = [-0.5,0.5];
  max_iter = 500;
  xopt = new(@resh, eps, max_iter);   
  xopt
endfunction
function a = new(f, eps, max_iter)
  x=-1;
  p0=1;
  i=0;
 while (abs(p0)>=eps)
    [p1,q1]=f(x);
    x=x-p1/q1;
   p0=p1;
   i=i+1;
 end
 i
 a=x;
endfunction
function[p,q]= resh(x)   % p= -5*x.^5+4*x.^4-12*x.^3+11*x.^2-2*x+1;
   p=-25*x.^4+16*x.^3-36*x.^2+22*x-2;
   q=-100*x.^3+48*x.^2-72*x+22;
endfunction

Delphi

// вычисляемая функция
function fx(x: Double): Double;
begin
  Result := x * x - 17;
end;

// производная функция от f(x)
function dfx(x: Double): Double;
begin
  Result := 2 * x;
end;

function solve(fx, dfx: TFunc<Double, Double>; x0: Double): Double;
const
  eps = 0.000001;
var
  x1: Double;
begin
  x1 := x0 - fx(x0) / dfx(x0); // первое приближение
  while (Abs(x1-x0) > eps) do begin // пока не достигнута точность 0.000001
    x0 := x1;
    x1 := x1 - fx(x1) / dfx(x1); // последующие приближения
  end;
  Result := x1;
end;

// Вызов
solve(fx, dfx,4);

C++

#include <iostream>
#include <cmath>

constexpr double eps = 0.000001;
double fx(double x) { return x * x - 17; } // вычисляемая функция
double dfx(double x) { return 2 * x; }     // производная функции

using function = double (*)(double x); // задание типа function

double solve(function fx, function dfx, double x0) {
  double x1  = x0 - fx(x0) / dfx(x0); // первое приближение
  while (fabs(x1 - x0) > eps) {       // пока не достигнута точность 0.000001
    x0 = x1;
    x1 = x0 - fx(x0) / dfx(x0);       // последующие приближения
  }
  return x1;
}

int main () {
  std::cout << solve(fx, dfx, 4) << "\n"; // вывод на экран
  return 0;
}

C

typedef double (*function)(double x);

double TangentsMethod(function f, function df, double xn, double eps) {
   double x1  = xn - f(xn)/df(xn);
   double x0 = xn;
   while(abs(x0-x1) > eps) {
      x0 = x1;
      x1 = x1 - f(x1)/df(x1);
   }
   return x1;
}

//Выбор начального приближения
xn = MyFunction(A)*My2Derivative(A) > 0 ? B : A;

double MyFunction(double x) { return (pow(x, 5) - x - 0.2); } //Ваша функция
double MyDerivative(double x) { return (5*pow(x, 4) - 1); } //Первая производная
double My2Derivative(double x) { return (20*pow(x, 3)); } //Вторая производная

//Пример вызова функции
double x = TangentsMethod(MyFunction, MyDerivative, xn, 0.1)

Haskell

import Data.List ( iterate' )

main :: IO ()
main = print $ solve (\ x -> x * x - 17) ( * 2) 4

-- Функция solve универсальна для всех вещественных типов значения которых можно сравнивать.
solve = esolve 0.000001

esolve epsilon func deriv x0 = fst . head $ dropWhile pred pairs
  where
    pred (xn, xn1) = (abs $ xn - xn1) > epsilon -- Функция pred определяет достигнута ли необходимая точность.
    next xn = xn - func xn / deriv xn -- Функция next вычисляет новое приближение.
    iters   = iterate' next x0        -- Бесконечный список итераций.
    pairs   = zip iters (tail iters)  -- Бесконечный список пар итераций вида: [(x0, x1), (x1, x2) ..].

Fortran

!   Main program
    REAL*8:: Xbeg, F, D1F, error  ! Имена переменных в главной программе и подпрограмме могут отличаться    
    INTEGER  Niter, Ncalc         ! Xbeg - начальное значение, F - функция, D1F - её производная, error - остаточная ошибка
        ***                       ! Niter - заданное число итераций, Ncalc - число выполненных итераций до достижения погрешности
    CALL NEWTON(Xbeg, Niter, F, D1F, Ncalc, error)
        ***
C======================================================
    SUBROUTINE NEWTON(X0, Nmax, Func, D1Func, Nevl, rer) ! Простейший вариант устойчиво работающей программы для нахождения корня  без второй производной                                                               

	REAL*8:: X0, X1, XB, q, Func, D1Func, rer, eps       ! Итог вычисления будет записан в переменную Х0
	INTEGER  Nmax, Nevl
	
	IF(Nmax*(1000-Nmax).LE.0) Nmax=1000                 ! Защита от дурака
	                Nevl=1; XB=X0
	DO I=1, Nmax
	IF(Func(X0).EQ.0.)            EXIT
    IF(D1Func(X0).EQ.0.)          THEN
    Print *, 'Error from NEWTON: D1Func=', D1Func(X0), '  X=', X0, '  I=', I
                                  EXIT
    END IF                   
		
	X1=X0-Func(X0)/D1Func(X0) 
	
	q=abs(D1Func(X0));  q=abs(1.-q)/q
	eps=MAX(rer, epsilon(X0))                           ! epsilon(X0) - машинная точность; выбирается, если rer=0.
	IF(abs(X0-X1).LE.q*eps)       EXIT
	
	X0=X1
	END DO
    
    IF(abs(Func(X0)).GE.abs(Func(XB))) PAUSE 'Error from NEWTON: Change the X0!'

	If(I.ne.Nmax+1) Nevl=I
	If(I.eq.Nmax+1) Nevl=Nmax
	
	END SUBROUTINE

PascalABC.Net

##
function solve(fx: real-> real; dfx: real-> real; x0: real): real;

const
  eps = 10e-12;
begin
  var x1 := x0 - fx(x0) / dfx(x0); // первое приближение
  while (Abs(x1 - x0) > eps) do // пока не достигнута точность 
    (x0, x1) := (x1, x1 - fx(x1) / dfx(x1)); // последующие приближения
  Result := x1;
end;

print(solve(x -> x * x - 17, x -> 2 * x, 4)); //Вызов функции и вывод.

Литература

См. также

Примечания

Шаблон:Примечания

Ссылки

Шаблон:Родственные проекты

Шаблон:Методы оптимизации

↑ Шаблон:Cite web

[1] Шаблон:Cite web

[1]