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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ \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]\displaystyle{ \overline{x^*} = {(x^*_1,x^*_2,\ldots,x^*_n)}^T }[/math], который, при подстановке в исходную систему, превращает каждое уравнение в верное числовое равенство.

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

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

[math]\displaystyle{ 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]\displaystyle{ F(\cdot) }[/math] при выборе начального приближения [math]\displaystyle{ x^{(0)} }[/math] из достаточно малой окрестности решения [math]\displaystyle{ \overline{x^*} }[/math] имеет место сходимость последовательности [math]\displaystyle{ \{x^{(k)}\} }[/math]. При дополнительном предположении [math]\displaystyle{ F(\cdot) \in C^2 }[/math] имеет место квадратичная сходимость метода.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ O(n^2) + O(n^3) + O(n) = O(n^3) }[/math]

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

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

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

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

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

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

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

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

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

Литература

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

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

<references \>