Метод Ньютона для систем нелинейных уравнений

Основные авторы статьи: К.О.Шохин и А.А.Лебедев,
О.Р.Гирняк и Д.А.Васильков (раздел 1.7).

Метод Ньютона решения систем нелинейных уравнений
Последовательный алгоритм
Последовательная сложность	$O(L \cdot n^3)$, если для решения СЛАУ использовать метод Гаусса, алгоритм сойдется за $L$ итераций, вычисление значения каждой функции системы и их производных в точке имеет оценку сложности $O(n)$
Объём входных данных	$n$ функций от $n$ переменных, $n^2$ частных производных этих функций по каждой переменной, $n$-мерный вектор(начальное приближение решения), $\varepsilon$ - требуемая точность, ${\rm max\_iter}$ - максимальное число итераций
Объём выходных данных	$n$-мерный вектор

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

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

Метод Ньютона решения систем нелинейных уравнений является обобщением метода Ньютона решения нелинейных уравнений, который основан на идее линеаризации. Пусть $F(x) : {\mathbb R}^1 \to {\mathbb R}^1$ - дифференцируемая функция и необходимо решить уравнение $F(x) = 0$.

Взяв некоторое $x_0$ в качестве начального приближения решения, мы можем построить линейную аппроксимацию

$F(x)$ в окрестности $x_0 : F(x_0+h) \approx F(x_0)+F^\prime(x_0)h$ и решить получающееся линейное уравнение $F(x_0 )+F^\prime (x_0 )h =0$.

Таким образом получаем итеративный метод : $x_{k+1} = x_k - {F^\prime(x_k)}^{-1}F(x_k), \ k = 0,1,\ldots$ .

Данный метод был предложен Ньютоном в 1669 году. Более точно, Ньютон оперировал только с полиномами; в выражении для $F(x+h)$ он отбрасывал члены более высокого порядка по h , чем линейные. Ученик Ньютона Рафсон в 1690 г. предложил общую форму метода (т. е. не предполагалось что $F(x)$ обязательно полином и использовалось понятие производной), поэтому часто говорят о методе Ньютона—Рафсона. Дальнейшее развитие исследований связано с именами таких известных математиков, как Фурье, Коши и другие. Например, Фурье доказал в 1818 г., что метод сходится квадратично в окрестности корня, а Коши (1829, 1847) предложил многомерное обобщение метода и использовал метод для доказательства существования решения уравнения.

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

Пусть дана система из $n$ нелинейных уравнений с $n$ неизвестными.

