Участник:Chist/Метод Ньютона

Чистяков И.А., метод Ньютона

Свойства и структура алгоритма

Общее описание алгоритма

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

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

Математическое описание алгоритма

Пусть задана некоторая функция [math]\displaystyle{ \ f(x) }[/math], её производная [math]\displaystyle{ \ f'(x) }[/math] и сетка на комплексной плоскости вида

[math]\displaystyle{ z_0 = \left( x_0 + \alpha \dfrac{x_1 - x_0}{N} \right) + i \cdot \left( y_0 + \beta \dfrac{y_1 - y_0}{M} \right), \qquad \alpha = \overline{1, N},\ \beta = \overline{1, M}, }[/math]

где [math]\displaystyle{ \ i }[/math] — мнимая единица.

Получим решения уравнения [math]\displaystyle{ \ f(z) = 0 }[/math] с помощью итерационного процесса, определяемого формулой:

[math]\displaystyle{ z_{k+1} = z_k - \dfrac{f(z_k)}{f'(z_k)}, }[/math]

где в качестве начального приближения [math]\displaystyle{ \ z_0 }[/math] перебираются всевозможные точки рассмотренной сетки.

Теоретическое обоснование

Пусть [math]\displaystyle{ \ F(x) }[/math] — оператор, отображающий линейное нормированное пространство [math]\displaystyle{ H }[/math] на линейное нормированное пространство [math]\displaystyle{ Y }[/math]. Линейный оператор [math]\displaystyle{ P }[/math], действующий из пространства [math]\displaystyle{ H }[/math] в пространство [math]\displaystyle{ Y }[/math], назовём производной оператора [math]\displaystyle{ F(x) }[/math] в точке [math]\displaystyle{ x }[/math], если

[math]\displaystyle{ \| F(x + \eta) − F(x) − P\eta \|_Y = o(\| \eta \|_H). }[/math]

Будем обозначать такой оператор [math]\displaystyle{ P }[/math] как [math]\displaystyle{ F′(x) }[/math].

Пусть [math]\displaystyle{ \ X }[/math] — решение уравнения [math]\displaystyle{ F(X) = 0 }[/math], [math]\displaystyle{ \ \Omega_a = \{ x: \| x - X \| \lt a \} }[/math]. Пусть также при некоторых [math]\displaystyle{ a \gt 0,\ a_1 \ge 0,\ a_2 \le \infty }[/math] выполнены условия:

[math]\displaystyle{ \label{requirement_one} \| (F'(X))^{-1} \|_{Y} \le a_1 \text{ при } x \in \Omega_a, }[/math]

[math]\displaystyle{ \| F(u_1) - F(u_2) - F'(u_2)(u_1 - u_2) \|_{Y} \le a_2 \|u_2 - u_1 \|_{H}^{2} \text{ при } u_1, u_2 \in \Omega_a. }[/math]

Обозначим [math]\displaystyle{ c = a_1 a_2 }[/math], [math]\displaystyle{ b = \min \{a, c^{-1}\} }[/math].

Метод Ньютона применим на основании следующей теоремы:

Теорема. Если выполнены указанные условия и начальное приближение [math]\displaystyle{ x_0 }[/math] принадлежит [math]\displaystyle{ \Omega_b }[/math], то итерационный процесс Ньютона

[math]\displaystyle{ x_{n + 1} = x_{n} - (F'(x_n))^{-1} F(x_n) }[/math]

сходится с оценкой погрешности

[math]\displaystyle{ \| x_n - X \|_{H} \le c^{-1} \left( c \| x_0 - X \|_{H} \right)^{2^n}. }[/math]

Доказательство этой теоремы можно найти в литературе.

Вычислительное ядро алгоритма

Вычислительную сложность алгоритма представляют операции вычисления значения функций в заданной точке, а так же операции вычитания и деления. При этом итерационный процесс запускается [math]\displaystyle{ N \cdot M }[/math] раз, а количество итераций в каждом процессе, вообще говоря, не определено; наиболее популярными условиями останова являются следующие:

[math]\displaystyle{ \| z_{n + 1} - z_{n} \| \lt \varepsilon, }[/math]

[math]\displaystyle{ | f(z_{n}) | \lt \varepsilon. }[/math]

Однако, для верхней оценки числа операций удобно зафиксировать число итераций в каждом процессе. Это так же позволит завершить вычисления в случае, если итерационный процесс не сходится.

Макроструктура алгоритма

Для обеспечения параллелизма сперва исходная сетка разбивается на [math]\displaystyle{ n }[/math] по возможности равных участков (например, на полосы длины [math]\displaystyle{ M }[/math]) — запускается [math]\displaystyle{ n }[/math] процессов. В каждом участке последовательно запускаются итерационные процессы для каждого узла сетки, в ходе которых:

• 1. вычисляются значения функций [math]\displaystyle{ \ f(z_k) }[/math], [math]\displaystyle{ \ f'(z_k) }[/math] (предполагается, что функции заданы аналитически);
• 2. вычисляется следующая точка:

[math]\displaystyle{ z_{k+1} = z_k - \dfrac{f(z_k)}{f'(z_k)}; }[/math]

• 3. если условие останова не выполнено, переходим к шагу 1.

Поскольку полученные в каждом процессе результаты необходимо обработать (проверить, являются ли полученные числа корнями; также желательно оставить только уникальные решения), нужно переслать полученные данные в какой-то один процесс. Проверку корней можно осуществлять на том процессоре, где они были получены, но это вызовет дополнительные сложности при работе с памятью во время последующей пересылки.

Схема реализации последовательного алгоритма

Последовательный алгоритм является частным случаем описанного алгоритма при [math]\displaystyle{ n = 1 }[/math]. При этом не возникает затрат на пересылку информации о найденных корнях к одному процессу, что позволяет более эффективно работать с памятью.

Последовательная сложность алгоритма

Для определения сложности алгоритма, необходимо выбрать класс функций, для которого будет написана программа (предполагается, что функции [math]\displaystyle{ \ f(x) }[/math], [math]\displaystyle{ \ f'(x) }[/math] заданы аналитически). Проведём исследование для многочленов степени [math]\displaystyle{ p }[/math], каждый из которых определяется вектором из [math]\displaystyle{ p + 1 }[/math] коэффициентов. Пусть выбран критерий останова, при котором каждое использование метода Ньютона включает в себя [math]\displaystyle{ iter\_num }[/math] операций.

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

Для вычисления значения многочлена степени [math]\displaystyle{ p }[/math] в точке с помощью тривиального алгоритма необходимо [math]\displaystyle{ \frac{p(p+1)}{2} + p }[/math] операций умножения. Воспользовавшись схемой Горнера, можно сократить число умножений до [math]\displaystyle{ p }[/math]:

[math]\displaystyle{ P(x) = ( \ldots (a_p x + a_{p-1}) x + a_{p-2})x + \ldots + a_1)x + a_0. }[/math]

В каждой итерации метода нужно вычислить не только значение многочлена в точке, но и значение его производной в этой точке. Это добавляет ещё [math]\displaystyle{ p - 1 }[/math] операцию. Наконец, в каждой итерации производится ровно одно деление.

Учитывая количество узлов сетки и заданное количество итераций, получим общее число операций умножения деления:

[math]\displaystyle{ p + N \cdot M \cdot iter\_num \cdot (p + (p - 1) + 1) = p + 2 \cdot N \cdot M \cdot iter\_num \cdot p, }[/math]

где первое слагаемое соответствует числу операций для вычисления коэффициентов производной многочлена.

Информационный граф

Информационный граф алгоритма представлен на рисунке 1. Каждой зелёной вершине соответствует одна итерация алгоритма (при этом запущено [math]\displaystyle{ n }[/math] различных процессов), далее все полученные результаты пересылаются в один процесс для последующей обработки и вывода (синяя вершина).

Рисунок 1. Информационный граф алгоритма.

Ресурс параллелизма алгоритма

Алгоритм является строго последовательным, однако существует возможность распараллеливания благодаря перебору. Таким образом сложность алгоритма равна [math]\displaystyle{ O(\frac{N \cdot M}{n}) }[/math].

Входные и выходные данные алгоритма

Входные данные:

• [math]\displaystyle{ p + 1 }[/math] действительное число — коэффициенты многочлена.
• Действительные числа [math]\displaystyle{ x_0 }[/math], [math]\displaystyle{ y_0 }[/math], [math]\displaystyle{ x_1 }[/math], [math]\displaystyle{ y_1 }[/math], определяющие левый нижний [math]\displaystyle{ (x_0, y_0) }[/math] и правый верхний [math]\displaystyle{ (x_1, y_1) }[/math] узлы прямоугольной сетки.
• Натуральные числа [math]\displaystyle{ M }[/math] и [math]\displaystyle{ N }[/math], задающие количество узлов сетки по горизонтали и вертикали
• Натуральное число [math]\displaystyle{ iter\_num }[/math], равное количеству итераций при каждом запуске метода Ньютона

Выходные данные:

• список всех найденных комплексных чисел [math]\displaystyle{ z }[/math], удовлетворяющих условию [math]\displaystyle{ |f(z)| \lt \varepsilon }[/math].

Свойства алгоритма

Как хорошо видно, соотношение последовательной и параллельной сложности является линейным.

При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных, так же является линейной.

Очевидно, что алгоритм является детерминированным.

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

Программная реализация алгоритма

Возможные способы и особенности параллельной реализации алгоритма

Основная идея параллелизма заключается в разбиении сетки на [math]\displaystyle{ n }[/math] участков. Идеи параллелизма можно также применять при вычислении значения многочлена в точке, однако это становится нерациональным при использовании схемы Горнера. Стоит отметить, что, за исключением перебора по сетке, в алгоритме нет каких-либо циклов.

Масштабируемость алгоритма и его реализации

На рисунке 2 показана зависимость времени работы программы от количества процессоров и степени многочлена (рассматривались уравнения вида [math]\displaystyle{ z^p - 1 = 0 }[/math]).

Рисунок 2. Масштабируемость алгоритма.

Как видно из рисунка, время выполнения пропорционально степени полинома и обратно пропорционально числу процессоров, что согласуется с теоретическими результатами. На рисунке 3 показана зависимость времени работы от числа процессоров при фиксированной степени полинома [math]\displaystyle{ p = 4 }[/math]. Можно заметить, что в данном случае полученные результаты достаточно хорошо аппроксимируются функцией [math]\displaystyle{ y(x) = 50 / x }[/math].

Рисунок 3. Зависимость от числа процессоров.

Приведённые выше результаты были получены на суперкомпьютере “Ломоносов”. Компиляция осуществлялась с помощью команды mpicc библиотеки OpenMPI, версия компилятора — gcc 1.8.4. При этом использовался ключ -std=c99, необходимый для использования описанной в языке Си реализации комплексных чисел.

Существующие реализации алгоритма

• Последовательные реализации
  • ALIAS C++. Язык реализации — C++. Распространяется бесплатно, исходные коды и примеры использования можно скачать на сайте.
  • Numerical Recipes. Язык реализации — C++. Исходный код можно найти в книге "Numerical recipes" . Бесплатно доступны предыдущие издания для языков C, Fortran77, Fortran90.
  • Numerical Mathematics — NewtonLib. Язык реализации - C, Fortran. Исходные коды доступны в качестве приложения к книге Peter Deuflhards "Newton Methods for Nonlinear Problems — Affine Invariance and Adaptive Algorithms".

• Параллельные реализации
  • Sundials. Язык реализации — C, также есть интерфейс для использования в Fortran-программах. Распространяется по лицензии BSD.
  • PETSc. Язык реализации — C, есть интерфейс для Java и Python. Распространяется по лицензии BSD.

Литература

<references \>