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

Метод Ньютона
Последовательный алгоритм
Последовательная сложность	$O(n^3)$ в случае решения СЛАУ методом Гаусса
Объём входных данных	$F(x)$ - $n$ $n$-мерных функций (так же дополнительно могут быть даны производные) + $n$-мерный вектор $x^0$ - начальное приближение
Объём выходных данных	$n$-мерный вектор

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

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

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

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

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

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

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

$F(x) = 0, F(x), x \in \mathbb{R}^{n}, (1)$

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

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

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

Здесь $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}$

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

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

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

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

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

$||x^{k+1}-\bar{x}|| \le \omega||x^{k}-\bar{x}||^2$.

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

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

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

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

$\frac{\partial F(x^k)}{\partial x} \cdot \Delta x^k = -F(x^k)$,

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

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

$\frac{\partial F(x^{k_0})}{\partial x} \cdot \Delta x^k = -F(x^k)$,

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

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

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

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

• Численное вычисление Якобиана(в случае, если производные не даны) $\frac{\partial F(x^k)}{\partial x}$
• Решение СЛАУ(например, методом Гаусса) $\frac{\partial F(x^k)}{\partial x} \cdot \Delta x^k = -F(x^k)$

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

Как уже было описано в вычислительном ядре, данный метод сводится к решению СЛАУ.

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

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

---

$x^k = x^0$

DO

Вычисляем $\frac{\partial F(x^k)}{\partial x}$ (если производные не даны, то их можно вычислить численно)

Вычисляем $F(x^k)$

Решаем СЛАУ(например, методом Гаусса) $\frac{\partial F(x^k)}{\partial x} \cdot \Delta x^k = -F(x^k)$

$x^{k} = x^k + \Delta x^k$

WHILE $||\Delta x^k|| \lt \varepsilon$

---

Выходные данные: численное решение $x^k$

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

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

Выходные данные: численное решение $x^k$

2 Программная реализация алгоритма

2.2 Существующие реализации алгоритма

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

3 Литература

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

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

<references \>