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

Метод Ньютона решения систем нелинейных уравнений
Последовательный алгоритм
Объём выходных данных	$n$-мерный вектор

Авторы: Шохин К.О., Лебедев А.А.

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

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

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

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

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

Таким образом получаем итеративный метод : $x_{k+1} = x_k - {F^'(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) : R^n \to R, i = 1, \ldots ,n$ - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой

области $G \subset 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}} = -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}} = -F(x^{(k)})$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Пусть функция $F(x)$ непрерывно дифференцируема в открытом выпуклом множестве $G \subset R^n$. Предположим, что существуют $\overline{x^*} \in 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)$.

3 Литература

Тыртышников Е. Е. Методы численного анализа — М., Академия, 2007. - 320 c.

Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.

<references \>