Участник:Арутюнов А.В.

Решение системы нелинейных уравнений методом Ньютона
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(n^3) }[/math] - одна итерация
Объём входных данных	[math]\displaystyle{ n^2 + n }[/math] функций от [math]\displaystyle{ n }[/math] переменных, [math]\displaystyle{ n }[/math] мерный вектор, [math]\displaystyle{ x^0 }[/math] - начальное приближение, ε - точность решения.
Объём выходных данных	[math]\displaystyle{ n }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ }[/math]

Автор описания: Арутюнов А.В., Жилкин А.С.

Решение системы нелинейных уравнений методом Ньютона

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

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

Это итерационный численный метод нахождения корня (нуля) заданной функции. Модификацией метода является метод хорд и касательных.

Классический метод Ньютона или касательных заключается в том, что если [math]\displaystyle{ x_n }[/math] — некоторое приближение к корню [math]\displaystyle{ x_* }[/math] уравнения [math]\displaystyle{ f(x)=0 }[/math], [math]\displaystyle{ f(x) \in C^1 }[/math], то следующее приближение определяется как корень касательной [math]\displaystyle{ f(x) }[/math] к функции, проведенной в точке [math]\displaystyle{ x_n }[/math].

Метод был предложен Исааком Ньютоном в 1669 году. Поиск решения осуществляется путём построения последовательных приближений и основан на принципах простой итерации.

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

Идея метода заключается в линеаризации уравнений системы

