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

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

Авторы: Гирняк О.Р., Васильков Д.А.

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

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

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

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

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

Рассмотрим систему нелинейных уравнений

[math]F(x) = 0, F(x), x \in \mathbb{R}^{n}, (1)[/math]

и предположим, что существует вектор [math]\bar{x} \in D \subset \mathbb{R}^{n}[/math], являющийся решением системы (1).

Будем считать, что [math]F(x) = (f_1(x), f_2(x), ... f_n(x))^T[/math], причем [math]f_i(\cdot) \in C^1(D) \forall i[/math]

Разложим [math]F(x)[/math] в окрестности точки [math]\bar{x}: F(x) = F(x^0) + F'(x^0)(x-x^0) + o(||x-x^0||)[/math].

Здесь [math] F'(x) = \frac{\partial{F(x)}}{\partial{x}} = \begin{pmatrix} \frac{\partial{f_1(x_1)}}{\partial{x_1}} & \frac{\partial{f_1(x_2)}}{\partial{x_2}} & \cdots & \frac{\partial{f_1(x_n)}}{\partial{x_n}}\\ \frac{\partial{f_2(x_1)}}{\partial{x_1}} & \frac{\partial{f_2(x_2)}}{\partial{x_2}} & \cdots & \frac{\partial{f_2(x_n)}}{\partial{x_n}}\\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial{f_n(x_1)}}{\partial{x_1}} & \frac{\partial{f_n(x_2)}}{\partial{x_2}} & \cdots & \frac{\partial{f_n(x_n)}}{\partial{x_n}} \end{pmatrix} [/math]

называется матрицей Якоби, а её определитель – якобианом системы (1).

Исходное уравнение заменим следующим: [math]F(x^0) + F'(x^0)(x-x^0) = 0[/math]. Считая матрицу Якоби [math]F'(x^0)[/math] неособой, разрешим это уравнение относительно [math]x: x = x^0 - [F'(x)]^{-1}F(x^0)[/math]. И вообще положим

[math]x^{k+1} = x^{k} - [F'(x^k)]^{-1}F(x^k)[/math].

При сделанных относительно [math]F(\cdot)[/math] предположениях имеет место сходимость последовательности [math]x^{k}[/math] к решению системы со скоростью геометрической прогрессии при условии, что начальное приближение [math]x^0[/math] выбрано из достаточно малой окрестности решения [math]\bar{x}[/math].

При дополнительном предположении [math]F(\cdot) \in C^2[/math] имеет место квадратичная сходимость метода, т.е.

[math]||x^{k+1}-\bar{x}|| \le \omega||x^{k}-\bar{x}||^2[/math].

Сформулируем теорему.

Теорема. Пусть в некоторой окрестности решения [math]\bar{x}[/math] системы (1) функции [math]f_i(\cdot) \in C^2[/math] и якобиан системы отличен от нуля в этой окрестности. Тогда существует [math]\delta[/math]-окрестность точки [math]\bar{x}[/math] такая, что при любом выборе начального приближения [math]x^0[/math] из этой окрестности последовательность [math]{x_k}[/math] не выходит из неё и имеет место квадратичная сходимость этой последовательности.

Замечание 1. В качестве критерия окончания процесса итераций обычно берут условие: [math]||x^{k+1}-x^k|| \lt \varepsilon[/math].

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

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

откуда и находим искомую поправку [math]\Delta x^k[/math], а затем и следующее приближение [math]x^{k+1} = x^k + \Delta x[/math] к решению [math]\bar{x}[/math]. Очевидно, что это значительно сокращает количество арифметических операций для построения очередного приближения.

Замечание 3. Начиная с некоторого шага [math]k_0[/math] решают стационарную СЛАУ

[math]\frac{\partial F(x^{k_0})}{\partial x} \cdot \Delta x^k = -F(x^k)[/math],

Данное видоизменение носит название модифицированный метод Ньютона.

Замечание 4. (О выборе начального приближения). Пусть вектор-функция [math]\Phi(\lambda, x)[/math] такова, что [math]\Phi(1, x) = F(x)[/math], а система [math]\Phi(0, x) = 0[/math] может быть решена. Тогда разбивая [math][0,1][/math] на [math]N[/math] частей решают методом Ньютона набор из [math]N[/math] систем [math]\Phi(i/N, x) = F(x), i = 1,... N[/math], принимая для каждой следующей системы в качестве начального приближения решение предыдущей системы.

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

Входные данные:F(x) - n n-мерных функций; n-мерный вектор x^0 - начальное приближение. eps - точность

[math]x^k = x^0[/math]

Пока |x^k+1-x^k| > eps:

   Вычисляем [math]\frac{\partial F(x^k)}{\partial x} [/math] (если производные не даны, то их можно вычислить численно)
   Решаем СЛАУ(например, методом Гаусса) [math]\frac{\partial F(x^k)}{\partial x} \cdot \Delta x^k = -F(x^k)[/math]
   [math]x^{k+1} = x^k + \Delta x^k[/math]

Выходные данные - [math]x^k+1[/math]