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

Решение системы нелинейных уравнений методом Ньютона
Последовательный алгоритм
Последовательная сложность	[math]O(n^2 [/math] - одна итерация
Объём входных данных	[math]n * m + 1[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math][/math]
Ширина ярусно-параллельной формы	[math][/math]

Автор описания: Арутюнов А.В.

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

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

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

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

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

Уравнение касательной к функции f(x) в точке имеет вид:

[math]f^(x_j)=y-f(x_n)/(x-x_n)[/math]

В уравнении касательной положим y=0 и [math]x=x_n+1[/math].

Тогда алгоритм последовательных вычислений в методе Ньютона состоит в следующем:

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

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

Пусть необходимо найти решение системы:

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

где [math]f = (f_1, ..., f_n) : \mathbb{R}^n \rightarrow \mathbb{R}^n[/math] - непрерывно дифференцируемое отображение в окрестности решения.

При начальном приближении [math]x_0[/math] и сделанных предположениях [math](k+1)[/math]-я итерация метода будет выглядеть следующим образом: [math] x_{k+1} = x_k - [f'(x_k)]^{-1}f(x_k)Δ [/math]

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

[math]\sum_{n=0}^\infty\frac{x^n}{n!}[/math]

[math] x_{k+1} = x_k - [f'(x_k)]^{-1}f(x_k) [/math]

Пусть необходимо найти решение системы:

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

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

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

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

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

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

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

Значения поправок используются для оценки достигнутой точности решения. Если максимальная по абсолютной величине поправка меньше заданной точности ε, расчет завершается. Таким образом, условие окончания расчета:

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

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

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

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

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

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

[math]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]\Delta X^{(k)}= \left\{\begin{matrix} \Delta (x_1)^{(k)} \\ \Delta (x_2)^{(k)} \\ ... \\ \Delta (x_n)^{(k)} \end{matrix}\right. [/math]

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

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

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

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

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

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

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

Основная вычислительная нагрузка приходится на 1) Решение СЛАУ:[math]F(X^{(k)}) = \frac{\partial F(x^{(k)})}{\partial x} \Delta x^{(k)}[/math]

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

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

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

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

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

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

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

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

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

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

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

8)Проверяется условие завершения итерационного процесса δ<=ε

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

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

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

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

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

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

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

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