Участник:N Zakharov/Метод Ньютона для решения систем нелинейных уравнений

Основные авторы описания: Александр Чернышев, Никита Захаров

Метод Ньютона
Последовательный алгоритм
Последовательная сложность	[math]O(kn^3)[/math]
Объём входных данных	[math]2n+1[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]000[/math]
Ширина ярусно-параллельной формы	[math]000[/math]

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

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

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

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

Метод был описан Исааком Ньютоном в рукописи «Об анализе уравнениями бесконечных рядов» (лат. «De analysi per aequationes numero terminorum infinitas»), адресованной в 1669 году Барроу, и в работе «Метод флюксий и бесконечные ряды» (лат. «De metodis fluxionum et serierum infinitarum») или «Аналитическая геометрия» (лат. «Geometria analytica») в собраниях трудов Ньютона, которая была написана в 1671 году. В своих работах Ньютон вводит такие понятия, как разложение функции в ряд, бесконечно малые и флюксии (производные в нынешнем понимании). Указанные работы были изданы значительно позднее: первая вышла в свет в 1711 году благодаря Уильяму Джонсону, вторая была издана Джоном Кользоном в 1736 году уже после смерти создателя. Однако описание метода существенно отличалось от его нынешнего изложения: Ньютон применял свой метод исключительно к полиномам. Он вычислял не последовательные приближения [math]x_{n}[/math], а последовательность полиномов и в результате получал приближённое решение [math]x[/math].

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

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

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

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

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