$\left\{\begin{matrix} f_1(x_1, \ldots, x_n) = 0, \\ f_2(x_1, \ldots, x_n) = 0, \\ \vdots \\ f_n(x_1, \ldots, x_n) = 0. \end{matrix}\right.$ , где $f_i(x_1, \ldots,x_n) : {\mathbb R}^n \to {\mathbb R}, \ i = 1,\ldots,n$ - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой области $G \subset {\mathbb R}^n$.

Запишем ее в векторном виде:

$\overline{x} = {(x_1,x_2,\ldots,x_n)}^T, F(x) ={[f_1(x),f_2(x),\ldots,f_n(x)]}^T, F(x)=0$

Требуется найти такой вектор $\overline{x^*} = {(x^*_1,x^*_2,\ldots,x^*_n)}^T$, который, при подстановке в исходную систему, превращает каждое уравнение в верное числовое равенство.

При таком подходе формула для нахождения решения является естественным обобщением формулы одномерного итеративного метода:

$x^{(k+1)} = x^{(k)} - W^{-1}(x^{(k)})\cdot F(x^{(k)}) , k=0,1,2,\ldots$, где

$W = \begin{pmatrix} \frac{\partial{f_1(x_1)}}{\partial{x_1}} & \cdots & \frac{\partial{f_1(x_n)}}{\partial{x_n}} \\ \vdots & \ddots & \vdots \\ \frac{\partial{f_n(x_1)}}{\partial{x_1}} & \cdots & \frac{\partial{f_n(x_n)}}{\partial{x_n}} \end{pmatrix}$ – матрица Якоби.

В рассмотренных предположениях относительно функции $F(\cdot)$ при выборе начального приближения $x^{(0)}$ из достаточно малой окрестности решения $\overline{x^*}$ имеет место сходимость последовательности $\{x^{(k)}\}$. При дополнительном предположении $F(\cdot) \in C^2$ имеет место квадратичная сходимость метода.

В качестве критерия окончания процесса итераций обычно берут условие $\left \| x^{(k+1)} - x^{(k)} \right \| \lt \varepsilon$, где $\varepsilon$ - требуемая точность решения.

Основная сложность метода Ньютона заключается в обращении матрицы Якоби. Вводя обозначение $\Delta x^{(k)} = x^{(k+1)} - x^{(k)}$ получаем СЛАУ для вычисления $\Delta x^{(k)}:$ $\frac{\partial{F(x^{(k)})}}{\partial{x}} \Delta x^{(k)} = -F(x^{(k)})$

Тогда $x^{(k+1)} = x^{(k)}+ \Delta x^{(k)}.$

Часто метод Ньютона модифицируют следующим образом. По ходу вычислений или заранее выбирают возрастающую последовательность чисел $n_0=0, n_1,\ldots$

При $n_i \le k \lt n_{i+1}$ вычисление $\Delta x^{(k)}$ осуществляют по следующей формуле:

$\frac{\partial{F(x^{n_i})}}{\partial{x}} \Delta x^{(k)} = -F(x^{(k)}).$

Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.

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

Основная вычислительная нагрузка алгоритма заключается в решении СЛАУ:

$\frac{\partial F(x^{(k)})}{\partial x}\Delta x^{(k)} = -F(x^{(k)}) \qquad (1)$

Для нахождения $\Delta x^{(k)}$, по которому вычисляется значение вектора $\overline{x}$ на очередной итерации: $x^{(k+1)} = x^{(k)} + \Delta x^{(k)}.$

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

Как уже было сказано, основную часть каждой итерации метода составляет нахождение матрицы Якоби и решение СЛАУ $(1)$ для нахождения $\Delta x^{(k)}.$

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

Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:

1. $n$ - число уравнений в СЛАУ
2. Max - предельное число итераций
3. $\varepsilon$ - точность вычислений
4. $x^{(0)}$ - начальное приближение

Рис.1 Блок-схема последовательного алгоритма

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

Пусть вычисление значения первой производной функции $f_i(x_1,x_2,\ldots,x_n)$ в точке $x^{(k)}$ имеет оценку сложности $D_i$, решение СЛАУ $(1)$ имеет оценку сложности $S$ (в эту оценку также включаются операции по нахождению элементов матрицы Якоби), а оценка сложности сложения векторов – $V$. Тогда оценка сложности одной итерации представляется в виде:

$\sum^{n}_{i=1} {D_i} + S + V$

Количество итераций определяется скоростью сходимости алгоритма. Скорость сходимости метода Ньютона в общем случае является квадратичной. Однако успех или неуспех применения алгоритма во многом определяется выбором начального приближения решения. Если начальное приближение выбрано достаточно хорошо и матрица системы линейных уравнений на каждой итерации хорошо обусловлена и имеет обратную матрицу, то метод Ньютона сходится к единственному в данной окрестности решению. На практике критерием работоспособности метода является число итераций: если оно оказывается большим (для большинства задач >100), то начальное приближение выбрано плохо.

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

• Вычисление значения каждой функции в заданной точке требует $O(n)$ операций.
• Вычисление значения производной каждой функции в заданной точке требует $O(n)$ операций.

В таких предположениях указанное выражение для вычисления сложности примет следующий вид:

$O(n^2) + O(n^3) + O(n) = O(n^3)$

Предположим также, что выполнены все условия для квадратичной сходимости метода, и для решения системы потребуется небольшое (<100) количество итераций.

Пусть алгоритм сойдется за $L$ итераций, тогда итоговая сложность будет равна $O(L \cdot n^3)$.

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

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

На рис. 2 изображена структура информационного графа алгоритма. Блок IN представляет собой выбор начального приближения, необходимой точности. Каждый блок IT соответствует одной итерации алгоритма, блок-схема итерации алгоритма представлена на рис. 3. Блок OUT соответствует успешному окончанию алгоритма.

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

Рис.3 Блок-схема итерации алгоритма

Информационный граф метода Ньютона может быть также представлен в виде следующего рисунка:

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

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

Хотя сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания отдельных этапов. Основной ресурс параллелизма заключен в решении СЛАУ $(1)$. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.

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

Пусть задача решается на $p$ процессорах, тогда сложность решения СЛАУ имеет оценку $O(\frac{n^3}{p})$.

Теперь каждому процессору нужно вычислить не $n$, а $\frac{n}{p}$ функций, поэтому для вычисления элементов матрицы Якоби и правой части СЛАУ потребуется $O(\frac{n^2}{p})$ операций.

Таким образом, получаем следующую оценку параллельной сложности одной итерации алгоритма: $O(\frac{n^2}{p}) + O(\frac{n^3}{p}) + O(n) = O(\frac{n^3}{p})$.

Пусть алгоритм сходится за $L$ итераций, тогда итоговая сложность будет равна $O(L \cdot \frac{n^3}{p})$.

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

При определении объема входных и выходных данных алгоритма будем придерживаться следующих предположений:

• Рассматривается 64-битная система, таким образом каждый адрес занимает 64 бита.
• Все числа, с которыми оперирует алгоритм также занимают 64 бита.
• Каждая функция будет представлена своим адресом, таким образом вклад функций в объем входных данных будет равен $n$ чисел.
• Производная каждой функции также представлена своим адресом, и вклад производных в объем входных данных будет равен $n^2$ чисел.

Входными данными алгоритма являются:

1. Функции, образующие систему уравнений.
2. Производные функций, образующих систему уравнений, если их можно задать аналитически.
3. Точность, с которой требуется найти решение.
4. Максимальное число итераций, которое может быть проделано.
5. Начальное приближение $x^{(0)}$.

Объем входных данных алгоритма равен $n ^2 + 2\cdot n + 2$ в случае, когда на вход алгоритму передаются частные производные ($n^2$ адресов частных производных $n$ функций по $n$ переменным, $n$ адресов функций, образующих систему уравнений, $n$-мерный вектор начального приближения, требуемая точность и максимальное количество итераций), когда частные производные приближенно вычисляются непосредственно в программе, объем входных данных равен $2\cdot n + 2.$

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

Объем выходных данных: $n.$

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

Соотношение последовательной и параллельной сложностей зависит от того, какие алгоритмы применяются для решения СЛАУ $(1)$. Результат работы алгоритма сильно зависит от выбора начального приближения решения. При удачном выборе алгоритм достаточно быстро сходится к искомому решению. Глобальная сходимость для многих задач отсутствует.

Существует теорема о достаточном условии сходимости метода Ньютона:

Пусть функция $F(x)$ непрерывно дифференцируема в открытом выпуклом множестве $G \subset {\mathbb R}^n$. Предположим, что существуют $\overline{x^*} \in {\mathbb R}^n$ и $r,\beta \gt 0$, такие что $N(\overline{x^*},r) \subset G , F(\overline{x^*}) = 0$ , и существует $W^{-1}(\overline{x^*})$, причем $\left \| W^{-1}(\overline{x^*})\right \| \le \beta$ и $W(x) \in Lip_{\gamma} (N(\overline{x^*},r))$. Тогда существует $\varepsilon \gt 0$ такое, что для всех $x^{(0)} \in N(\overline{x^*}, \varepsilon)$ последовательность $x^{(1)},x^{(2)},\ldots$ порождаемая итерационным процессом сходится к $\overline{x^*}$ и удовлетворяет неравенству $\left \|x^{(k+1)} - \overline{x^*}\right \| \le \beta\cdot\gamma\cdot{\left \| x^{(k)}- \overline{x^*}\right \|}^2$.

Использовались следующие обозначения:

• $N(x,r)$ - открытая окрестность радиуса $r$ с центром в точке $x$;

• $W(x) \in Lip_{\gamma}(N(x,r))$ означает, что $W(x)$ непрерывна по Липшицу с константой $\gamma$ в области $N(x,r)$.

Оставаясь в предположениях, описанных при вычислении объема входных данных, имеем объем входных данных равный $n^2 + 2\cdot n + 2 = O(n^2)$ чисел.

Оценку числа операций возьмем для реализации метода Ньютона с помощью метода Гаусса: $O(n^3)$. Тогда вычислительная мощность есть $O(n)$.