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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

$\left\{\begin{matrix}f_1(x_1, ..., x_n) = 0 \\ ... \\ f_n(x_1, ..., x_n) = 0, \end{matrix}\right.$

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

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

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

$\sum_{n=0}^\infty\frac{x^n}{n!}$

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

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

$\left\{\begin{matrix}f_1(x_1, x_2, ..., x_n) = 0 \\ ... \\ f_n(x_1, x_2, ..., x_n) = 0, \end{matrix}\right.$

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

$\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.$

что позволяет свести исходную задачу СНУ к

многократному решению системы линейных уравнений.

Пусть известно приближение $(x_i)^k$ решения системы нелинейных уравнений $(x_i)^*$.Введем в рассмотрение поправку $\Delta x_i$ как разницу между решением и его приближением: $\Delta x_i=(x_i)^* — (x_i)^k=\gt (x_i)^*=(x_i)^k + \Delta x_i, i=\overline(1,n)$

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

$\left\{\begin{matrix}f_1((x_1)^k +Δx_1, (x_2)^k+Δx_2, ..., (x_n)^k + \Delta x_n)=0, \\ f_2((x_1)^k +\Delta x_1, (x_2)^k+Δx_2, ..., (x_n)^k + \Delta x_n)=0, \\ ... \\ f_n((x_1)^k +\Delta x_1, (x_2)^k+Δx_2, ..., (x_n)^k + \Delta x_n)=0, \end{matrix}\right.$

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

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

$\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)$

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

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

$\delta =\min_{i=\overline(1,n)|\Delta {x_i}^k|}$

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

$\delta =\frac{1}{n}\sum_{i=1}\n|\Delta x_i|\lt \varepsilon$

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

$W(\Delta X^k)*X^k=-F(X^k)$

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

$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.$

$\Delta X^k= \left\{\begin{matrix} \Delta (x_1)^k \\ \Delta (x_2)^k \\ ... \\ \Delta (x_n)^k \end{matrix}\right.$

$F(x)$ - вектор-функция

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

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

$W^{-1}$, где $W^{-1}$

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

$X^{(k+1)}=X^k=W^{-1}(X^k)*F(X^k)$

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

Основная вычислительная нагрузка приходится на 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.3 Макроструктура алгоритма

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

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

1)Задаётся размерность системы n, требуемая точность , начальное приближённое решение $X=(x_i)_n$ 2)Вычисляются элементы матрицы Якоби $W=\left( {f_i\over x_i} \right)_{n,n}$ 3)Вычисляется обратная матрица $W^{-1}$ 4)Вычисляются вектор функция $F=(f_i)_n, f_i=f_i(x_1, x_2,..., x_n), i=1,...,n$ 5)Вычисляются вектор поправок $\Delta X=W_{-1}*F$ 6)Уточняется решение $X_{n+1}=X_n+\Delta X$ 7)Оценивается достигнутая точность $δ=\max_{i=1,n} \Delta x_i^k$ 8)Проверяется условие завершения итерационного процесса δ<=ε