[math]\displaystyle{ \left\{\begin{matrix}f_1(x_1, ..., x_n) = 0 \\ f_2(x_1,...x_n)=0, \\ ..., \\ f_n(x_1, ..., x_n) = 0, \end{matrix}\right. }[/math]

Что позволяет свести исходную задачу СНУ(система нелинейных уравнений) к многократному решению системы линейных уравнений.

Пусть известно приближение [math]\displaystyle{ (x_i)^{(k)} }[/math] решения системы нелинейных уравнений [math]\displaystyle{ (x_i)^* }[/math].Введем в рассмотрение поправку [math]\displaystyle{ \Delta x_i }[/math] как разницу между решением и его приближением: [math]\displaystyle{ \Delta x_i = (x_i)^* -(x_i)^{(k)} \Rightarrow (x_i)^*=(x_i)^{(k)} + \Delta x_i, i=\overline{(1,n)} }[/math]

Подставим полученное выражение для [math]\displaystyle{ (x_i)^* }[/math] в исходную систему.

[math]\displaystyle{ \left\{\begin{matrix} f_1((x_1)^{(k)} + \Delta x_1, (x_2)^{(k)}+ \Delta x_2, ..., (x_n)^{(k)} + \Delta x_n) = 0, \\ f_2((x_1)^{(k)} + \Delta x_1, (x_2)^{(k)} + \Delta x_2, ..., (x_n)^{(k)} + \Delta x_n) = 0, \\ ... \\ f_n((x_1)^{(k)} + \Delta x_1, (x_2)^{(k)} + \Delta x_2, ..., (x_n)^{(k)} + \Delta x_n) = 0, \end{matrix}\right. }[/math]

Неизвестными в этой системе нашейных уравнений являются поправки [math]\displaystyle{ \Delta x_i }[/math]. Для определения [math]\displaystyle{ \Delta x_i }[/math] нужно решить эту систему. Но решить эту задачу так же сложно, как и исходную. Однако эту систему можно линеаризовать и, решив её, получить приближённые значения поправок [math]\displaystyle{ \Delta x_i }[/math] для нашего приближения, т.е. [math]\displaystyle{ \Delta (x_i)^{(k)} }[/math]. Эти поправки не позволяют сразу получить точное решение [math]\displaystyle{ (x_i)^* }[/math], но дают возможность приблизиться к решению, - получить новое приближение решения [math]\displaystyle{ ((x_i)^{(k+1)} = ((x_i)^{(k)} + ( \Delta x_i)^{(k)}, i = \overline{(1,n)} }[/math]

Для линеаризации системы следует разложить функцию [math]\displaystyle{ f_i }[/math] в ряды Тейлора в окрестности [math]\displaystyle{ (x_i)^k }[/math], ограничиваясь первыми дифференциалами. Полученная система имеет вид:

[math]\displaystyle{ \sum_{i=1}^n\frac{ \partial f_i({x_1}^{(k)}, {x_2}^{(k)}, ..., {x_n}^{(k)})}{ \partial x_i} \Delta {x_i}^{(k)}= -f_j({x_1}^{(k)}, {x_2}^{(k)}, ..., {x_n}^{(k)}), j=\overline{(1,n)} }[/math]

Все коэффициенты этого уравнения можно вычислить, используя последнее приближение решения [math]\displaystyle{ (x_i)^{(k)} }[/math]. Для решения системы линейных уравнений при [math]\displaystyle{ n=2,3 }[/math] можно использовать формулы Крамера, при большей размерности системы [math]\displaystyle{ n }[/math] - метод исключения Гаусса.

Значения поправок используются для оценки достигнутой точности решения. Если максимальная по абсолютной величине поправка меньше заданной точности [math]\displaystyle{ \varepsilon }[/math], расчет завершается. Таким образом, условие окончания расчета:

[math]\displaystyle{ \delta =\min_{ i=\overline{(1,n)}}|\Delta {x_i}^{(k)}| }[/math]

Можно использовать и среднее значение модулей поправок:

[math]\displaystyle{ \delta = \frac{1}{n}\sum_{i=1}^n |\Delta x_i|\lt \varepsilon }[/math]

В матричной форме систему можно записать как:

[math]\displaystyle{ W(\Delta X^k)*X^{(k)} = -F(X^k) }[/math]

где [math]\displaystyle{ W(x) }[/math] - матрица Якоби(производных):

[math]\displaystyle{ W(x)=(\frac{\partial f_j}{\partial x_i})_{n,n}= \left\{\begin{matrix}(\frac{\partial f_1}{\partial x_1} \frac{\partial f_1}{\partial x_2} ... \frac{\partial f_1}{\partial x_n}) \\ (... ... ... ...) \\( \frac{\partial f_n}{\partial x_1} \frac{\partial f_n}{\partial x_2} ... \frac{\partial f_n}{\partial x_n} ) \end{matrix}\right. }[/math]

[math]\displaystyle{ \Delta X^{(k)}= \left\{\begin{matrix} \Delta (x_1)^{(k)} \\ \Delta (x_2)^{(k)} \\ ... \\ \Delta (x_n)^{(k)} \end{matrix}\right. }[/math]

[math]\displaystyle{ F(x) }[/math] - вектор-функция

[math]\displaystyle{ W(X^{(k)}) }[/math] - матрица Якоби, вычисленная для очередного приближения [math]\displaystyle{ F(X^{(k)}) }[/math] - вектор-функция, вычисленная для очередного приближения

Выразим вектор поправок [math]\displaystyle{ X^{(k)}=-W^{-1}(X^{(k)})*F(X^{(k)}) }[/math] :

[math]\displaystyle{ W^{-1} }[/math], где [math]\displaystyle{ W^{-1} }[/math]

Окончательная формула последовательных приближений метода Ньютона решения СНУ в матричной форме имеет вид:

[math]\displaystyle{ X^{(k+1)}=X^{(k)}=W^{-1}(X^{(k)})*F(X^{(k)}) }[/math]

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

Основная вычислительная нагрузка приходится на

1) Решение СЛАУ: [math]\displaystyle{ F(X^{(k)})=\frac{\partial F(x^{(k)})}{\partial x} \Delta x^{(k)} }[/math]

2)Численное вычисление Якобиана(если производные не даны): [math]\displaystyle{ \frac{\partial F(x^{(k)})}{\partial x} }[/math]

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

Данный алгоритм использует два основных метода решения в каждой итерации, это нахождением матрицы Якоби и решение СЛАУ.

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

1)Задаётся размерность системы [math]\displaystyle{ n }[/math], требуемая точность, начальное приближённое решение [math]\displaystyle{ X=(x_i)_n }[/math]

2)Вычисляются элементы матрицы Якоби [math]\displaystyle{ W=\left( {\partial f_i\over \partial x_i} \right)_{n,n} }[/math]

3)Вычисляется обратная матрица [math]\displaystyle{ W^{-1} }[/math]

4)Вычисляются вектор функция [math]\displaystyle{ F=(f_i)_n, f_i=f_i(x_1, x_2,..., x_n), i=1,...,n }[/math]