[math] \left\{\begin{matrix}f_1(x_1, ..., x_n) = 0 \\ ... \\ f_n(x_1, ..., x_n) = 0, \end{matrix}\right. [/math]

где [math]f = (f_1, ..., f_n) : \mathbb{R}^n \rightarrow \mathbb{R}^n[/math] - непрерывно дифференцируемое отображение в окрестности решения.

При начальном приближении [math]x_0[/math] и сделанных предположениях [math](k+1)[/math]-я итерация метода будет выглядеть следующим образом: [math] x_{k+1} = x_k - [f'(x_k)]^{-1}f(x_k) [/math]

Проекционные итерационные методы основываются на тех же принципах, что и методы, которые используются для задач поиска собственных значений. На самом деле, многие методы Крылова напрямую вытекают из алгоритмов Арнольди. Приближенное решение можно найти по формуле [math]\tilde{x} = x_0 + V_my[/math], а после подставления условий Петрова-Гарелкина оно принимает вид [math]\tilde{x} = x_0 + V (W^TAV)^{-1}W^Tr_0[/math], где [math]r_0 = b - Ax_0[/math], [math]V = \{v_1, v_2, ..., v_m\}[/math] - матрица, чьи столбцы являются ортонормированным базисом пространства, а матрица [math]W[/math] находится из соотношения [math]W^HV = I[/math].
[math]V^T_mAV_m = H_m[/math]
[math]G(y) = \left \| b-Ax \right \| = \left \| b-A(x_0+V_my) \right \| = \left \| r_0 -AV_my \right \|[/math]
Решение задачи минимизации
Наша цель состоит в нахождении вектора [math]y_n[/math], минимизирующего задачу [math] \| \tilde{H}_n y_n - \beta e_1 \|, \, [/math] где [math]\tilde{H}_n[/math] матрица размера (n+1) на n.
В результате применения QR-декомпозиции к матрице [math]\tilde{H}_n[/math] находим ортогональную матрицу [math]\Omega_n[/math] размера (n+1) на (n+1) и верхнюю треугольную матрицу [math]\tilde{R}_n[/math] размера (n+1) на n.
При этом в матрице [math]\tilde{R}_n[/math] последняя строка будет нулевой, то есть ее можно представить в следующем виде

[math] \tilde{R}_n = \begin{bmatrix} R_n \\ 0 \end{bmatrix}. [/math]

Перепишем нашу исходную задачу [math] \| \tilde{H}_n y_n - \beta e_1 \| = \| \Omega_n (\tilde{H}_n y_n - \beta e_1) \| = \| \tilde{R}_n y_n - \beta \Omega_n e_1 \|. [/math]

Положим [math]\beta\Omega_ne_1=\tilde{g}_n[/math],

где [math] \tilde{g}_n = \begin{bmatrix} g_n \\ \gamma_n \end{bmatrix} [/math]

Вектор [math]y_n[/math] находится из следующего соотношения [math] y_n = R_n^{-1} g_n. [/math]

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

Описание алгоритма Арнольди
В результате мы получаем матрицу [math]H^{{m+1}\times{m}}[/math]
Пусть [math]r_0[/math] - произвольный входной вектор [math]v_1[/math] = [math]r_0/\begin{Vmatrix}r_0\end{Vmatrix}[/math]

 Для [math]m = 1,2,...,k [/math]
     [math]\omega_m[/math] = [math]Av_k[/math]
     Для [math]j = 1,2, ..,m [/math]
         [math]h_{jm}=(\omega_m, v_j)[/math]
         [math]\omega_m=\omega_m-h_{jm}v_j[/math]
     Конец
     [math]h_{m+1,m}=\begin{Vmatrix}\omega_m\end{Vmatrix}[/math]
     [math]v_{m+1}=\omega_m/h_{m+1,m}[/math]
 Конец

Не будем подробно останавливаться на описании решения задачи оптимизации.
Однако, важно отметить, что расчет матриц [math]\Omega_m, R_m[/math] и вектора [math]g_m[/math] можно упростить, используя формулы для перехода от m к m+1.
Модифицируя алгоритм Арнольди соответствующим образом, мы можем ощутимо снизить затраты, необходимые для расчета матрицы [math]R_m[/math] и вектора [math]g_m[/math], необходимые для расчета вектора минимизации [math]y_m[/math].

В результате, основной интерес для нас представляет решение СЛАУ [math]R_my_m=g_m[/math]. Так как матрица [math]R_m[/math] является верхнетреугольной, то вычисление вектора [math]y_m[/math] совпадает с обратным ходом метода Гаусса.

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

В данном алгоритме можно выделить 2 базовых макрооперации:

Умножение матрицы на вектор

Скалярное произведение векторов

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

Выбор начального приближения [math]x_k[/math]

Вычисление вектора [math]F(x_k)[/math]

WHILE (условие остановки не выполнено)

Вычисление якобиана [math]J(x_k) = F'(x_k)[/math]

Решение линейной системы уравнений алгоритмом GMRES [math]J(x_k)s_k = -F(x_k)[/math]

[math]x_k = x_k + s_k[/math]

Вычисление вектора [math]F(x_k)[/math]

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

Будем оценивать сложность одной итерации метода Ньютона. Для начала оценим сложность решения СЛАУ Вычислительная сложность GMRES алгоритма зависит от сложности метода, используемого для решения задачи минимизации. Предположим, что мы нашли решение СЛАУ за [math]k[/math] итераций. В результате выполнения алгоритма Арнольди мы получаем матрицу размера [math](k+1)*k[/math]. QR-факторизация данной матрицы будет стоить нам порядка [math]O(k^2)[/math] арифметических операций. Однако, при условии, что мы будем обновлять QR-факторизацию на каждом шаге, что возможно(ссылка на википедию на тот метод), то сложность решения задачи минимизации составит даже [math]O(k)[/math]. Сложность алгоритма Арнольди при этом составляет [math]O(k*n)[/math] арифметических операций плюс сложность [math]k[/math] раз умножения матрицы на вектор - [math]O(k*n^2)[/math] Сложность оставшихся операций из алгоритма Ньютона (сложение векторов + подсчет нормы) составляет [math]O(2*n) [/math] В результате мы получаем, что на одну итерацию метода Ньютона нужно [math]O(k*n^2)[/math] операций

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

Рис. 1. Общая схема

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

Вообще говоря, метод Ньютона является последовательным методом, но на определенных шагах, таких как решение СЛАУ, приближенное вычисление Якобиана, можно применить идеи параллельных вычислений.

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

Входные данные:
[math]n^2 + n[/math] функций от [math]n[/math] переменных
[math]n[/math] вещественных чисел
Выходные данные:
[math]n[/math] вещественных чисел

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

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

Из существующих реализаций важно отметить пакет PETSc и его компонент SNES, позволяющий решать системы нелинейных уравнений методом Ньютона
https://www.mcs.anl.gov/petsc/

3 Литература

<references \>