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

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

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

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

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

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

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

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

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

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

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

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

[math] \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. [/math] , где [math]f_i(x_1, \ldots,x_n) : R^n \to R, i = 1, \ldots ,n[/math] - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой

области [math]G \subset R^n[/math].

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

[math]\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[/math]

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

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

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

[math] 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}[/math] – матрица Якоби.

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

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

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

Тогда [math]x^{(k+1)} = x^{(k)}+ \Delta x^{(k)}[/math].

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

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

[math]\frac{\partial{F(x^{n_i})}}{\partial{x}} = -F(x^{(k)})[/math]

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

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

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

[math]\frac{\partial F(x^{(k)})}{\partial x}\Delta x^{(k)} = -F(x^{(k)}) \qquad(*)[/math]

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

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

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

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

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

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

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

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

[math]\sum^{n}_{i=1} {D_i} + S + V[/math]

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

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

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

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

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

[math]O(n^2) + O(n^3) + O(n) = O(n^3)[/math]

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

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

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

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

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

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

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

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

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

3 Литература

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

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

<references \>