5)Вычисляются вектор поправок [math]\displaystyle{ \Delta X=W_{-1}*F }[/math]

6)Уточняется решение [math]\displaystyle{ X_{n+1}=X_n+\Delta X }[/math]

7)Оценивается достигнутая точность [math]\displaystyle{ \delta = \max_{i=1,n} \Delta x_i^k }[/math]

8)Проверяется условие завершения итерационного процесса [math]\displaystyle{ \delta\lt =\varepsilon }[/math]

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

Сперва стоит отметить, что итоговая сложность алгоритма зависит от того, насколько быстро он будет сходиться. Это в свою очередь зависит от заданного начального приближения [math]\displaystyle{ x^0 }[/math] и от условия остановки алгоритма [math]\displaystyle{ \varepsilon }[/math]. Можно однако вычислить сложность одного шага итерации алгоритма.

Предполагаем, что у нас нет значений производных заданных функций

[math]\displaystyle{ f_1(.), f_2(.), ..., f_n(.) }[/math].

Тогда используя формулу центральной разности для производной :

[math]\displaystyle{ f_i'(x) = (f(x+h) - f(x-h))/2h }[/math],

мы находим приближённое значение производной в интересующей нас точке с точностью порядка [math]\displaystyle{ O(h^2) }[/math] за 5 операций.

Учитывая, что в Якобиане содержится [math]\displaystyle{ n^2 }[/math] элементов - производные каждой функции по каждой переменной, - то для нахождения Якобиана нам суммарно требуется

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

операций. Сложность вычислений обратной матрицы

[math]\displaystyle{ W^{-1} }[/math]

зависит от выбранного алгоритма решения полученной СЛАУ. Будем использовать метод Гаусса. В таком случае сложность составит

[math]\displaystyle{ O(n^3) }[/math].

Таким образом сложность вычислительного ядра алгоритма составляет

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

для одного шага алгоритма.

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

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

Сложность вычисление элементов матрицы Якоби - производных [math]\displaystyle{ \frac{\partial F(x^{(k)})}{\partial x} }[/math] в случае, если они не заданы будет - [math]\displaystyle{ O(n^2/p) }[/math].

Решение СЛАУ [math]\displaystyle{ F(X^{(k)}) = \frac{\partial F(x^{(k)})}{\partial x} \Delta x^{(k)} }[/math] одним из параллельных методов [2] имеет сложность [math]\displaystyle{ O(n^3/p) }[/math].

Пусть алгоритм имеет [math]\displaystyle{ N }[/math] итераций, тогда итоговая сложность будет: [math]\displaystyle{ O(N \cdot \frac{n^3}{p}) }[/math]

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

Входные данные: [math]\displaystyle{ n^2 + n }[/math] функций от [math]\displaystyle{ n }[/math] переменных, [math]\displaystyle{ n }[/math] мерный вектор, [math]\displaystyle{ x^0 }[/math] - начальное приближение, [math]\displaystyle{ \varepsilon }[/math] - точность решения.

Выходные данные: [math]\displaystyle{ n }[/math] вещественных чисел

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

Для данного метода трудно универсально оценить соотношение последовательной и параллельной сложностей алгоритма, поскольку это зависит непосредственно от вида задаваемых нелинейных функций, способа решения СЛАУ и выбора начального приближения.